R特異点を持つlmモデルのロバスト標準誤差（vcovHC）を計算する

Question

Rでは、特異点のために係数が落ちたときにvcovHC（）を使ってロバスト標準誤差を計算するにはどうすればよいですか？標準のlm関数は、実際に見積もられたすべての係数に対して正常な標準誤差を計算するように見えますが、vcovHC（）は「パンのエラー％*％肉：適合しない引数」というエラーをスローします。

（私が使用している実際のデータはもう少し複雑ですが、実際には2つの異なる固定効果を使ったモデルですが、私は単純に取り除くことができない地元の特異点に走ります。 2番目の固定効果には150レベル、2番目の固定レベルには142レベル、10ブロックでデータが収集されたことに起因する9つの特異点があります）

ここに

Call: 
lm(formula = one ~ two + three + Jan + Feb + Mar + Apr + May + 
Jun + Jul + Aug + Sep + Oct + Nov + Dec, data = dat) 

Residuals: 
    Min  1Q Median  3Q  Max 
-130.12 -60.95 0.08 61.05 137.35 

Coefficients: (1 not defined because of singularities) 
       Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1169.74313 57.36807 20.390 <2e-16 *** 
two   -0.07963 0.06720 -1.185 0.237  
three   -0.04053 0.06686 -0.606 0.545  
Jan   8.10336 22.05552 0.367 0.714  
Feb   0.44025 22.11275 0.020 0.984  
Mar   19.65066 22.02454 0.892 0.373  
Apr   -13.19779 22.02886 -0.599 0.550  
May   15.39534 22.10445 0.696 0.487  
Jun   -12.50227 22.07013 -0.566 0.572  
Jul   -20.58648 22.06772 -0.933 0.352  
Aug   -0.72223 22.36923 -0.032 0.974  
Sep   12.42204 22.09296 0.562 0.574  
Oct   25.14836 22.04324 1.141 0.255  
Nov   18.13337 22.08717 0.821 0.413  
Dec     NA   NA  NA  NA  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 69.63 on 226 degrees of freedom 
Multiple R-squared: 0.04878, Adjusted R-squared: -0.005939 
F-statistic: 0.8914 on 13 and 226 DF, p-value: 0.5629 

> model$se <- vcovHC(model) 
Error in bread. %*% meat. : non-conformable arguments 

エラーを再現するために最小限のコードが抜粋されています。

library(sandwich) 
set.seed(101) 
dat<-data.frame(one=c(sample(1000:1239)), 
       two=c(sample(200:439)), 
       three=c(sample(600:839)), 
       Jan=c(rep(1,20),rep(0,220)), 
       Feb=c(rep(0,20),rep(1,20),rep(0,200)), 
       Mar=c(rep(0,40),rep(1,20),rep(0,180)), 
       Apr=c(rep(0,60),rep(1,20),rep(0,160)), 
       May=c(rep(0,80),rep(1,20),rep(0,140)), 
       Jun=c(rep(0,100),rep(1,20),rep(0,120)), 
       Jul=c(rep(0,120),rep(1,20),rep(0,100)), 
       Aug=c(rep(0,140),rep(1,20),rep(0,80)), 
       Sep=c(rep(0,160),rep(1,20),rep(0,60)), 
       Oct=c(rep(0,180),rep(1,20),rep(0,40)), 
       Nov=c(rep(0,200),rep(1,20),rep(0,20)), 
       Dec=c(rep(0,220),rep(1,20))) 
model <- lm(one ~ two + three + Jan + Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec, data=dat) 
summary(model) 
model$se <- vcovHC(model) 

Answer 1

特異点を持つモデルは決して良好ではなく、固定する必要があります。あなたのケースでは、12ヶ月間に12の係数がありますが、グローバルな切片もあります！したがって実際には13個の係数があり、実際のパラメータは12個しか推定できません。あなたが実際にしたいことは世界的なインターセプトを無効にすることである - ので、あなたはより多くの月に固有の切片のようなものがあります。そして、

> model <- lm(one ~ 0 + two + three + Jan + Feb + Mar + Apr + May + Jun + Jul + Aug + Sep + Oct + Nov + Dec, data=dat) 
> summary(model) 

Call: 
lm(formula = one ~ 0 + two + three + Jan + Feb + Mar + Apr + 
    May + Jun + Jul + Aug + Sep + Oct + Nov + Dec, data = dat) 

Residuals: 
    Min  1Q Median  3Q  Max 
-133.817 -55.636 3.329 56.768 126.772 

Coefficients: 
     Estimate Std. Error t value Pr(>|t|)  
two  -0.09670 0.06621 -1.460 0.146  
three 0.02446 0.06666 0.367 0.714  
Jan 1130.05812 52.79625 21.404 <2e-16 *** 
Feb 1121.32904 55.18864 20.318 <2e-16 *** 
Mar 1143.50310 53.59603 21.336 <2e-16 *** 
Apr 1143.95365 54.99724 20.800 <2e-16 *** 
May 1136.36429 53.38218 21.287 <2e-16 *** 
Jun 1129.86010 53.85865 20.978 <2e-16 *** 
Jul 1105.10045 54.94940 20.111 <2e-16 *** 
Aug 1147.47152 54.57201 21.027 <2e-16 *** 
Sep 1139.42205 53.58611 21.263 <2e-16 *** 
Oct 1117.75075 55.35703 20.192 <2e-16 *** 
Nov 1129.20208 53.54934 21.087 <2e-16 *** 
Dec 1149.55556 53.52499 21.477 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 69.81 on 226 degrees of freedom 
Multiple R-squared: 0.9964, Adjusted R-squared: 0.9961 
F-statistic: 4409 on 14 and 226 DF, p-value: < 2.2e-16 

をあなたがvcovHCで何か問題があってはならないので、それは通常のモデルです。