estimatr::lm_robust() 估计与 lm() 不同的系数

Question

使用

stats::lm()

运行面板模型（5475 个单位，13 年）后，我一直在尝试使用

estimatr::lm_robust()

估计异方差稳健标准误差。事实证明，每个函数中估计的一些系数（非常）不同 - 这不应该发生（只有 se、t.test 和 p 值应该改变）。

# df <- load(url("https://github.com/victorrssx/victorrssx/blob/main/df.RData?raw=true"))
df <- readRDS(url("https://github.com/victorrssx/victorrssx/raw/main/df.RDS")) %>% 
      dplyr::mutate(ano = as.factor(ano))


lm_mod <- stats::lm(tx_hom_tot ~ pos2018 + res_ou:pos2018 + ti:pos2018 + res_ou:ti:pos2018 + ano, data = df)
summary(lm_mod)

    Coefficients: (1 not defined because of singularities)
                      Estimate Std. Error t value             Pr(>|t|)    
    (Intercept)         26.827      0.382   70.25 < 0.0000000000000002 ***
    pos2018              2.927      0.538    5.44     0.00000005226909 ***
    ano2011              0.220      0.540    0.41              0.68385    
    ano2012              1.182      0.536    2.20              0.02754 *  
    ano2013              0.892      0.537    1.66              0.09640 .  
    ano2014              1.960      0.534    3.67              0.00025 ***
    ano2015              2.985      0.533    5.60     0.00000002131005 ***
    ano2016              5.556      0.530   10.48 < 0.0000000000000002 ***
    ano2017              6.835      0.528   12.95 < 0.0000000000000002 ***
    ano2018              4.841      0.531    9.12 < 0.0000000000000002 ***
    ano2019             -1.633      0.520   -3.14              0.00169 ** 
    ano2020              0.280      0.518    0.54              0.58944    
    ano2021             -0.272      0.517   -0.53              0.59925    
    ano2022                 NA         NA      NA                   NA    
    pos2018:res_ou       0.862      0.643    1.34              0.18005    
    pos2018:ti           2.874      0.629    4.57     0.00000498294524 ***
    pos2018:res_ou:ti    9.824      1.365    7.20     0.00000000000063 ***
            
lm_mod <- estimatr::lm_robust(tx_hom_tot ~ pos2018 + res_ou:pos2018 + ti:pos2018 + res_ou:ti:pos2018 + ano, data = df, se_type = "HC1")
summary(lm_mod)

Standard error type:  HC1 

Coefficients:
                            Estimate             Std. Error     t value  Pr(>|t|)
(Intercept)                   26.930              10407.875  0.00258749 0.9979355
pos2018            2842202543499.796 324275566161917504.000  0.00000876 0.9999930
ano2011                        0.222                196.496  0.00112810 0.9990999
ano2012                        1.179                330.843  0.00356436 0.9971561
ano2013                        0.890                397.685  0.00223694 0.9982152
ano2014                        1.957                125.512  0.01558929 0.9875621
ano2015                        2.982                396.852  0.00751471 0.9940042
ano2016                        5.553                470.073  0.01181333 0.9905746
ano2017                        6.832                355.406  0.01922433 0.9846622
ano2018                        4.838                355.006  0.01362761 0.9891271
ano2019           -2842202543498.518 264769891052461504.000 -0.00001073 0.9999914
ano2020           -2842202543497.091 374441170828414656.000 -0.00000759 0.9999939
ano2021           -2842202543497.189 324275566161908160.000 -0.00000876 0.9999930
ano2022           -2842202543496.979 324275566161925760.000 -0.00000876 0.9999930
pos2018:res_ou                 0.863                  5.743  0.15017708 0.8806255
pos2018:ti                     2.874                  0.677  4.24755884 0.0000216
pos2018:res_ou:ti              9.852               2739.434  0.00359654 0.9971304

具体来说，问题在

pos2018

、

ano2019

、

ano2020

、

ano2021

上严重发生，这些都是虚拟的。第一个在 2018 年之后的年份设置为 1，否则设置为 0；其余的是特定于年份的虚拟变量。即使设置

fixed_effects = ~ ano

也会产生奇怪的结果。

lm_mod <- estimatr::lm_robust(tx_hom_tot ~ pos2018 + res_ou:pos2018 + ti:pos2018 + res_ou:ti:pos2018, data = df, se_type = "HC1", fixed_effects = ~ ano)
summary(lm_mod) 

Standard error type:  HC1 

Coefficients: (1 not defined because the design matrix is rank deficient)
                  Estimate Std. Error t value   Pr(>|t|) CI Lower CI Upper    DF
pos2018                 NA         NA      NA         NA       NA       NA    NA
pos2018:res_ou       0.862      0.707    1.22 0.22251985   -0.523     2.25 50693
pos2018:ti           2.874      0.640    4.49 0.00000709    1.620     4.13 50693
pos2018:res_ou:ti    9.824      2.252    4.36 0.00001292    5.409    14.24 50693

使用

lmtest::coeftest(lm_mod, vcov. = plm::vcovHC(lm_mod, type = "HC1"))

时，似乎会产生正确的结果。那么为什么

lm_robust

会发生这种情况呢？难道是由于函数对数据中

NA

的解释？

Answer 1

您可以在 作者的 github 提交问题。这并不妨碍您同时使用经典方法计算 HC1 误差：

> lmtest::coeftest(lm_mod, vcov.=sandwich::vcovHC(lm_mod, type='HC1'))

t test of coefficients:

                  Estimate Std. Error t value  Pr(>|t|)    
(Intercept)       26.82729    0.33009 81.2731 < 2.2e-16 ***
pos2018            2.92704    0.51619  5.6705 1.432e-08 ***
ano2011            0.21990    0.46329  0.4747  0.635035    
ano2012            1.18210    0.47655  2.4805  0.013122 *  
ano2013            0.89249    0.47448  1.8810  0.059979 .  
ano2014            1.95953    0.47255  4.1467 3.378e-05 ***
ano2015            2.98510    0.48176  6.1962 5.827e-10 ***
ano2016            5.55602    0.50851 10.9262 < 2.2e-16 ***
ano2017            6.83532    0.52874 12.9275 < 2.2e-16 ***
ano2018            4.84076    0.52573  9.2076 < 2.2e-16 ***
ano2019           -1.63344    0.53712 -3.0411  0.002358 ** 
ano2020            0.27967    0.55194  0.5067  0.612367    
ano2021           -0.27170    0.56632 -0.4798  0.631404    
pos2018:res_ou     0.86248    0.70703  1.2199  0.222520    
pos2018:ti         2.87381    0.63982  4.4916 7.085e-06 ***
pos2018:res_ou:ti  9.82393    2.25226  4.3618 1.292e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

请注意，您应该先运行

summary(lm_mod)

，否则您将错过

coeftest

未显示的重要信息。