为什么我的带有贬值数据 ( lm() ) 的 fe 不能从 plm()、lfe() 和 lsdv 中再现系数?

Why does my fe with demeaned data ( lm() ) not reproduce coefficients from plm(), lfe() and lsdv?

我正在尝试使用不同的包和技术重现面板数据的固定效应系数:(1) plm(), (2) lfe(), (3) dummy-lsdv with lm(),以及 (4) 贬低-fe lm().

我的数据集包含 1581 个观察值和 13 个变量。它是来自 527 名受访者(var = respondent)的 3 波(var = wave)调查数据。我有一个 DV (y) 和 10 个 IV(x1 到 x10)。

数据集是这样的:

  respondent  wave      y    x1    x2    x3    x4    x5    x6    x7    x8    x9   x10
       <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1          1     1  1         1     2    NA    NA     2   1     1.5  NA      NA     2
2          1     2 NA         2    NA     0     0     0   1     4     4       1     3
3          1     3 NA         4     5    NA    NA    NA  NA     8    NA      NA     1
4          2     1  0.931     3     3     2     2     2   4     7.5   7.5    NA     3
5          2     2  0.986     4    NA    NA     2     2   4.5   6.5   5       3     4
6          2     3  0.986     4     3     2     2     2   3     3     3       2     3

我的问题: 当我使用 (1) plm()、(2) lfe() 和 (3) dummy-lsdv 执行固定效应回归时使用 lm(),模型总是 return 相同的系数。但是,当我使用 (4) 贬值数据和 lm() 包执行固定效应回归时,我得到了不同的系数。这让我很困惑,我想知道:为什么?

这是我的代码:

1. plm():

输入:

library(plm)
model_plm <- plm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
             data = dataset,
             index=c("respondent","wave"),
             model = "within",
             effect = 'individual')
summary(model_plm)

输出:

Unbalanced Panel: n = 228, T = 1-2, N = 316

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.3240866 -0.0048416  0.0000000  0.0048416  0.3240866 

Coefficients:
      Estimate Std. Error t-value Pr(>|t|)  
x1  -0.0216484  0.0167614 -1.2916  0.20032  
x2   0.0178114  0.0141219  1.2613  0.21097  
x3  -0.0145262  0.0103954 -1.3974  0.16627  
x4  -0.0061660  0.0133069 -0.4634  0.64439  
x5   0.0174401  0.0144256  1.2090  0.23032  
x6  -0.0053556  0.0067210 -0.7968  0.42796  
x7   0.0065517  0.0097627  0.6711  0.50415  
x8  -0.0151375  0.0081992 -1.8462  0.06865 .
x9   0.0235351  0.0092612  2.5412  0.01303 *
x10  0.0235181  0.0228927  1.0273  0.30745  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2。 lfe():

输入:

library(lfe)
model_lfe <- felm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 | respondent, data = dataset)
summary(model_lfe)

输出:

Coefficients:
     Estimate Std. Error t value Pr(>|t|)  
x1  -0.021648   0.016761  -1.292   0.2003  
x2   0.017811   0.014122   1.261   0.2110  
x3  -0.014526   0.010395  -1.397   0.1663  
x4  -0.006166   0.013307  -0.463   0.6444  
x5   0.017440   0.014426   1.209   0.2303  
x6  -0.005356   0.006721  -0.797   0.4280  
x7   0.006552   0.009763   0.671   0.5041  
x8  -0.015138   0.008199  -1.846   0.0687 .
x9   0.023535   0.009261   2.541   0.0130 *
x10  0.023518   0.022893   1.027   0.3074  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

3。 LSDV lm():

输入:

model_lsdv <- lm(y ~ as_factor(respondent) + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, data = dataset)
options(max.print=2000)
summary(model_lsdv)

输出:

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.9499746  0.1505806   6.309 1.57e-08 ***
[...]    
x1                    -0.0216484  0.0167614  -1.292  0.20032    
x2                     0.0178114  0.0141219   1.261  0.21097    
x3                    -0.0145262  0.0103954  -1.397  0.16627    
x4                    -0.0061660  0.0133069  -0.463  0.64439    
x5                     0.0174401  0.0144256   1.209  0.23032    
x6                    -0.0053556  0.0067210  -0.797  0.42796    
x7                     0.0065517  0.0097627   0.671  0.50415    
x8                    -0.0151375  0.0081992  -1.846  0.06865 .  
x9                     0.0235351  0.0092612   2.541  0.01303 *  
x10                    0.0235181  0.0228927   1.027  0.30745    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.使用 lm():

贬低 FE

输入:

dataset_demeaned <- with(dataset, data.frame(respondent = respondent,
                                             wave = wave,
                                             y = y - ave(y, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x1 = x1 - ave(x1, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x2 = x2 - ave(x2, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x3 = x3 - ave(x3, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x4 = x4 - ave(x4, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x5 = x5 - ave(x5, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x6 = x6 - ave(x6, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x7 = x7 - ave(x7, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x8 = x8 - ave(x8, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x9 = x9 - ave(x9, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x10 = x10 - ave(x10, respondent, FUN=function(x) mean(x, na.rm=T))
                                            )
                        )

model_dmd <- lm(y ~ 0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, data = dataset_demeaned)
summary(model_dmd)

输出:

Coefficients:
     Estimate Std. Error t value Pr(>|t|)   
x1  -0.006223   0.008220  -0.757  0.44957   
x2   0.013181   0.007880   1.673  0.09543 . 
x3  -0.012807   0.005484  -2.335  0.02018 * 
x4  -0.006431   0.006311  -1.019  0.30900   
x5   0.015455   0.005941   2.602  0.00973 **
x6  -0.001429   0.003402  -0.420  0.67483   
x7   0.004362   0.004698   0.929  0.35387   
x8  -0.009336   0.004366  -2.139  0.03326 * 
x9   0.015731   0.005267   2.987  0.00305 **
x10  0.007631   0.010922   0.699  0.48529   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

更多信息:
我已经执行了这些检查:

所以我最好的猜测是问题不是来自贬低的过程,而是来自lm()函数。也许面板是 unbalanced 的事实在这里发挥了作用?

非常感谢您的任何建议和解释!


解决方案:

感谢@G.Grothendieck,我可以post解决这里问题。 (4) Demeaned FE with lm() 的正确代码应该是:

输入:

# Delete all rows with NAs
dataset <- na.omit(dataset)

# Demean the rows that are left behind
dataset_demeaned <- with(dataset, data.frame(respondent = respondent,
                                             wave = wave,
                                             y = y - ave(y, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x1 = x1 - ave(x1, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x2 = x2 - ave(x2, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x3 = x3 - ave(x3, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x4 = x4 - ave(x4, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x5 = x5 - ave(x5, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x6 = x6 - ave(x6, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x7 = x7 - ave(x7, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x8 = x8 - ave(x8, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x9 = x9 - ave(x9, respondent, FUN=function(x) mean(x, na.rm=T)),
                                             x10 = x10 - ave(x10, respondent, FUN=function(x) mean(x, na.rm=T))
                                             )
                         )

model_dmd <- lm(y ~ 0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, data = dataset_demeaned)
summary(model_dmd)

输出:

Coefficients:
     Estimate Std. Error t value Pr(>|t|)    
x1  -0.021648   0.008462  -2.558 0.011004 *  
x2   0.017811   0.007130   2.498 0.013009 *  
x3  -0.014526   0.005248  -2.768 0.005989 ** 
x4  -0.006166   0.006718  -0.918 0.359452    
x5   0.017440   0.007283   2.395 0.017240 *  
x6  -0.005356   0.003393  -1.578 0.115530    
x7   0.006552   0.004929   1.329 0.184768    
x8  -0.015138   0.004140  -3.657 0.000301 ***
x9   0.023535   0.004676   5.033 8.24e-07 ***
x10  0.023518   0.011558   2.035 0.042734 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

通过分别贬低每一列,这是在不一致地处理 NA。必须使用或不使用每一行。不能将一行用于一个变量而不用于另一个变量。