Obtenga coeficientes estimados por máxima probabilidad en una tabla de observación de estrellas

83

Stargazer produce tablas de látex muy agradables para objetos lm (y otros). Supongamos que encajo en un modelo por máxima probabilidad. Me gustaría que Stargazer produjera una tabla parecida a una película para mis estimaciones. ¿Cómo puedo hacer esto?

Aunque es un poco hacky, una forma podría ser crear un objeto lm "falso" que contenga mis estimaciones; creo que esto funcionaría siempre que el resumen (my.fake.lm.object) funcione. ¿Es eso fácilmente factible?

Un ejemplo:

library(stargazer)

N <- 200
df <- data.frame(x=runif(N, 0, 50))
df$y <- 10 + 2 * df$x + 4 * rt(N, 4)  # True params
plot(df$x, df$y)

model1 <- lm(y ~ x, data=df)
stargazer(model1, title="A Model")  # I'd like to produce a similar table for the model below

ll <- function(params) {
    ## Log likelihood for y ~ x + student's t errors
    params <- as.list(params)
    return(sum(dt((df$y - params$const - params$beta*df$x) / params$scale, df=params$degrees.freedom, log=TRUE) -
               log(params$scale)))
}

model2 <- optim(par=c(const=5, beta=1, scale=3, degrees.freedom=5), lower=c(-Inf, -Inf, 0.1, 0.1),
                fn=ll, method="L-BFGS-B", control=list(fnscale=-1), hessian=TRUE)
model2.coefs <- data.frame(coefficient=names(model2$par), value=as.numeric(model2$par),
                           se=as.numeric(sqrt(diag(solve(-model2$hessian)))))

stargazer(model2.coefs, title="Another Model", summary=FALSE)  # Works, but how can I mimic what stargazer does with lm objects?

Para ser más precisos: con objetos lm, Stargazer imprime muy bien la variable dependiente en la parte superior de la tabla, incluye SE entre paréntesis debajo de las estimaciones correspondientes y tiene el R ^ 2 y el número de observaciones en la parte inferior de la tabla. ¿Existe una manera (n fácil) de obtener el mismo comportamiento con un modelo "personalizado" estimado por máxima verosimilitud, como el anterior?

Aquí están mis débiles intentos de disfrazar mi salida óptima como un objeto lm:

model2.lm <- list()  # Mimic an lm object
class(model2.lm) <- c(class(model2.lm), "lm")
model2.lm$rank <- model1$rank  # Problematic?
model2.lm$coefficients <- model2$par
names(model2.lm$coefficients)[1:2] <- names(model1$coefficients)
model2.lm$fitted.values <- model2$par["const"] + model2$par["beta"]*df$x
model2.lm$residuals <- df$y - model2.lm$fitted.values
model2.lm$model <- df
model2.lm$terms <- model1$terms  # Problematic?
summary(model2.lm)  # Not working
r optimization lm stargazer Adrian
fuente
6
He intentado algo similar con el texregpaquete. Debido a la pereza, terminé sobrescribiendo los coeficientes y errores estándar de un modelo diferente, lo que me dio el resultado deseado. En su caso, podría, por ejemplo, sobrescribir los coeficientes y errores estándar de model1. Si bien esta no es una solución sofisticada, debería funcionar. Ni que decir tiene, tengo curiosidad para ver si algunas soluciones mejores vienen ...
coffeinjunky
1
Puede echar un vistazo a la función Stargazer que hace el trabajo pesado stargazer:::.stargazer.wrap. Parece un contenedor con muchas otras funciones además del código que formatea las tablas. Y parece que evalúa bastantes componentes para lm(y glm) que harían muy difícil disfrazar sus optim()resultados.
andybega
3
En texreg, sería suficiente crear un texregobjeto utilizando la createTexregfunción. Básicamente, solo entrega los coeficientes, SE, etc. Vea ?createTexreg. El texregobjeto puede ser alimentado en el texreg, htmlreg, screenreg, y plotregfunciones. Alternativamente, la sección 6 del artículo de JSS describe cómo escribir y registrar métodos para nuevos tipos de modelos en caso de que desee reciclar la misma plantilla más adelante.
Philip Leifeld

Respuestas:

2

Estaba teniendo este problema y lo superé mediante el uso de las funciones coef se, y omitdentro de Stargazer ...

stargazer(regressions, ...
                     coef = list(... list of coefs...),
                     se = list(... list of standard errors...),
                     omit = c(sequence),
                     covariate.labels = c("new names"),
                     dep.var.labels.include = FALSE,
                     notes.append=FALSE), file="")
Curva
fuente
1

Primero debe crear una instancia de un lmobjeto ficticio y luego disfrazarlo:

#...
model2.lm = lm(y ~ ., data.frame(y=runif(5), beta=runif(5), scale=runif(5), degrees.freedom=runif(5)))
model2.lm$coefficients <- model2$par
model2.lm$fitted.values <- model2$par["const"] + model2$par["beta"]*df$x
model2.lm$residuals <- df$y - model2.lm$fitted.values
stargazer(model2.lm, se = list(model2.coefs$se), summary=FALSE, type='text')

# ===============================================
#                         Dependent variable:    
#                     ---------------------------
#                                  y             
# -----------------------------------------------
# const                        10.127***         
#                               (0.680)          
#                                                
# beta                         1.995***          
#                               (0.024)          
#                                                
# scale                        3.836***          
#                               (0.393)          
#                                                
# degrees.freedom              3.682***          
#                               (1.187)          
#                                                
# -----------------------------------------------
# Observations                    200            
# R2                             0.965           
# Adjusted R2                    0.858           
# Residual Std. Error       75.581 (df = 1)      
# F Statistic              9.076 (df = 3; 1)     
# ===============================================
# Note:               *p<0.1; **p<0.05; ***p<0.01

(y luego, por supuesto, asegúrese de que las estadísticas de resumen restantes sean correctas)

luffe
fuente
0

No sé qué tan comprometido está con el uso de stargazer, pero puede intentar usar los paquetes broom y xtable, el problema es que no le dará los errores estándar para el modelo optim 

library(broom)
library(xtable)
xtable(tidy(model1))
xtable(tidy(model2))
Andrelrms
fuente