零膨脹負二項模型的集羣標準誤差

Question

我想計算零膨脹負二項模型的集羣標準誤差。默認情況下，zeroinfl（從pscl包）返回使用由optim返回的Hessian矩陣，例如衍生標準誤差：

library(pscl) 
data("bioChemists", package = "pscl") 
dim(bioChemists) 
head(bioChemists) 
## default start values 
fm1 <- zeroinfl(art ~ ., data = bioChemists, dist = "negbin")) 
summary(fm1) 

有一種使用非對稱的/對稱距離矩陣觀測之間或使用中的一個的方式變量（例如玩具數據集中的kid5）來計算聚集標準誤差？

我發現this從an answer at stackexchange，但我不知道如何/是否可以使用零膨脹模型。在Stata的rzinb中的等價物可能是cluster clustvar，其在vce下：http://www.stata.com/manuals13/rzinb.pdf。

任何想法？

Answer 1

R-Forge上的sandwich軟件包的開發版本已經擴展爲允許面向對象的集羣協方差計算。這也支持零膨脹迴歸模型。您可以通過R-Forge安裝Devel版本：

install.packages("sandwich", repos = "http://R-Forge.R-project.org") 

然後加載所有必需的軟件包。該lmtest包用於coeftest()功能到其中的協方差矩陣的估計可以插入。

library("pscl") 
library("sandwich") 
library("lmtest") 

您使用的插圖模型如下。

data("bioChemists", package = "pscl") 
fm1 <- zeroinfl(art ~ ., data = bioChemists, dist = "negbin") 

的coeftest()功能默認返回相同的邊際瓦爾德測試，summary()。

coeftest(fm1) 
## t test of coefficients: 
## 
##      Estimate Std. Error t value Pr(>|t|)  
## count_(Intercept) 0.41674653 0.14359655 2.9022 0.003796 ** 
## count_femWomen -0.19550683 0.07559256 -2.5863 0.009856 ** 
## count_marMarried 0.09758263 0.08445195 1.1555 0.248199  
## count_kid5  -0.15173246 0.05420606 -2.7992 0.005233 ** 
## count_phd   -0.00070013 0.03626966 -0.0193 0.984603  
## count_ment   0.02478620 0.00349267 7.0966 2.587e-12 *** 
## zero_(Intercept) -0.19168829 1.32281889 -0.1449 0.884815  
## zero_femWomen  0.63593320 0.84891762 0.7491 0.453986  
## zero_marMarried -1.49946849 0.93867060 -1.5974 0.110518  
## zero_kid5   0.62842720 0.44278263 1.4193 0.156166  
## zero_phd   -0.03771474 0.30800817 -0.1224 0.902572  
## zero_ment   -0.88229322 0.31622813 -2.7901 0.005381 ** 
## --- 
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

這可以被容易地擴展使用vcovCL()函數採用一個羣集的協方差矩陣的估計。這裏，kid5變量按照您的建議使用。 （注意，如果別人正在讀這個：kid5的用法只是爲了展示「工作」，但在這個應用程序中並沒有真正意義。）

coeftest(fm1, vcov = vcovCL(fm1, cluster = bioChemists$kid5)) 
## t test of coefficients: 
## 
##      Estimate Std. Error t value Pr(>|t|)  
## count_(Intercept) 0.41674653 0.17009748 2.4500 0.01447 * 
## count_femWomen -0.19550683 0.01701325 -11.4914 < 2.2e-16 *** 
## count_marMarried 0.09758263 0.02401883 4.0628 5.272e-05 *** 
## count_kid5  -0.15173246 0.03612916 -4.1997 2.938e-05 *** 
## count_phd   -0.00070013 0.04852615 -0.0144 0.98849  
## count_ment   0.02478620 0.00263208 9.4170 < 2.2e-16 *** 
## zero_(Intercept) -0.19168829 0.51865043 -0.3696 0.71177  
## zero_femWomen  0.63593320 0.87775846 0.7245 0.46895  
## zero_marMarried -1.49946849 1.03481783 -1.4490 0.14768  
## zero_kid5   0.62842720 0.35073624 1.7917 0.07351 . 
## zero_phd   -0.03771474 0.13873870 -0.2718 0.78581  
## zero_ment   -0.88229322 0.07481264 -11.7934 < 2.2e-16 *** 
## --- 
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1