errore standard robusto e in cluster in R per regressione probit e logit

Question

Originariamente, io principalmente voglio eseguire un modello probit/logit con errore standard in cluster in R che è abbastanza intuitivo in Stata. Mi sono imbattuto con la risposta qui Logistic regression with robust clustered standard errors in R. Pertanto, ho provato a confrontare il risultato di Stata e di R con il robusto errore standard e l'errore standard del cluster. Ma ho notato che le uscite per entrambi gli errori standard tra software non sono esattamente le stesse. Tuttavia, se utilizzo il metodo suggerito qui https://diffuseprior.wordpress.com/2012/06/15/standard-robust-and-clustered-standard-errors-computed-in-r/. Posso ottenere l'output esatto sia da R che da Stata per la regressione lineare. Pertanto, temo che il codice che ho scritto in R non sia corretto e quale comando utilizzare se voglio eseguire un modello probit invece di un modello logit. O se esistono alternative eleganti per risolvere questo problema? Grazie.

codice R

codice

## 1. linear regression 
library(rms) 
# model<-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris) 
summary(model) 
fit=ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris) 
fit 
robcov(fit) #robust standard error 
robcov(fit, cluster=iris$Species) #clustered standard error 


## 2. logistic regression 
##demo data generation 
set.seed(1234) 
subj<-rep(1:20,each=4) 
con1<-rep(c(1,0),40) 
con2<-rep(c(1,1,0,0),20) 
effect<-rbinom(80,1,0.34) 
data<-data.frame(subj,con1,con2,effect) 
library(foreign);write.dta(data,'demo_data.dta') 

library(rms) 
fit=lrm(effect ~ con1 + con2, x=T, y=T, data=data) 
fit 
robcov(fit) ##robust standard error 
robcov(fit, cluster=data$subj) ## clustered standard error 

Stata

## 1. linear regression 
webuse iris 
reg seplen sepwid petlen petwid 
reg seplen sepwid petlen petwid,r 
reg seplen sepwid petlen petwid,cluster(iris) 


## 2. logistic regression 

use demo_data,clear 
logit effect con1 con2 
logit effect con1 con2,r 
logit effect con1 con2,cluster(subj) 

Answer 1

preferisco il pacchetto sandwich per calcolare gli errori standard robusti. Una ragione è la sua eccellente documentazione. Vedere vignette("sandwich") che mostra chiaramente tutti i valori predefiniti e le opzioni disponibili e il numero corresponding article che spiega come utilizzare ?sandwich con personalizzato bread e meat per casi speciali.

Possiamo usare sandwich per capire la differenza tra le opzioni che hai postato. La differenza sarà probabilmente il grado di correzione della libertà. Ecco un confronto per il semplice regressione lineare:

library(rms) 
library(sandwich) 

fitlm <-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris) 

#Your Blog Post: 
X <- model.matrix(fitlm) 
n <- dim(X)[1]; k <- dim(X)[2]; dfc <- n/(n-k)  
u <- matrix(resid(fitlm)) 
meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X 
Blog <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X)))) 

# rms fits: 
fitols <- ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris) 
Harrell <- sqrt(diag(robcov(fitols, method="huber")$var)) 
Harrell_2 <- sqrt(diag(robcov(fitols, method="efron")$var)) 

# Variations available in sandwich:  
variations <- c("const", "HC0", "HC1", "HC2","HC3", "HC4", "HC4m", "HC5") 
Zeileis <- t(sapply(variations, function(x) sqrt(diag(vcovHC(fitlm, type = x))))) 
rbind(Zeileis, Harrell, Harrell_2, Blog) 

      (Intercept) Sepal.Width Petal.Length Petal.Width 
const  0.2507771 0.06664739 0.05671929 0.1275479 
HC0   0.2228915 0.05965267 0.06134461 0.1421440 
HC1   0.2259241 0.06046431 0.06217926 0.1440781 
HC2   0.2275785 0.06087143 0.06277905 0.1454783 
HC3   0.2324199 0.06212735 0.06426019 0.1489170 
HC4   0.2323253 0.06196108 0.06430852 0.1488708 
HC4m  0.2339698 0.06253635 0.06482791 0.1502751 
HC5   0.2274557 0.06077326 0.06279005 0.1454329 
Harrell  0.2228915 0.05965267 0.06134461 0.1421440 
Harrell_2 0.2324199 0.06212735 0.06426019 0.1489170 
Blog  0.2259241 0.06046431 0.06217926 0.1440781 

1. Il risultato dalla voce del blog è equivalente a HC1. Se il post di blog è simile all'output Stata, Stata utilizza HC1.
2. La funzione di Frank Harrel produce risultati simili a HC0. Per quanto ho capito, questa è stata la prima soluzione proposta e quando si guarda attraverso vignette(sandwich) o gli articoli menzionati in ?sandwich::vcovHC, altri metodi hanno proprietà leggermente migliori. Differiscono nel loro grado di regolazione della libertà. Si noti inoltre che la chiamata a robcov(., method = "efron") è simile a HC3.

In ogni caso, se si desidera un output identico, utilizzare HC1 o semplicemente regolare la matrice varianza-covarianza in modo appropriato. Dopotutto, dopo aver osservato vignette(sandwich) per le differenze tra le diverse versioni, si vede che è necessario ridimensionare con una costante per passare da HC1 a HC0, che non dovrebbe essere troppo difficile. A proposito, si noti che HC3 o HC4 sono in genere preferiti a causa delle migliori proprietà dei campioni piccoli e del loro comportamento in presenza di osservazioni influenti. Pertanto, probabilmente si desidera modificare le impostazioni predefinite in Stata.

È possibile utilizzare queste matrici di varianza-covarianza fornendole a funzioni appropriate, come ad esempio lmtest::coeftest o car::linearHypothesis. Per esempio:

library(lmtest) 
coeftest(fitlm, vcov=vcovHC(fitlm, "HC1")) 

t test of coefficients: 

       Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.855997 0.225924 8.2151 1.038e-13 *** 
Sepal.Width 0.650837 0.060464 10.7640 < 2.2e-16 *** 
Petal.Length 0.709132 0.062179 11.4046 < 2.2e-16 *** 
Petal.Width -0.556483 0.144078 -3.8624 0.0001683 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

per errori standard cluster robusta, si dovrà regolare la carne del sandwich (vedi ?sandwich) o cercare un function farlo.Ci sono già severalsourcesexplaining in excruciatingdetailhowtodoitwithappropriatecodesorfunctions. Non c'è motivo per me di reinventare la ruota qui, quindi lo salterò.

C'è anche un relativamente nuovo e comodo errore di calcolo del pacchetto di errori standard di cluster per modelli lineari e modelli lineari generalizzati. Vedi here.