Логит биномиальная регрессия с кластерными стандартными ошибками

Question

Я пытаюсь повторить оценку glm из stata:

sysuse auto
logit foreign weight mpg, cluster(rep78)


Logistic regression                               Number of obs   =         69
                                                  Wald chi2(2)    =      31.57
                                                  Prob > chi2     =     0.0000
Log pseudolikelihood = -22.677963                 Pseudo R2       =     0.4652

                                  (Std. Err. adjusted for 5 clusters in rep78)
------------------------------------------------------------------------------
         |               Robust
 foreign |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
 -------------+----------------------------------------------------------------
  weight |  -.0042281   .0008022    -5.27   0.000    -.0058004   -.0026558
     mpg |  -.1524015   .0271285    -5.62   0.000    -.2055724   -.0992306
   _cons |   14.21689   2.031826     7.00   0.000     10.23458    18.19919

Я пробовал разные ответы, которые нашел в Интернете:

data <- read.dta("http://www.stata-press.com/data/r9/auto.dta")
data <- data[,c("foreign", "weight", "mpg", "rep78")]
data <- na.omit(data)
logit.model <- glm(foreign ~ weight + mpg ,  data, family = binomial(logit))
clx <- function(fm, dfcw, cluster){
     library(sandwich);library(lmtest)
     M <- length(unique(cluster))   
     N <- length(cluster)           
     K <- fm$rank                        
     dfc <- (M/(M-1))*((N-1)/(N-K))  
     uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
     vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
     coeftest(fm, vcovCL) 
}
clx(logit.model, 1, data$rep78)

z test of coefficients:

               Estimate  Std. Error z value  Pr(>|z|)    
(Intercept) 14.21688944  2.06235711  6.8935 5.443e-12 ***
weight      -0.00422810  0.00081428 -5.1924 2.076e-07 ***
mpg         -0.15240148  0.02753611 -5.5346 3.119e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

или же:

logit.model <- lrm(foreign ~ weight + mpg, x=T, y=T, data=data)
robcov(logit.model, cluster=data$rep78)

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept 14.2169 1.8173  7.82  <0.0001 
weight    -0.0042 0.0007 -5.89  <0.0001 
mpg       -0.1524 0.0243 -6.28  <0.0001 

или же:

logit.model <- glm(foreign ~ weight + mpg ,  data, family = binomial(logit))
G <- length(unique(data$rep78))
N <- length(data$rep78)
dfa <- (G/(G - 1)) * (N - 1)/logit.model$df.residual
c_vcov <- dfa * vcovHC(logit.model, type = "HC1", cluster = "group", adjust = T)
coeftest(logit.model, vcov = c_vcov)

z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept) 14.2168894  5.0830412  2.7969 0.0051591 ** 
weight      -0.0042281  0.0010915 -3.8736 0.0001072 ***
mpg         -0.1524015  0.1127248 -1.3520 0.1763823  

Однако ни с одним из предыдущих я не получаю точно таких же стандартных ошибок. Мне было интересно, если я делаю что-то не так или есть другой пакет, который я мог бы использовать, чтобы получить точно такие же результаты, как в Stata.

Answer 1

Я был в состоянии повторить результаты в stata, выполнив:

clrobustse <- function(fit.model, clusterid) {
    rank=fit.model$rank
    N.obs <- length(clusterid)            
    N.clust <- length(unique(clusterid))  
    dfc <- N.clust/(N.clust-1)                    
    vcv <- vcov(fit.model)
    estfn <- estfun(fit.model)
    uj <- apply(estfn, 2, function(x) tapply(x, clusterid, sum))
    N.VCV <- N.obs * vcv
    ujuj.nobs  <- crossprod(uj)/N.obs
    vcovCL <- dfc*(1/N.obs * (N.VCV %*% ujuj.nobs %*% N.VCV))
    coeftest(fit.model, vcov=vcovCL)
}
clrobustse(logit.model, data$rep78)