r - 如何在 R（新更新）中对纵向温度序列执行分段/样条回归？

问题描述

这里我有温度时间序列面板数据，我打算为它运行分段回归或三次样条回归。因此，首先我快速研究了分段回归概念及其在 R in 中的基本实现，SO初步了解了如何继续我的工作流程。在我的第一次尝试中，我尝试通过splines::ns在splines包中使用来运行样条回归，但我没有得到正确的条形图。对我来说，使用基线回归、分段回归或样条回归都可以。

这是我的面板数据规范的一般图片：下面显示的第一行是我的因变量，以自然对数项和自变量表示：平均温度、总降水量和 11 个温度箱和每个箱宽（AKA，箱的窗口) 为 3 摄氏度。(<-6, -6~-3,-3~0,...>21)。

可重现的例子：

以下是使用实际温度时间序列面板数据模拟的可重现数据：

set.seed(1) # make following random data same for everyone
dat <- data.frame(index=rep(c("dex111", "dex112", "dex113", "dex114", "dex115"), 
                          each=30),
                year=1980:2009,
                region= rep(c("Berlin", "Stuttgart", "Böblingen", 
                              "Wartburgkreis", "Eisenach"), each=30),
                ln_gdp_percapita=rep(sample.int(40, 30), 5), 
                ln_gva_agr_perworker=rep(sample.int(45, 30), 5),
                temperature=rep(sample.int(50, 30), 5), 
                precipitation=rep(sample.int(60, 30), 5), 
                bin1=rep(sample.int(32, 30), 5), 
                bin2=rep(sample.int(34, 30), 5), 
                bin3=rep(sample.int(36, 30), 5),
                bin4=rep(sample.int(38, 30), 5), 
                bin5=rep(sample.int(40, 30), 5), 
                bin6=rep(sample.int(42, 30), 5),
                bin7=rep(sample.int(44, 30), 5), 
                bin8=rep(sample.int(46, 30), 5), 
                bin9=rep(sample.int(48, 30), 5),
                bin10=rep(sample.int(50, 30), 5), 
                bin11=rep(sample.int(52, 30), 5))

请注意，除了极端温度值外，每个 bin 都有等分的温度区间，因此每个 bin 给出了落在各自温度区间内的天数。

更新 2：回归规范：

这是我的回归规范：

其中地区由索引i，年份由索引t。y_it是产出的量度 y_it∈ {ln GDP per capita, ln GVA per capita (by six sectors respectively)}，μ_i是一组地区固定效应，这些效应解释了地区之间未观察到的恒定差异。θ_t是一组灵活地解释共同趋势的年份固定效应。T_it^m is the number of days in the districti and yeart` 在第 m 个温度箱中具有一日平均温度。每个内部温度箱宽 3℃。当我对其运行样条回归时，我需要添加两种固定方式（按年份固定和按地区固定）。

新更新 1：

在这里，我想完全重新定义我的意图。最近我发现了非常有趣的 R 包，plm它适用于面板数据。这是我的新解决方案，使用plm效果很好：

library(plm)
pdf <- pdata.frame(dat, index = c("region", "year"))
model.b <- plm(ln_gdp_percapita ~ bin1+bin2+bin3+bin4+bin5+bin6+bin7+bin8+bin9+bin10+bin11, data = pdf, model = "pooling", effect = "twoways")

library(lmtest)    
coeftest(model.b)
res <- summary(model.b, cluster=c("c"))  ## add standard clustered error on it

新更新3 ：

summary(model.b, cluster=c("c"))$coefficients  # only render coefficient estimates table

新更新 2：我的输出：

    > coeftest(model.b)

t test of coefficients:

         Estimate  Std. Error t value  Pr(>|t|)    
bin1   1.7773e-04  4.8242e-04  0.3684 0.7125716    
bin2   2.4031e-03  4.3999e-04  5.4617 4.823e-08 ***
bin3   7.9238e-04  3.9733e-04  1.9943 0.0461478 *  
bin4  -2.0406e-05  3.7496e-04 -0.0544 0.9566001    
bin5   9.9911e-04  3.6386e-04  2.7459 0.0060451 ** 
bin6   6.0026e-05  3.4915e-04  0.1719 0.8635032    
bin7   2.5621e-04  3.0243e-04  0.8472 0.3969170    
bin8  -9.5919e-04  2.7136e-04 -3.5347 0.0004099 ***
bin9  -1.8195e-04  2.5906e-04 -0.7023 0.4824958    
bin10 -5.2064e-04  2.7006e-04 -1.9279 0.0538948 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

所需的散点图：

下面是我想要实现的散点图。它只是一个模拟散点图，灵感来自 NBER 工作论文第 32 页，标题为“温度对生产力和要素重新分配的影响：来自 50 万个中国制造工厂的证据 -此处提供非门控版本”，并且可以通过以下方式在整个文件中固定页面方向从命令行运行以下命令：
pdftk w23991.pdf cat 1-31 32-37east 38-40 41east 42-44 45east 46 output w23991-oriented.pdf

所需的散点图：

在该图中，黑点线是估计回归（基线或限制样条回归）系数，点蓝线是基于聚类标准误差的 95% 置信区间。

我刚刚联系了论文的作者，他们只是Excel用来获取那个情节。基本上，他们只是使用Estimate95% 置信区间数据的右侧和左侧来生成图。我知道那种情节Excel非常容易，但我有兴趣在R. 那可行吗？任何想法？

我想要一种更加程序化的方法来通过使用R而不是使用来呈现情节Excel。有什么聪明的举动吗？

标签： rregressionnon-linear-regressionplm

解决方案

前言：我完全不熟悉这个问题背后的统计数据。以下内容可能对开始使用有帮助ggplot2。让我知道你的想法。

set.seed(1) # make following random data same for everyone
dat <- data.frame(index=rep(c("dex111", "dex112", "dex113", "dex114", "dex115"), 
                              each=30),
                    year=1980:2009,
                    region= rep(c("Berlin", "Stuttgart", "Böblingen", 
                                  "Wartburgkreis", "Eisenach"), each=30),
                    ln_gdp_percapita=rep(sample.int(40, 30), 5), 
                    ln_gva_agr_perworker=rep(sample.int(45, 30), 5),
                    temperature=rep(sample.int(50, 30), 5), 
                    precipitation=rep(sample.int(60, 30), 5), 
                    bin1=rep(sample.int(32, 30), 5), 
                    bin2=rep(sample.int(34, 30), 5), 
                    bin3=rep(sample.int(36, 30), 5),
                    bin4=rep(sample.int(38, 30), 5), 
                    bin5=rep(sample.int(40, 30), 5), 
                    bin6=rep(sample.int(42, 30), 5),
                    bin7=rep(sample.int(44, 30), 5), 
                    bin8=rep(sample.int(46, 30), 5), 
                    bin9=rep(sample.int(48, 30), 5),
                    bin10=rep(sample.int(50, 30), 5), 
                    bin11=rep(sample.int(52, 30), 5))

library(plm)
pdf <- pdata.frame(dat, index=c("region", "year"))
model.b <- plm(ln_gdp_percapita ~ 
               bin1+bin2+bin3+bin4+bin5+bin6+bin7+bin8+bin9+bin10+bin11,
                   data=pdf, model="pooling", effect="twoways")
pdf$ln_gdp_percapita_predicted <- plm:::predict.plm(model.b, pdf)

library(ggplot2)
x <- ggplot(pdf, aes(y=ln_gdp_percapita_predicted, x=temperature))+
            geom_point()+
            geom_smooth(method=lm, formula=y~x, se=TRUE, level=.95)+ # see ?geom_smooth
            ylab("ln_gdp_percapita_predicted")+
            ggtitle("ln_gdp_percapita modeled as temperature")

ggsave("scatter_plot_2.png")
x

参考：R：使用 plm 和 pglm 绘制面板模型预测

更新：

制作一个情节res（??coefplot有关更多信息，请参阅）：

res <- plm:::summary.plm(model.b, cluster=c("c"))

library(coefplot)
coefplot::coefplot(res)
ggsave("model.b.coefplot.png")

r - 如何在 R（新更新）中对纵向温度序列执行分段/样条回归？

问题描述

解决方案

推荐阅读