Репликация Stata Probit с надежными ошибками в R

Я пытаюсь воспроизвести вывод Stata в R. Я использую набор данных дела. У меня проблемы с воспроизведением пробит-функции с надежными стандартными ошибками.

Код Stata выглядит так:

probit affair male age yrsmarr kids relig educ ratemarr, r

Я начал с:

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

Затем я пробовал различные настройки с помощью пакета sandwich, например:

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

Или (со всеми типами от HC0 до HC5):

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

Или это, как было предложено здесь (нужно ли мне вводить что-то другое для object?):

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

Ни одна из этих попыток не привела к одинаковым стандартным ошибкам или z-значениям в выходных данных статистики.

Надеемся на конструктивные идеи!

Заранее спасибо!


r stata replicate robust
person Semprini    schedule 14.05.2015    source источник
comment
Взгляните на пример 5 здесь и параграф справа выше . Кроме того, если у вас есть гетероскедастические ошибки, этот подход последовательно оценивает стандартные ошибки предвзятых и противоречивых параметров. Многие думают, что это глупо.   -  person Dimitriy V. Masterov    schedule 15.05.2015
comment
Может быть, вы можете выложить полный код репликации вместе с выводом? В настоящее время мне не совсем ясно, какую версию данных вы использовали и каковы результаты в Stata и R соответственно.   -  person Achim Zeileis    schedule 15.05.2015
comment
Спасибо @Dimitriy V. Masterov за публикацию ваших результатов. Так что это не просто фактор, связанный с регулировкой степеней свободы. Код R / sandwich действительно идентичен (просто используются разные результаты make.link), поэтому я немного удивлен, что эта стратегия работает для репликации logit, но не для пробита. Не знаю, как такое могло случиться ...   -  person Achim Zeileis    schedule 15.05.2015
comment
@AchimZeileis Метод sandwich1 работает для стандартных ошибок в логит-регрессии с использованием тех же данных, но мне не удалось получить с ним правильный chi2. Есть идеи, почему это может быть?   -  person Semprini    schedule 16.05.2015
comment
Если вы имеете в виду статистику хи-квадрат из теста Вальда, то вы можете использовать функцию waldtest() из пакета lmtest или функцию linearHypothesis() из пакета car. Оба позволяют подключать необязательные vcov аргументы.   -  person Achim Zeileis    schedule 16.05.2015

Ответы (3)


arrow_upward
3
arrow_downward

Для людей, которые собираются запрыгнуть на эту повозку, вот код, демонстрирующий проблему (данные здесь):

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 дает:

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 дает:

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
------------------------------------------------------------------------------

Приложение:

Различие в ковариационной оценке коэффициентов связано с разными алгоритмами подбора. В R команда glm использует итеративный метод наименьших квадратов, а команда Stata probit использует метод машинного обучения, основанный на алгоритме Ньютона-Рафсона. Вы можете сопоставить то, что делает R, с glm в Stata с опцией irls:

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

Это дает:

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
------------------------------------------------------------------------------

Они будут близки, но не идентичны. Я не уверен, как заставить R использовать что-то вроде NR без большой работы.

person Dimitriy V. Masterov    schedule 14.05.2015
comment
Спасибо, что еще раз проиллюстрировали это! Поскольку у меня нет лицензии stata и есть только физическая распечатка, я не мог попробовать поэкспериментировать с данными на Stata. Кажется, что , r использует разные стандартные ошибки для пробита и логита, но у меня есть только базовые знания о Stata, поэтому я не могу это понять - person Semprini; 15.05.2015

arrow_upward
2
arrow_downward

Я использую матричный подход, как подробно описано здесь (стр.57) чтобы сопоставить результаты R со Stata. Однако я пока не смог точно сопоставить результат. Думаю, небольшая разница может быть из-за разницы в баллах. Очки в R совпадают с Stata только до 4 знаков после запятой.

Статистика

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

Это дает те же стандартные ошибки, что и:

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

Это соответствует:

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

Вывод: разница в результатах между R и Stata связана с разницей в баллах (соответствует только до 4 знаков после запятой).

person Metrics    schedule 16.05.2015
comment
Интересная догадка! К сожалению, я думаю, что это выходит за рамки моего понимания R, чтобы иметь возможность исправить это. - person Semprini; 16.05.2015

arrow_upward
2
arrow_downward

Чтобы замкнуть круг в этом обсуждении, можно сопоставить исходный вывод Stata в R, используя sampleSelection::probit для оценки и пакет sandwich (я использовал версию 2.5) для вычисления надежных стандартных ошибок. Функция probit использует максимальное правдоподобие, как и ее аналог в Stata.

Как и в исходном сообщении, код Stata

probit affair male age yrsmarr kids relig educ ratemarr, robust

который дает

------------------------------------------------------------------------------
             |               Robust
      affair |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        male |   .1888175   .1319271     1.43   0.152    -.0697548    .4473898
         age |  -.0243996   .0111237    -2.19   0.028    -.0462017   -.0025975
     yrsmarr |    .054608   .0189634     2.88   0.004     .0174405    .0917755
        kids |   .2080754   .1662426     1.25   0.211     -.117754    .5339049
       relig |  -.1860854   .0532405    -3.50   0.000    -.2904348    -.081736
        educ |   .0155052   .0263552     0.59   0.556    -.0361501    .0671605
    ratemarr |  -.2727101   .0533916    -5.11   0.000    -.3773558   -.1680644
       _cons |     .76416    .534335     1.43   0.153    -.2831173    1.811437
------------------------------------------------------------------------------

Код R, который дает те же результаты, -

library(AER)
library(sampleSelection)
data(Affairs)
Affairs$affair = Affairs$affairs > 0
Affairs$male = Affairs$gender == 'male'
reg = probit(affair ~ male + age + yearsmarried + children + religiousness +
           education + rating, data=Affairs)
print(coeftest(reg, vcovCL), digits=6)

Это дает

                Estimate Std. Error  t value   Pr(>|t|)    
(Intercept)    0.7641600  0.5343350  1.43011  0.1532109    
maleTRUE       0.1888175  0.1319271  1.43123  0.1528921    
age           -0.0243996  0.0111237 -2.19347  0.0286608 *  
yearsmarried   0.0546080  0.0189634  2.87966  0.0041248 ** 
childrenyes    0.2080755  0.1662426  1.25164  0.2111955    
religiousness -0.1860854  0.0532405 -3.49519  0.0005091 ***
education      0.0155052  0.0263552  0.58832  0.5565446    
rating        -0.2727101  0.0533916 -5.10773 4.4012e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Используя эти функции, обе вычисляют пробит-оценки максимального правдоподобия, и обе вычисляют устойчивые стандартные ошибки. В качестве отступления: снимаю шляпу перед авторами пакета sandwich, который (IMO) действительно очистил вычисления стандартных ошибок в R.

person Robert McDonald    schedule 24.08.2018