Errori standard robusti diversi della regressione della logica in Stata e R

Question

Sto tentando di replicare una regressione logit da Stata a R. In Stata, utilizzo l'opzione "robusta" per ottenere il robusto errore standard (errore standard coerente eterodosso). Sono in grado di replicare esattamente gli stessi coefficienti di Stata, ma non sono in grado di avere lo stesso robusto errore standard con il pacchetto "sandwich".

Ho provato alcuni esempi di regressione lineare OLS; sembra che gli stimatori sandwich di R e Stata mi diano lo stesso robusto errore standard per OLS. Qualcuno sa come Stata calcola lo stimatore sandwich per la regressione non lineare, nel mio caso la regressione logit?

Grazie!

Codici allegati: in 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

Il predefinito cosiddetti "robusti" errori standard in Stata corrispondono a quanto sandwich() dal pacchetto omonimo calcola. L'unica differenza è come viene eseguita la regolazione del campione finito. Nella funzione sandwich(...) non viene eseguita alcuna regolazione del campione finito per impostazione predefinita, vale a dire, il sandwich è diviso per 1/n dove n è il numero di osservazioni. In alternativa, è possibile utilizzare sandwich(..., adjust = TRUE) che divide per 1/(n - k) dove k è il numero di regressori. E Stata divide per 1/(n - 1).

Ovviamente, questi non differiscono affatto. E a parte alcuni casi speciali (ad esempio, la regressione lineare OLS) non esiste alcun argomento per 1/(n - k) o 1/(n - 1) per funzionare "correttamente" in campioni finiti (ad esempio imparzialità). Almeno non al meglio delle mie conoscenze.

Quindi, per ottenere gli stessi risultati come in Stata si può fare fare:

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

Questo produce

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 

E solo per la cronaca: In caso di risposta binaria, questi "robusto" gli errori standard non sono robusti contro nulla. A condizione che il modello sia specificato correttamente, sono coerenti e va bene utilizzarli ma non si guardano da eventuali errori di identificazione nel modello. Poiché l'ipotesi di base per gli errori standard del sandwich da eseguire è che l'equazione del modello (o più precisamente la funzione corrispondente del punteggio) sia specificata correttamente mentre il resto del modello può essere errato. Tuttavia, in una regressione binaria non vi è spazio per l'errata specificazione perché l'equazione del modello consiste solo della media (= probabilità) e la probabilità è la media e 1 - media, rispettivamente. Ciò è in contrasto con la regressione dei dati lineare o di conteggio dove possono esserci eteroschedasticità, sovradosaggio, ecc.