如何使用lfe::felm创建多变量和固定效果的回归？

Question

我有一个大的数据集，想要回归变量A，B，C，D和E上的变量X，还包括W，Y，Z年的固定效应。我也想使用变量C的自然日志(加一个)。

我该怎么做呢？

我的直觉是用猫

# install.packages("lfe")
library(lfe)
regress <- felm(formula= X ~ A, B, C, D, E + W + Y + Z)
regress

Answer 1

在我看来，你有像这样的假人的数据：

head(dat, 3)
#   id year2020 year2021 year2022         y         x1       x2
# 1  1        1        0        0 107.42623  32.903003 298.8692
# 2  2        1        0        0  32.90695 -13.552756 187.4316
# 3  3        1        0        0 123.78364   8.715082 507.0717

您将需要创建一个来自年份假人的因素变量，如下所示：

ycols <- c('year2020', 'year2021', 'year2022')
dat$year <- gsub('year', '', ycols)[apply(dat[ycols], 1, which.max)]
dat <- dat[setdiff(names(dat), ycols)]

head(dat, 3)
#   id         y         x1       x2 year
# 1  1 107.42623  32.903003 298.8692 2020
# 2  2  32.90695 -13.552756 187.4316 2020
# 3  3 123.78364   8.715082 507.0717 2020

然后以这种方式在公式中添加固定效果(请参阅help('felm'))，在这里您也可以直接使用log和任何加法。

library(lfe)

est1 <- felm(y ~ x1 + log(x2 + .001) | id + year, dat)
summary(est1)
# Call:
#   felm(formula = y ~ x1 + log(x2 + 0.001) | id + year, data = dat) 
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -104.680  -14.638   -0.788   13.477  150.973 
# 
# Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
# x1               1.03560    0.09927   10.43   <2e-16 ***
# log(x2 + 0.001) 71.01762    2.88983   24.57   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 33.47 on 196 degrees of freedom
# Multiple R-squared(full model): 0.8559   Adjusted R-squared: 0.7802 
# Multiple R-squared(proj model): 0.7894   Adjusted R-squared: 0.6787 
# F-statistic(full model): 11.3 on 103 and 196 DF, p-value: < 2.2e-16 
# F-statistic(proj model): 367.3 on 2 and 196 DF, p-value: < 2.2e-16 

您可以使用LSDV确认结果：

est2 <- lm(y ~ 0 + x1 + log(x2 + 0.001) + factor(id) + factor(year), dat)
summary(est2)$coe[1:2, ] |> signif(5)
#                  Estimate Std. Error t value   Pr(>|t|)
# x1                1.0356   0.099272  10.432 1.5098e-20
# log(x2 + 0.001)  71.0180   2.889800  24.575 9.0691e-62

数据：

nid <- 100; nyr <- 3
set.seed(42)
dat <- expand.grid(id=factor(seq_len(nid)), year=factor(2019+seq_len(nyr)))
dat <- within(dat, {
  x1 <- rnorm(nid*nyr, 0, 24)
  x2 <- rgamma(nid*nyr, scale=200, shape=2)
  y <- x1 + .25*x2 + rnorm(nlevels(id)) + rnorm(nlevels(year)) +
    rnorm(nid*nyr, 0, 12)
})