2012-02-17 48 views
7

Trong R, làm thế nào tôi có thể tính các lỗi tiêu chuẩn mạnh mẽ bằng vcovHC() khi một số hệ số bị giảm do số ít? Hàm lm chuẩn dường như làm tốt việc tính toán các lỗi chuẩn thông thường cho tất cả các hệ số thực sự được ước tính, nhưng vcovHC() ném một lỗi: "Lỗi trong bánh mì.% *% Thịt.: Đối số không phù hợp".R tính toán các lỗi tiêu chuẩn mạnh (vcovHC) cho mô hình lm với số ít

(Dữ liệu thực tế tôi đang sử dụng phức tạp hơn một chút. Trên thực tế, nó là một mô hình sử dụng hai hiệu ứng cố định khác nhau và tôi chạy vào những điểm kỳ dị địa phương mà tôi không thể loại bỏ. Vì hai hiệu ứng cố định tôi đang sử dụng yếu tố đầu tiên có 150 cấp độ, thứ hai có 142 cấp độ và có tổng cộng 9 điểm riêng lẻ do dữ liệu được thu thập trong mười khối.)

Ở đây là đầu ra của tôi:

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 

Đây là một mã tối thiểu được nén để tạo lại lỗi.

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

Mã của bạn dường như hoạt động. Bạn có thể muốn loại bỏ một cách rõ ràng một trong các tháng (hoặc chặn) khỏi mô hình, để tránh sự kỳ dị. –

+0

Thật không may là quan điểm của tôi: Tôi không thể loại bỏ sự kỳ dị. Đây chỉ là một tập dữ liệu ví dụ đơn giản mà tôi đã đăng. Trong tập dữ liệu đó, tôi đồng ý: bạn có thể chỉ cần xóa tháng 12 khỏi hồi quy, do đó loại bỏ sự kỳ dị và sau đó vcovHC() sẽ hoạt động. Trong dữ liệu thực tế của tôi, nó phức tạp hơn một chút khi điểm kỳ dị xuất phát từ hai hiệu ứng cố định với nhiều cấp độ (tương ứng là 150 và 142). Tôi đã không tìm thấy cách để loại bỏ các điểm kỳ dị trong dữ liệu đó. – Chris

+3

@Chris: Bạn vẫn nhận được lỗi này? Sau khi thay đổi 'Dec' thành' c (rep (0,240)) 'để tạo ra một điểm kỳ dị, một lời gọi tới' vcovHC (model) 'thành công mà không có lỗi bạn lưu ý. Trong các thay đổi cho bánh sandwich 2.2-9: 'lm/mlm/glm mô hình với các thông số bí danh đã không được xử lý một cách chính xác (dẫn đến lỗi trong bánh sandwich/vcovHC vv), cố định ngay bây giờ.' Có lẽ đó cố định nó? – jthetzel

Trả lời

0

Mô hình có điểm kỳ dị không bao giờ tốt và chúng phải được sửa. Trong trường hợp của bạn, bạn có 12 hệ số trong 12 tháng, nhưng cũng chặn toàn cầu! Vì vậy, bạn có 13 hệ số cho chỉ 12 tham số thực được ước tính. Những gì bạn thực sự muốn là không cho phép đánh chặn toàn cầu - vì vậy bạn sẽ có một cái gì đó giống như tháng cụ thể đánh chặn:

> 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 

Sau đó, nó là một mô hình bình thường, do đó bạn không nên có bất kỳ vấn đề với vcovHC.

Các vấn đề liên quan