Replica di Stata Probit con errori gravi in ​​R

Question

Sto tentando di replicare un output Stata in R. Sto usando il set di dati affairs. Ho difficoltà a replicare la funzione probit con errori standard affidabili.

codice

Lo Stata sembra che:

probit affair male age yrsmarr kids relig educ ratemarr, r

ho iniziato con:

probit1 <- glm(affair ~ male + age + yrsmarr + kids + relig + educ + ratemarr, 
      family = binomial (link = "probit"), data = mydata) 

Poi ho provato varie regolazioni con il pacchetto sandwich, come ad esempio:

myProbit <- function(probit1, vcov = sandwich(..., adjust = TRUE)) { 
      print(coeftest(probit1, vcov = sandwich(probit1, adjust = TRUE))) 
} 

Oppure (con tutti i tipi HC0-HC5):

myProbit <- function(probit1, vcov = sandwich) { 
      print(coeftest(probit1, vcovHC(probit1, type = "HC0")) 
} 

O questo, come è stato suggerito here (devo inserire qualcosa di diverso per object):?

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

Nessuno di questi tentativi ha portato agli stessi errori standard o valori z dall'uscita stato.

Sperando in idee costruttive!

Grazie in anticipo!

Answer 1

Per le persone che stanno pensando di saltare su questo carro, qui è un codice che dimostra il problema (i dati here):

clear 
set more off 
capture ssc install bcuse 
capture ssc install rsource 
bcuse affairs 

saveold affairs, version(12) replace 

rsource, terminator(XXX) 
    library("foreign") 
    library("lmtest") 
    library("sandwich") 
    mydata<-read.dta("affairs.dta") 
    probit1<-glm(affair ~ male + age + yrsmarr + kids + relig + educ + ratemarr, family = binomial (link = "probit"), data = mydata) 
    sandwich1 <- function(object,...) sandwich(object) * nobs(object)/(nobs(object) - 1) 
    coeftest(probit1,vcov = sandwich1) 
XXX 

probit affair male age yrsmarr kids relig educ ratemarr, robust cformat(%9.6f) nolog 

R dà:

z test of coefficients: 

      Estimate Std. Error z value Pr(>|z|)  
(Intercept) 0.764157 0.546692 1.3978 0.1621780  
male   0.188816 0.133260 1.4169 0.1565119  
age   -0.024400 0.011423 -2.1361 0.0326725 * 
yrsmarr  0.054608 0.019025 2.8703 0.0041014 ** 
kids   0.208072 0.168222 1.2369 0.2161261  
relig  -0.186085 0.053968 -3.4480 0.0005647 *** 
educ   0.015506 0.026389 0.5876 0.5568012  
ratemarr -0.272711 0.053668 -5.0814 3.746e-07 *** 
--- 
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Stata rendimenti:

Probit regression        Number of obs  =  601 
               Wald chi2(7)  =  54.93 
               Prob > chi2  =  0.0000 
Log pseudolikelihood = -305.2525    Pseudo R2   =  0.0961 

------------------------------------------------------------------------------ 
      |    Robust 
     affair |  Coef. Std. Err.  z P>|z|  [95% Conf. Interval] 
-------------+---------------------------------------------------------------- 
     male | 0.188817 0.131927  1.43 0.152 -0.069755 0.447390 
     age | -0.024400 0.011124 -2.19 0.028 -0.046202 -0.002597 
    yrsmarr | 0.054608 0.018963  2.88 0.004  0.017441 0.091775 
     kids | 0.208075 0.166243  1.25 0.211 -0.117754 0.533905 
     relig | -0.186085 0.053240 -3.50 0.000 -0.290435 -0.081736 
     educ | 0.015505 0.026355  0.59 0.556 -0.036150 0.067161 
    ratemarr | -0.272710 0.053392 -5.11 0.000 -0.377356 -0.168064 
     _cons | 0.764160 0.534335  1.43 0.153 -0.283117 1.811437 
------------------------------------------------------------------------------ 

---

Addendum:

La differenza di covarianza stima dei coefficienti è dovuta alle diversi algoritmi di adattamento. In R, il comando glm utilizza il metodo iterativo dei minimi quadrati, mentre lo probit di Stata utilizza un metodo ML basato sull'algoritmo di Newton-Raphson. Si può corrispondere a ciò che sta facendo con R glm in Stata con irls opzione:

glm affair male age yrsmarr kids relig educ ratemarr, irls family(binomial) link(probit) robust 

Questo produce:

Generalized linear models       No. of obs  =  601 
Optimization  : MQL Fisher scoring    Residual df  =  593 
        (IRLS EIM)      Scale parameter =   1 
Deviance   = 610.5049916     (1/df) Deviance = 1.029519 
Pearson   = 619.0405832     (1/df) Pearson = 1.043913 

Variance function: V(u) = u*(1-u)     [Bernoulli] 
Link function : g(u) = invnorm(u)    [Probit] 

                BIC    = -3183.862 

------------------------------------------------------------------------------ 
      |    Semirobust 
     affair |  Coef. Std. Err.  z P>|z|  [95% Conf. Interval] 
-------------+---------------------------------------------------------------- 
     male | 0.188817 0.133260  1.42 0.157 -0.072367 0.450002 
     age | -0.024400 0.011422 -2.14 0.033 -0.046787 -0.002012 
    yrsmarr | 0.054608 0.019025  2.87 0.004  0.017319 0.091897 
     kids | 0.208075 0.168222  1.24 0.216 -0.121634 0.537785 
     relig | -0.186085 0.053968 -3.45 0.001 -0.291862 -0.080309 
     educ | 0.015505 0.026389  0.59 0.557 -0.036216 0.067226 
    ratemarr | -0.272710 0.053668 -5.08 0.000 -0.377898 -0.167522 
     _cons | 0.764160 0.546693  1.40 0.162 -0.307338 1.835657 
------------------------------------------------------------------------------ 

Questi sarà vicino, anche se non identici. Non sono sicuro di come ottenere R per utilizzare qualcosa come NR senza un sacco di lavoro.

Answer 2

Sto utilizzando l'approccio a matrice come descritto nel dettaglio here (p.57) per abbinare i risultati R con Stata. Tuttavia, non ho potuto ancora confrontare esattamente il risultato. Penso che la piccola differenza potrebbe essere dovuta alla differenza nei punteggi. I punteggi in R corrispondono a Stata solo fino a 4 cifre decimali.

Stata

clear all 
bcuse affairs 

probit affair male age yrsmarr kids relig educ ratemarr 
mat var_nr=e(V) 
predict double u, score 
matrix accum s = male age yrsmarr kids relig educ ratemarr [iweight=u^2*601/600] //n=601,n-1=600 
matrix rv = var_nr*s*var_nr 
mat diagrv=vecdiag(rv) 
matmap diagrv rse,m(sqrt(@)) //install matmap 
mat list rse //standard errors 

Questo ti dà gli stessi errori standard come:

qui probit affair male age yrsmarr kids relig educ ratemarr,r 



rse[1,8] 
     affair: affair: affair: affair: affair: affair: affair: affair: 
     male  age yrsmarr  kids  relig  educ ratemarr  _cons 
r1 .13192707 .01112372 .01896336 .16624258 .05324046 .02635524 .05339163 .53433495 

R:

library(AER) # Affairs data 
data(Affairs) 
mydata<-Affairs 
mydata$affairs<-with(mydata,ifelse(affairs>0,1,affairs)) # convert to 1 and 0 
probit1<-glm(affairs ~ gender+ age + yearsmarried + children + religiousness+education + rating,family = binomial(link = "probit"),data = mydata) 
u<-subset(estfun(probit1),select="(Intercept)") #scores: perfectly matches to 4 decimals with Stata: difference may be due to this step 
w0<-u%*%t(u)*(601/600) #(n/n-1) 
iweight<-matrix(0,nrow=601,ncol=601) #perfectly matches to 4 decimals with Stata 
diag(iweight)<-diag(w0) 
x<-model.matrix(probit1) 
s<-t(x)%*%iweight%*%x #doesn't match with Stata : 
rv<-vcov(probit1)%*%s%*%vcov(probit1) 
rse<-sqrt(diag(rv)) # standard errors 
    rse 
    (Intercept) gendermale   age yearsmarried childrenyes religiousness  education  rating 
    0.54669177 0.13325951 0.01142258 0.01902537 0.16822161 0.05396841 0.02638902 0.05366828 

Questo corrisponde con:

012.351.

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

Conclusione: la differenza nei risultati tra R e Stata è dovuta alla differenza nei punteggi (corrisponde solo fino a 4 posizioni decimali).