2012-02-17 9 views
7

आर में, मैं vcovHC() का उपयोग करके मजबूत मानक त्रुटियों की गणना कैसे कर सकता हूं जब कुछ गुणांक एकवचन के कारण गिराए जाते हैं? मानक एलएम फ़ंक्शन वास्तव में अनुमानित सभी गुणांक के लिए सामान्य मानक त्रुटियों की गणना करने लगता है, लेकिन vcovHC() एक त्रुटि फेंकता है: "रोटी में त्रुटि।% *% मांस।: गैर-अनुरूपनीय तर्क"।आर एकवचन के साथ एलएम मॉडल के लिए मजबूत मानक त्रुटियों (vcovHC) की गणना

(वास्तविक डेटा जो मैं उपयोग कर रहा हूं वह थोड़ा अधिक जटिल है। असल में, यह दो अलग-अलग निश्चित प्रभावों का उपयोग कर एक मॉडल है और मैं स्थानीय एकवचन में भाग लेता हूं जिसे मैं आसानी से छुटकारा नहीं पा सकता हूं। कम से कम मुझे नहीं पता दो फिक्स्ड इफेक्ट्स के लिए मैं पहले कारक का उपयोग कर 150 स्तरों का उपयोग कर रहा हूं, दूसरे में 142 स्तर हैं और कुल 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) 
+0

आपका कोड काम करने लगता है। एकवचन से बचने के लिए, आप मॉडल से एक महीने (या अवरोध) को स्पष्ट रूप से निकालना चाहते हैं। –

+0

दुर्भाग्यवश मेरा मुद्दा है: मैं एकवचन को हटा नहीं सकता। यह सिर्फ एक साधारण उदाहरण डेटा सेट है जिसे मैंने पोस्ट किया था। उस डेटासेट में, मैं सहमत हूं: आप आसानी से रिग्रेशन से दिसंबर को हटा सकते हैं, इस प्रकार एकवचन से छुटकारा पा सकते हैं और फिर vcovHC() काम करेगा। मेरे वास्तविक डेटा में यह थोड़ा और जटिल है जहां एकवचनता कई निश्चित प्रभावों से क्रमशः दो स्तरों (150 और 142) के साथ होती है। मुझे उस डेटा में एकवचन से छुटकारा पाने का कोई तरीका नहीं मिला है। – Chris

+3

@ क्रिस: क्या आपको अभी भी यह त्रुटि मिलती है? एकवचन को प्रेरित करने के लिए 'दिसंबर' से 'सी (प्रतिनिधि (0,240)) को बदलने के बाद, 'vcovHC (मॉडल)' के लिए एक कॉल आपके द्वारा नोट की गई त्रुटि के बिना सफल हो जाती है। सैंडविच 2.2-9 के लिए चेंजलॉग में: एलआईएड/एमएलएम/एलआईएम मॉडल एलियाज्ड पैरामीटर के साथ सही ढंग से संभाला नहीं गया था (अब सैंडविच/वीसीओवीएचसी इत्यादि में त्रुटियों का कारण बनता है), अब तय किया गया है। शायद यह तय है? – jthetzel

उत्तर

0

मॉडल अच्छा कभी नहीं होते हैं और वे किया जाना चाहिए। आपके मामले में, आपके पास 12 महीने के लिए 12 गुणांक हैं, लेकिन वैश्विक अवरोध भी है! तो आपके पास अनुमानित होने के लिए केवल 12 वास्तविक पैरामीटर के लिए वास्तव में 13 गुणांक हैं। क्या आप वास्तव में चाहते हैं वैश्विक अवरोधन निष्क्रिय करने के लिए है - तो आप अधिक महीने के विशिष्ट अवरोधन की तरह कुछ करना होगा:

> 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 के साथ कोई समस्या नहीं होनी चाहिए।

संबंधित मुद्दे