Lỗi chuẩn khác nhau của hồi quy logit trong Stata và R

Question

Tôi đang cố gắng sao chép hồi quy logit từ Stata sang R. Trong Stata tôi sử dụng tùy chọn "mạnh" để có lỗi tiêu chuẩn mạnh mẽ (lỗi tiêu chuẩn thống nhất dị thường). Tôi có thể sao chép chính xác các hệ số tương tự từ Stata, nhưng tôi không thể có cùng một lỗi tiêu chuẩn mạnh với gói "sandwich".

Tôi đã thử một số ví dụ về hồi quy tuyến tính OLS; có vẻ như các bộ ước lượng bánh sandwich của R và Stata cho tôi cùng một lỗi tiêu chuẩn mạnh mẽ cho OLS. Có ai biết làm thế nào Stata tính toán các ước tính bánh sandwich cho hồi quy phi tuyến tính, trong trường hợp của tôi hồi quy logit?

Cảm ơn bạn!

Codes đính kèm: trong R:

library(sandwich) 
library(lmtest)  
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")  
mydata$rank<-factor(mydata$rank)  
myfit<-glm(admit~gre+gpa+rank,data=mydata,family=binomial(link="logit"))  
summary(myfit)  
coeftest(myfit, vcov = sandwich)  
coeftest(myfit, vcov = vcovHC(myfit, "HC0"))  
coeftest(myfit, vcov = vcovHC(myfit))  
coeftest(myfit, vcov = vcovHC(myfit, "HC3"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC1"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC2"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC"))  
coeftest(myfit, vcov = vcovHC(myfit, "const"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC4"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC4m"))  
coeftest(myfit, vcov = vcovHC(myfit, "HC5"))  

Stata:

use http://www.ats.ucla.edu/stat/stata/dae/binary.dta, clear  
logit admit gre gpa i.rank, robust  

Answer 1

Giá trị mặc định cái gọi là "mạnh mẽ" sai số chuẩn trong Stata tương ứng với những gì sandwich() từ gói cùng tên tính toán. Sự khác biệt duy nhất là cách điều chỉnh mẫu hữu hạn được thực hiện. Trong các chức năng sandwich(...) không có điều chỉnh mẫu hữu hạn được thực hiện ở tất cả theo mặc định, tức là, bánh sandwich được chia cho 1/n trong đó n là số quan sát. Ngoài ra, sandwich(..., adjust = TRUE) có thể được sử dụng chia cho 1/(n - k) trong đó k là số hồi quy. Và Stata chia cho 1/(n - 1).

Tất nhiên, tất nhiên, điều này không khác biệt chút nào. Và ngoại trừ một vài trường hợp đặc biệt (ví dụ, hồi quy tuyến tính OLS) không có đối số cho 1/(n - k) hoặc 1/(n - 1) để hoạt động "đúng" trong các mẫu hữu hạn (ví dụ: không thiên vị). Ít nhất là không theo sự hiểu biết tốt nhất của tôi.

Vì vậy, để có được kết quả tương tự như trong Stata bạn có thể làm làm:

sandwich1 <- function(object, ...) sandwich(object) * nobs(object)/(nobs(object) - 1) 
coeftest(myfit, vcov = sandwich1) 

này mang lại

z test of coefficients: 

       Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.9899791 1.1380890 -3.5059 0.0004551 *** 
gre   0.0022644 0.0011027 2.0536 0.0400192 * 
gpa   0.8040375 0.3451359 2.3296 0.0198259 * 
rank2  -0.6754429 0.3144686 -2.1479 0.0317228 * 
rank3  -1.3402039 0.3445257 -3.8900 0.0001002 *** 
rank4  -1.5514637 0.4160544 -3.7290 0.0001922 *** 
--- 
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Và chỉ cho các hồ sơ: Trong trường hợp phản ứng nhị phân, các "mạnh mẽ" lỗi tiêu chuẩn không mạnh mẽ chống lại bất cứ điều gì. Miễn là mô hình được chỉ định chính xác, chúng nhất quán và không sao để sử dụng chúng nhưng chúng không bảo vệ chống lại bất kỳ lỗi chính xác nào trong mô hình. Bởi vì giả định cơ bản cho các lỗi chuẩn của sandwich là công thức mô hình (hoặc chính xác hơn là hàm số điểm tương ứng) được chỉ định chính xác trong khi phần còn lại của mô hình có thể bị thiếu chính xác. Tuy nhiên, trong hồi quy nhị phân không có chỗ cho sự thiếu chính xác bởi vì phương trình mô hình chỉ bao gồm trung bình (= xác suất) và khả năng là trung bình và 1 - trung bình, tương ứng. Điều này trái ngược với hồi quy dữ liệu tuyến tính hoặc đếm, có thể có tính không đồng nhất, quá mức, vv ..