Verschiedene robuste Standardfehler der Logit-Regression in Stata und R

Question

Ich versuche, eine Logit-Regression von Stata nach R zu replizieren. In Stata verwende ich die Option "robust", um den robusten Standardfehler (Heteroskedastizität-konsistenter Standardfehler) zu haben. Ich bin in der Lage, die genau gleichen Koeffizienten von Stata zu replizieren, aber ich kann nicht den gleichen robusten Standardfehler mit dem Paket "Sandwich" haben.

Ich habe einige OLS lineare Regressionsbeispiele versucht; es scheint, als ob die Sandwich-Schätzer von R und Stata mir den gleichen robusten Standardfehler für OLS geben. Weiß jemand, wie Stata den Sandwichschätzer für die nichtlineare Regression berechnet, in meinem Fall die Logit-Regression?

Vielen Dank!

Befestigt Codes: 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

Der Standard sogenannte "robust" Standardfehler in Stata entspricht das, was sandwich() aus dem Paket mit dem gleichen Namen berechnet. Der einzige Unterschied besteht darin, wie die Finite-Sample-Anpassung durchgeführt wird. In der sandwich(...)-Funktion wird standardmäßig keine finite Abtastwerteinstellung vorgenommen, d. H. Das Sandwich wird durch 1/n geteilt, wobei n die Anzahl der Beobachtungen ist. Alternativ kann sandwich(..., adjust = TRUE) verwendet werden, das durch 1/(n - k) teilt, wobei k die Anzahl der Regressoren ist. Und Stata teilt sich durch 1/(n - 1).

Natürlich unterscheiden sich diese asymptotisch überhaupt nicht. Und mit Ausnahme einiger weniger spezieller Fälle (z. B. lineare Regression durch OLS) gibt es kein Argument dafür, dass 1/(n - k) oder 1/(n - 1) in finiten Proben "korrekt" arbeiten (z. B. Unverzerrtheit). Zumindest nicht nach meinem besten Wissen.

So die gleichen Ergebnisse zu erhalten, wie es in Stata können Sie tun, tun:

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

Dies ergibt

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 

Und nur für das Protokoll: In dem binären Antwort Fall dieser "robust" Standardfehler sind gegen nichts robust. Vorausgesetzt, dass das Modell korrekt angegeben ist, sind sie konsistent und es ist in Ordnung, sie zu verwenden, aber sie schützen nicht vor einer falschen Angabe im Modell. Da die Grundannahme für die Sandwich-Standardfehler ist, dass die Modellgleichung (genauer gesagt die entsprechende Bewertungsfunktion) korrekt spezifiziert ist, während der Rest des Modells falsch spezifiziert sein kann. In einer binären Regression gibt es jedoch keinen Spielraum für eine Fehlspezifizierung, da die Modellgleichung nur aus dem Mittelwert (= Wahrscheinlichkeit) und der Wahrscheinlichkeit den Mittelwert bzw. 1 - Mittelwert besteht. Dies steht im Gegensatz zur linearen oder Zähldatenregression, bei der Heteroskedastizität, Überdispersion usw. auftreten können.