R: Wytrzymała diagnostyka SE i model w tabeli Stargazer

Question

Próbuję umieścić pewne wyjścia regresji 2SLS wygenerowane przez ivreg z pakietu AER w dokumencie Latex za pomocą pakietu stargazer. Mam jednak kilka problemów, których nie mogę rozwiązać sam.
1. Nie mogę wymyślić, jak wstawić diagnostykę modelu zgodnie z podsumowaniem ivreg. Mianowicie słabe testy instrumentów, Wu-Hausmann i Sargan Test. Chciałbym mieć je ze statystykami zwykle zgłaszanymi pod tabelą, takimi jak liczba obserwacji, R kwadrat i Resid. SE. Funkcja stargazer wydaje się nie mieć argumentu, w którym można podać listę z dodatkową diagnostyką. Nie dodałem tego do mojego przykładu, ponieważ szczerze mówiąc nie mam pojęcia, od czego zacząć.
2. Chcę wymieniać zwykłe standardowe błędy z solidnymi standardowymi błędami i jedynym sposobem, aby to zrobić, jest tworzenie obiektów z solidnymi standardowymi błędami i dodawanie ich w funkcji stargazer z se = list(). Umieszczam to w poniższym minimalnym przykładzie roboczym. Czy istnieje być może bardziej elegancki sposób na zakodowanie tego kodu, a może na ponowne jego odświeżenie i zapisanie go z solidnymi błędami standardowymi? Pomoc doceniona.

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

Oto jeden ze sposobów, aby robić to, co chcesz:

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

Który przyniesie:

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

Jak widać, ta umożliwia ręczne tym diagnostyki w poszczególnych modelach.

---

Można zautomatyzować to podejście poprzez stworzenie funkcji, które odbywają się na liście modeli (np list(summ.fit1, summ.fit2)) i wysyła przedmiotów wymaganych przez se lub add.lines argumentów.

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) 

Will zarówno podjąć w list z coeftest obiektów i otrzymując wektor SES jako oczekiwany przez se:

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

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

samo można zrobić w diagnostyce:

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) 

Both połączenia przyniosą:

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

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

Teraz wywołanie stargazer() staje się tak proste, jak ten, uzyskując identyczną moc jak wyżej:

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