R: Robusto sistema di diagnostica per SE e modello nella tabella degli stargazer

Question

Provo a inserire alcuni output di regressione 2SLS generati tramite ivreg dal pacchetto AER in un documento Latex utilizzando il pacchetto stargazer. Ho comunque un paio di problemi che non riesco a risolvere da solo.
1. Non riesco a capire come inserire la diagnostica del modello come fornita dal riepilogo di ivreg. Vale a dire test degli strumenti deboli, Wu-Hausmann e Sargan Test. Mi piacerebbe averli con le statistiche solitamente riportate sotto la tabella come numero di osservazioni, R al quadrato e Resid. SE. La funzione Stargazer non sembra avere un argomento in cui è possibile fornire una lista con diagnostica aggiuntiva. Non ho inserito questo nel mio esempio perché onestamente non ho idea di dove cominciare.
2. Voglio scambiare gli errori standard normali con errori standard robusti e l'unico modo per farlo è quello di produrre oggetti con errori standard robusti e aggiungerli nella funzione di stargazer con se = list(). Inserisco questo nell'esempio di lavoro minimo di seguito. C'è forse un modo più elegante di codificarlo o magari di ricostituire il modello e salvarlo con solidi errori standard? Aiuto apprezzato

library(AER) 
library(stargazer) 

y <- rnorm(100, 5, 10) 
x <- rnorm(100, 3, 15) 
z <- rnorm(100, 3, 7) 
a <- rnorm(100, 1, 7) 
b <- rnorm(100, 3, 5) 

# Fitting IV models 
fit1 <- ivreg(y ~ x + a | 
      a + z, 
      model = TRUE) 
fit2 <- ivreg(y ~ x + a | 
      a + b + z, 
      model = TRUE) 

# Here are the se's and the diagnostics i want 
summary(fit1, vcov = sandwich, diagnostics=T) 
summary(fit2, vcov = sandwich, diagnostics=T) 

# Getting robust se's, i think HC0 is the standard 
# used with "vcov=sandwich" from the above summary 
cov1  <- vcovHC(fit1, type = "HC0") 
robust1  <- sqrt(diag(cov1)) 
cov2  <- vcovHC(fit2, type = "HC0") 
robust2  <- sqrt(diag(cov1)) 

# Create latex table 
stargazer(fit1, fit2, type = "latex", se=list(robust1, robust2)) 

Answer 1

Ecco un modo per fare ciò che si vuole:

require(lmtest) 

rob.fit1  <- coeftest(fit1, function(x) vcovHC(x, type="HC0")) 
rob.fit2  <- coeftest(fit2, function(x) vcovHC(x, type="HC0")) 
summ.fit1 <- summary(fit1, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T) 
summ.fit2 <- summary(fit2, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T) 

stargazer(fit1, fit2, type = "text", 
      se = list(rob.fit1[,"Std. Error"], rob.fit2[,"Std. Error"]), 
      add.lines = list(c(rownames(summ.fit1$diagnostics)[1], 
          round(summ.fit1$diagnostics[1, "p-value"], 2), 
          round(summ.fit2$diagnostics[1, "p-value"], 2)), 
          c(rownames(summ.fit1$diagnostics)[2], 
          round(summ.fit1$diagnostics[2, "p-value"], 2), 
          round(summ.fit2$diagnostics[2, "p-value"], 2)))) 

che produrrà:

========================================================== 
            Dependent variable:  
           ---------------------------- 
              y    
            (1)   (2)  
---------------------------------------------------------- 
x         -1.222  -0.912  
           (1.672)  (1.002) 

a         -0.240  -0.208  
           (0.301)  (0.243) 

Constant       9.662   8.450** 
           (6.912)  (4.222) 

---------------------------------------------------------- 
Weak instruments     0.45   0.56  
Wu-Hausman       0.11   0.18  
Observations      100   100  
R2        -4.414  -2.458  
Adjusted R2      -4.526  -2.529  
Residual Std. Error (df = 97)  22.075  17.641  
========================================================== 
Note:       *p<0.1; **p<0.05; ***p<0.01 

Come si può vedere, questo permette compresi manualmente la diagnostica nei rispettivi modelli.

---

Si può automatizzare questo approccio creando una funzione che prende in un elenco di modelli (per esempio list(summ.fit1, summ.fit2)) ed emette gli oggetti necessari per se o add.lines argomenti.

gaze.coeft <- function(x, col="Std. Error"){ 
    stopifnot(is.list(x)) 
    out <- lapply(x, function(y){ 
     y[ , col] 
    }) 
    return(out) 
} 
gaze.coeft(list(rob.fit1, rob.fit2)) 
gaze.coeft(list(rob.fit1, rob.fit2), col=2) 

saranno entrambi prendere in list di coeftest oggetti, e cedere il vettore SE come previsto dalla se:

[[1]] 
(Intercept)   x   a 
    6.9124587 1.6716076 0.3011226 

[[2]] 
(Intercept)   x   a 
    4.2221491 1.0016012 0.2434801 

stesso può essere fatto per la diagnostica:

gaze.lines.ivreg.diagn <- function(x, col="p-value", row=1:3, digits=2){ 
    stopifnot(is.list(x)) 
    out <- lapply(x, function(y){ 
     stopifnot(class(y)=="summary.ivreg") 
     y$diagnostics[row, col, drop=FALSE] 
    }) 
    out <- as.list(data.frame(t(as.data.frame(out)), check.names = FALSE)) 
    for(i in 1:length(out)){ 
     out[[i]] <- c(names(out)[i], round(out[[i]], digits=digits)) 
    } 
    return(out) 
} 
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2) 
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), col=4, row=1:2, digits=2) 

Entrambi le chiamate produrranno:

$`Weak instruments` 
[1] "Weak instruments" "0.45"    "0.56"    

$`Wu-Hausman` 
[1] "Wu-Hausman" "0.11"  "0.18"  

Ora la chiamata stargazer() diventa semplice come questo, con produzione identico come sopra:

stargazer(fit1, fit2, type = "text", 
     se = gaze.coeft(list(rob.fit1, rob.fit2)), 
     add.lines = gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2))