Get coefficients estimated by maximum likelihood into a stargazer table

Question

Stargazer produces very nice latex tables for lm (and other) objects. Suppose I've fit a model by maximum likelihood. I'd like stargazer to produce a lm-like table for my estimates. How can I do this?

Although it's a bit hacky, one way might be to create a "fake" lm object containing my estimates -- I think this would work as long as summary(my.fake.lm.object) works. Is that easily doable?

An example:

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?

To be more precise: with lm objects, stargazer nicely prints the dependent variable at the top of the table, includes SEs in parentheses below the corresponding estimates, and has the R^2 and number of observations at the bottom of the table. Is there a(n easy) way to obtain the same behavior with a "custom" model estimated by maximum likelihood, as above?

Here are my feeble attempts at dressing up my optim output as a lm object:

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

I have attempted sth similar with the `texreg` package. Due to laziness, I ended up overwriting the coefficients and standard errors of a different model, which gave me the desired output. In your case, you could e.g. overwrite the coefficients and standard errors of `model1`. While this is not a sophisticated solution, it should work. Needless to say, I am curious to see if any better solutions come up... — coffeinjunky, Jun 29 '14 at 23:20
You can take a look at the stargazer function that does the heavy lifting with `stargazer:::.stargazer.wrap`. It looks like a container with a bunch of other functions in addition to the code that formats the tables. And it seems like it evaluates quite a few components for `lm` (and `glm`) that would make it very hard to dress up your `optim()` results. — andybega, Dec 19 '14 at 09:21
In ``texreg``, it would be sufficient to create a ``texreg`` object by using the ``createTexreg`` function. You basically just hand over the coefficients, SEs etc. See ``?createTexreg``. The ``texreg`` object can then be fed into the ``texreg``, ``htmlreg``, ``screenreg``, and ``plotreg`` functions. Alternatively, section 6 of the JSS article describes how to write and register methods for new model types in case you want to recycle the same template later on. — Philip Leifeld, Mar 14 '15 at 09:26

score 2 · Answer 1 · answered Dec 09 '15 at 00:37

I was just having this problem and overcame this through the use of the coef se, and omit functions within stargazer... e.g.

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

score 1 · Answer 2 · answered Dec 19 '16 at 12:32

You need to first instantiate a dummy lm object, then dress it up:

#...
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

(and then of course make sure the remaining summary stats are correct)

score 0 · Answer 3 · answered Jun 29 '15 at 18:54

I don't know how committed you are to using stargazer, but you can try using the broom and the xtable packages, the problem is that it won't give you the standard errors for the optim model

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

Get coefficients estimated by maximum likelihood into a stargazer table

3 Answers3

Linked