Различные робастные стандартные ошибки регрессии логита в Stata и R

Question

Я пытаюсь повторить регрессию logit из Stata в R. В Stata я использую опцию "робастный", чтобы получить устойчивую стандартную ошибку (стандартную ошибку, согласованную с гетероскедастичностью). Я могу воспроизвести точно такие же коэффициенты из Stata, но у меня не может быть такой же надежной стандартной ошибки с пакетом "сэндвич".

Я пробовал некоторые примеры линейной регрессии OLS; похоже, что сэндвич-оценки R и Stata дают мне одинаковую здравую стандартную ошибку для OLS. Кто-нибудь знает, как Stata вычисляет сэндвич-оценку для нелинейной регрессии, в моем случае - логит-регрессии?

Спасибо!

Коды прилагаются: в 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

Так называемые "робастные" стандартные ошибки в Stata по умолчанию соответствуют тому, что sandwich() из одноименной упаковки вычисляет. Единственное отличие состоит в том, как выполняется корректировка конечной выборки. в sandwich(...) Функция по умолчанию вообще не выполняет настройку конечной выборки, т. е. сэндвич делится на 1/n, где n - количество наблюдений. С другой стороны, sandwich(..., adjust = TRUE) может использоваться, который делит на 1/(n - k), где k - количество регрессоров. А Stata делится на 1/(n - 1).

Конечно, асимптотически они ничем не отличаются. И за исключением нескольких особых случаев (например, линейная регрессия OLS) нет аргумента для 1/(n - k) или 1 / (n - 1), чтобы "работать" правильно в конечных выборках (например, беспристрастность). По крайней мере, насколько мне известно.

Таким образом, чтобы получить те же результаты, что и в Stata, вы можете сделать:

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

Это дает

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

И просто для записи: в случае двоичного ответа эти "надежные" стандартные ошибки не являются устойчивыми ни к чему. При условии, что модель указана правильно, они согласованы, и их можно использовать, но они не защищают от каких-либо ошибок в модели. Поскольку основное допущение для стандартных ошибок сэндвича состоит в том, что уравнение модели (или, точнее, соответствующая функция оценки) правильно задано, в то время как остальная часть модели может быть задана неправильно. Однако в бинарной регрессии нет места для ошибочной спецификации, потому что модельное уравнение состоит только из среднего значения (= вероятность) и вероятности - это среднее значение и 1 - среднее значение, соответственно. Это контрастирует с линейной или счетной регрессией данных, где может быть гетероскедастичность, избыточная дисперсия и т. Д.