在使用带固定效应的plm函数时，在R中计算聚集标准误差

Question

我正在尝试在R的plm包中运行具有固定效应和model = 'within'的回归，同时具有聚类标准误。使用plm中的Cigar数据集，我正在运行：

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales + factor(state), model = 'within', data = Cigar)
coeftest(model, vcovHC(model, type = 'HC0', cluster = 'group'))

  Estimate Std. Error t value Pr(>|t|)    
sales  -1.21956    0.21136 -5.7701 9.84e-09

这与我使用Stata（将Cigar文件写为.dta）获得的结果略有不同：

use cigar

xtset state year

xtreg price sales, fe vce(cluster state)


price   Coef.   Std. Err.   t   P>t [95% Conf.  Interval]

sales   -1.219563   .2137726    -5.70   0.000   -1.650124   -.7890033

即标准误差和T统计量是不同的。我尝试使用不同的“类型”重新运行R代码，但没有一种方法能够给出与Stata相同的结果。我是否漏掉了什么？

- user5450835

---

如果您正在寻找在三元固定效应中的聚类标准误差，请查看我的相关答案。 - jay.sf

Answer 1

Stata使用特定的小样本校正方法，该方法已在plm 1.5中实现。

尝试一下：

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales + factor(state), model = 'within', data = Cigar)
coeftest(model, function(x) vcovHC(x, type = 'sss'))

这将产生：

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
sales  -1.2196     0.2137  -5.707 1.415e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这将给出相同的标准误估计，最多保留3位数字：

x <- coeftest(model, function(x) vcovHC(x, type = 'sss'))
x[ , "Std. Error"]
## [1] 0.2136951

- landroni

---

请参考以下链接，其中讨论了在Stata和R中进行logit回归的不同鲁棒标准误差的情况。 - Helix123

@Helix123 只是部分相关。虽然那个问题确实涉及到Stata小样本校正，但它的重点在于 sandwich::vcovHC 而不是 plm::vcovHC（这个问题）。 - landroni

Answer 2

Stata使用有限样本修正来减少由于群集数量有限而导致的误差向下偏差。这是方差-协方差矩阵上的乘法因子，$c=\frac{G}{G-1} \cdot \frac{N-1}{N-K}$，其中G是组数，N是观测数，K是参数数。我认为coeftest只使用$c'=\frac{N-1}{N-K}$，因为如果我将R的标准误按c的第一项的平方缩放，我得到的结果与Stata的标准误非常接近：

display 0.21136*(46/(46-1))^(.5)
.21369554

这是我如何在R中复制Stata的操作：

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales, model = 'within', data = Cigar)
G <- length(unique(Cigar$state))
c <- G/(G - 1)
coeftest(model,c * vcovHC(model, type = "HC1", cluster = "group"))

这将产生：

t test of coefficients:

       Estimate Std. Error  t value   Pr(>|t|)    
sales -1.219563   0.213773 -5.70496 1.4319e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

这与Stata的误差为0.2137726和t统计量为-5.70相符。

这段代码可能不是很理想，因为数据中的状态数量可能与回归中的状态数量不同，但我懒得去弄清楚如何获取正确的面板数量。

- dimitriy

---

@user5450835 现在问题已经解决，请选择此答案。 - dimitriy

相关：https://dev59.com/Z14c5IYBdhLWcg3w7-B4#27368468 - Helix123