2012-02-17 11 views
7

R, bazı katsayılar tekillikler yüzünden düştüğünde, vcovHC() kullanarak sağlam standart hataları nasıl hesaplayabilirim? Standart lm fonksiyonu, gerçekte tahmin edilen tüm katsayılar için normal standart hataların hesaplanması gibi görünmektedir, ancak vcovHC() bir hata atar: "Ekmekte hata.% *% Et.: Uygun olmayan argümanlar".R tekil değerlerle lm modeli için sağlam standart hataları (vcovHC) hesaplar

(Kullandığım gerçek veriler biraz daha karmaşıktır. Aslında, iki farklı sabit efekt kullanan bir modeldir ve basitçe kurtulamayacağım yerel tekilliklere koşuyorum. Birinci faktörün 150 seviyesine ulaştığım iki sabit etki için ikincisinin 142 seviyesi vardır ve verilerin on blokta toplanması sonucu ortaya çıkan toplam 9 tekillik vardır.)

Burada Benim çıktı:

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 

Hatayı yeniden oluşturmak için kesilmiş en küçük bir kod burada. Singularities ile

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) 
+0

Sizin kodunuz işe yarayacak gibi görünüyor. Tekillikten kaçınmak için aylardan birini (veya kesişme) modelden açıkça kaldırmak isteyebilirsiniz. –

+0

Maalesef benim amacım: Tekilliği kaldıramıyorum. Bu, gönderdiğim basit bir örnek veri kümesidir. Bu veri kümesinde, katılıyorum: basitçe Dec'i regresyondan çıkarabilir, böylece tekillikten kurtulabilir ve vcovHC() çalışabilir. Benim gerçek verilerimde, tekilliğin, birçok seviyeyle (sırasıyla 150 ve 142) iki sabit etkiden kaynaklandığı biraz daha karmaşıktır. Bu verilerdeki tekilliklerden kurtulmanın bir yolunu bulamadım. – Chris

+3

@Chris: Hala bu hatayı alıyor musunuz? Tekilliği indüklemek için 'Dec'' 'c (rep (0,240))' 'ı değiştirdikten sonra' vcovHC (model) 'çağrısına not ettiğiniz hata olmadan başarılı olur. Sandviç 2.2-9 için changelog: 'lm/mlm/glm modelleri ile aliased parametreleri doğru bir şekilde (sandviç/vcovHC vb. Hatalara yol açıyor) doğru şekilde ele alınmadı, düzeltildi. – jthetzel

cevap

0

Modelleri iyi asla ve onlar sabit olmalıdır. Sizin durumunuzda, 12 ay için 12 katsayı var, ama aynı zamanda küresel müdahale! Yani aslında sadece 12 gerçek parametrenin tahmin edilmesi için 13 katsayı var. Ne gerçekten istediğiniz küresel yakalamayı devre dışı bırakmaktır - Daha aylık özgü kesişmesine gibi bir şey olacak, böylece: Eğer vcovHC ile herhangi bir sorun olmamalıdır böylece

> 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 

Ardından, normal bir modeldir.