运行 plm 固定效应模型并添加一个因子虚拟变量(树方式固定效应)是否可以?
Is it ok to run a plm fixed effect model and add a factor dummy variable (tree way fixed effects)?
运行“plm”固定效应模型并在 R 中添加一个因子虚拟变量是否可以,如下所示?
“时间”、“公司”和“国家”这三个因素都是单独的指标,我想把它们一起固定。
我发现下面的规范更适合我的情况,而不是通过组合“公司”和“国家/地区”来制作两个指数。
这是可接受的格式吗?
plm(y ~ lag(x1, 1) + x2 + x3 + x4 + x5 + factor(Country), data=DATA,
index=c("Firm","Time"), model="within")
可以添加额外的因素。我们可以通过计算 LSDV 模型来证明这一点。作为初步说明,您当然需要稳健的标准误差,通常聚集在 highest aggregate level,即本例中的国家/地区。
注:以下用到R >= 4.1
LSDV
fit1 <-
lm(y ~ d + x1 + x2 + x3 + x4 + factor(id) + factor(time) + factor(country),
dat)
lmtest::coeftest(
fit1, vcov.=sandwich::vcovCL(fit1, cluster=dat$country, type='HC0')) |>
{\(.) .[!grepl('\(|factor', rownames(.)), ]}()
# Estimate Std. Error t value Pr(>|t|)
# d 10.1398727 0.3181993 31.8664223 4.518874e-191
# x1 1.1217514 1.6509390 0.6794627 4.968995e-01
# x2 3.4913273 2.7782157 1.2566797 2.089718e-01
# x3 0.6257981 3.3162148 0.1887085 8.503346e-01
# x4 0.1942742 0.8998307 0.2159008 8.290804e-01
添加factor(country)
后,我们用plm::plm
得到的估计量与LSDV相同:
plm::plm
fit2 <- plm::plm(y ~ d + x1 + x2 + x3 + x4 + factor(country),
index=c('id', 'time'), model='within', effect='twoways', dat)
summary(fit2, vcov=plm::vcovHC(fit2, cluster='group', type='HC1'))$coe
# Estimate Std. Error t-value Pr(>|t|)
# d 10.1398727 0.3232850 31.3651179 5.836597e-186
# x1 1.1217514 1.9440165 0.5770277 5.639660e-01
# x2 3.4913273 3.2646905 1.0694206 2.849701e-01
# x3 0.6257981 3.1189939 0.2006410 8.409935e-01
# x4 0.1942742 0.9250759 0.2100089 8.336756e-01
但是,cluster='group'
将引用 "id"
而不是 "country"
,因此 标准错误是错误的 。目前似乎无法通过 plm
的附加因素进行聚类,至少我什么都不知道。
或者,您可以使用 lfe::felm
来避免相对于 LSDV 大大减少的计算时间:
lfe::felm
summary(lfe::felm(y ~ d + x1 + x2 + x3 + x4 | id + time + country | 0 | country,
dat))$coe
# Estimate Cluster s.e. t value Pr(>|t|)
# d 10.1398727 0.3184067 31.8456637 1.826374e-33
# x1 1.1217514 1.6520151 0.6790201 5.004554e-01
# x2 3.4913273 2.7800267 1.2558611 2.153737e-01
# x3 0.6257981 3.3183765 0.1885856 8.512296e-01
# x4 0.1942742 0.9004173 0.2157602 8.301083e-01
为了比较,这是 Stata 的计算结果,标准误差与 LSDV 和 lfe::felm
:
非常相似
斯塔塔
. reghdfe y d x1 x2 x3 x4, absorb (country time id) vce(cluster country)
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
d | 10.13987 .3185313 31.83 0.000 9.49907 10.78068
x1 | 1.121751 1.652662 0.68 0.501 -2.202975 4.446478
x2 | 3.491327 2.781115 1.26 0.216 -2.103554 9.086209
x3 | .6257981 3.319675 0.19 0.851 -6.052528 7.304124
x4 | .1942742 .9007698 0.22 0.830 -1.617841 2.006389
_cons | 14.26801 23.65769 0.60 0.549 -33.32511 61.86114
模拟面板数据:
n1 <- 20; t1 <- 4; n2 <- 48
dat <- expand.grid(id=1:n1, time=1:t1, country=1:n2)
set.seed(42)
dat <- within(dat, {
id <- as.vector(apply(matrix(1:(n1*n2), n1), 2, rep, t1))
d <- runif(nrow(dat), 70, 80)
x1 <- sample(0:1, nrow(dat), replace=TRUE)
x2 <- runif(nrow(dat))
x3 <- runif(nrow(dat))
x4 <- rnorm(nrow(dat))
y <-
10*d + ## treatment effect
as.vector(replicate(n2, rep(runif(n1, 2, 5), t1))) + ## id FE
rep(runif(n1, 10, 12), each=t1) + ## time FE
rep(runif(n2, 10, 12), each=n1*t1) + ## country FE
- .7*x1 + 1.3*x2 + 2.4*x3 +
.5 * x4 + rnorm(nrow(dat), 0, 50)
})
readstata13::save.dta13(dat, 'panel.dta') ## for Stata
运行“plm”固定效应模型并在 R 中添加一个因子虚拟变量是否可以,如下所示?
“时间”、“公司”和“国家”这三个因素都是单独的指标,我想把它们一起固定。
我发现下面的规范更适合我的情况,而不是通过组合“公司”和“国家/地区”来制作两个指数。
这是可接受的格式吗?
plm(y ~ lag(x1, 1) + x2 + x3 + x4 + x5 + factor(Country), data=DATA,
index=c("Firm","Time"), model="within")
可以添加额外的因素。我们可以通过计算 LSDV 模型来证明这一点。作为初步说明,您当然需要稳健的标准误差,通常聚集在 highest aggregate level,即本例中的国家/地区。
注:以下用到R >= 4.1
LSDV
fit1 <-
lm(y ~ d + x1 + x2 + x3 + x4 + factor(id) + factor(time) + factor(country),
dat)
lmtest::coeftest(
fit1, vcov.=sandwich::vcovCL(fit1, cluster=dat$country, type='HC0')) |>
{\(.) .[!grepl('\(|factor', rownames(.)), ]}()
# Estimate Std. Error t value Pr(>|t|)
# d 10.1398727 0.3181993 31.8664223 4.518874e-191
# x1 1.1217514 1.6509390 0.6794627 4.968995e-01
# x2 3.4913273 2.7782157 1.2566797 2.089718e-01
# x3 0.6257981 3.3162148 0.1887085 8.503346e-01
# x4 0.1942742 0.8998307 0.2159008 8.290804e-01
添加factor(country)
后,我们用plm::plm
得到的估计量与LSDV相同:
plm::plm
fit2 <- plm::plm(y ~ d + x1 + x2 + x3 + x4 + factor(country),
index=c('id', 'time'), model='within', effect='twoways', dat)
summary(fit2, vcov=plm::vcovHC(fit2, cluster='group', type='HC1'))$coe
# Estimate Std. Error t-value Pr(>|t|)
# d 10.1398727 0.3232850 31.3651179 5.836597e-186
# x1 1.1217514 1.9440165 0.5770277 5.639660e-01
# x2 3.4913273 3.2646905 1.0694206 2.849701e-01
# x3 0.6257981 3.1189939 0.2006410 8.409935e-01
# x4 0.1942742 0.9250759 0.2100089 8.336756e-01
但是,cluster='group'
将引用 "id"
而不是 "country"
,因此 标准错误是错误的 。目前似乎无法通过 plm
的附加因素进行聚类,至少我什么都不知道。
或者,您可以使用 lfe::felm
来避免相对于 LSDV 大大减少的计算时间:
lfe::felm
summary(lfe::felm(y ~ d + x1 + x2 + x3 + x4 | id + time + country | 0 | country,
dat))$coe
# Estimate Cluster s.e. t value Pr(>|t|)
# d 10.1398727 0.3184067 31.8456637 1.826374e-33
# x1 1.1217514 1.6520151 0.6790201 5.004554e-01
# x2 3.4913273 2.7800267 1.2558611 2.153737e-01
# x3 0.6257981 3.3183765 0.1885856 8.512296e-01
# x4 0.1942742 0.9004173 0.2157602 8.301083e-01
为了比较,这是 Stata 的计算结果,标准误差与 LSDV 和 lfe::felm
:
斯塔塔
. reghdfe y d x1 x2 x3 x4, absorb (country time id) vce(cluster country)
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
d | 10.13987 .3185313 31.83 0.000 9.49907 10.78068
x1 | 1.121751 1.652662 0.68 0.501 -2.202975 4.446478
x2 | 3.491327 2.781115 1.26 0.216 -2.103554 9.086209
x3 | .6257981 3.319675 0.19 0.851 -6.052528 7.304124
x4 | .1942742 .9007698 0.22 0.830 -1.617841 2.006389
_cons | 14.26801 23.65769 0.60 0.549 -33.32511 61.86114
模拟面板数据:
n1 <- 20; t1 <- 4; n2 <- 48
dat <- expand.grid(id=1:n1, time=1:t1, country=1:n2)
set.seed(42)
dat <- within(dat, {
id <- as.vector(apply(matrix(1:(n1*n2), n1), 2, rep, t1))
d <- runif(nrow(dat), 70, 80)
x1 <- sample(0:1, nrow(dat), replace=TRUE)
x2 <- runif(nrow(dat))
x3 <- runif(nrow(dat))
x4 <- rnorm(nrow(dat))
y <-
10*d + ## treatment effect
as.vector(replicate(n2, rep(runif(n1, 2, 5), t1))) + ## id FE
rep(runif(n1, 10, 12), each=t1) + ## time FE
rep(runif(n2, 10, 12), each=n1*t1) + ## country FE
- .7*x1 + 1.3*x2 + 2.4*x3 +
.5 * x4 + rnorm(nrow(dat), 0, 50)
})
readstata13::save.dta13(dat, 'panel.dta') ## for Stata