Regression table for bayesglm?

Question

stargazer is a great tool to generate a regression table if you are not using bayesglm. For example, suppose I have the following data:

library(dplyr)
set.seed(9782)
N<-1000
df1 <- data.frame(v1=sample(c(0,1),N,replace = T),
                  v2=sample(c(0,1),N,replace = T),
                  Treatment=sample(c("A", "B", "C"), N, replace = T),
                  noise=rnorm(N)) %>% 
  mutate(Y=0.5*v1-0.7*v2+2*I(Treatment=="B")+3*I(Treatment=="C")+noise)

I can run lm and then create html (or text) output for my r markdown:

lm(data = df1, formula = Y~Treatment+v1+v2) %>% 
  stargazer::stargazer(type="html", style = "qje")

Is there a way to do something similar for bayesglm? In this case, pointEstimate has the coefficients and standardError the standard errors

library(arm)
fit <- bayesglm(data = df1, formula = Y~Treatment+v1+v2)
posteriorDraws <- coef(sim(fit, n.sims=5000))
pointEstimate <- colMeans(posteriorDraws)
standardError <- apply(posteriorDraws, 2, sd)

Answer 1

It looks like this does the trick:

library(texreg)
htmlreg(fit)

and for text:

screenreg(list(fit))

Answer 2

As @rawr points out in comments stargazer is -- while useful -- unfortunately written in a rather monolithic, hard-to-extend style. The broom package is (IMO) a nicely designed modular/object-oriented framework for extracting summary information from a large range of model types, but it's not oriented toward producing textual/tabular summaries. For those who like that sort of thing, it would be great if the stargazer front end could be grafted onto a broom back end, but it would be a lot of work. In the meantime, short of hacking the internal stargazer:::.stargazer.wrap function, this method (tricking stargazer into believing that fit is actually an lm() model) sort of works:

class(fit) <- c("glm","lm")
fit$call[1] <- quote(lm())
stargazer(fit)

It presents the coefficients and standard errors that are built into the fit object rather than the output of your posterior draws, but in this example at least those answers are extremely similar.

Answer 3

If you are OK with markdown, then the generic pander S3 method usually does a pretty good job:

> pander(fit, round = 4)

--------------------------------------------------------------
              Estimate   Std. Error   t value   Pr(>|t|) 
----------------- ---------- ------------ --------- ----------
 **(Intercept)**    0.0864      0.0763      1.131     0.2581  

 **TreatmentB**     1.951       0.0826      23.62       0     

 **TreatmentC**     3.007       0.0802      37.49       0     

     **v1**         0.4555      0.0659      6.915       0     

     **v2**        -0.6845      0.0659     -10.39       0     
--------------------------------------------------------------

Table: Fitting generalized (gaussian/identity) linear model: Y ~ Treatment + v1 + v2

Please note that I forced the numbers to be rounded in this example, as some P values were extremely low, so the default digits and other global options would result in an extremely wide table. But there are some othe params you might want to use, eg:

> pander(summary(fit), round = 4, add.significance.stars = TRUE, move.intercept = TRUE, summary = TRUE, split.cells = Inf)

----------------------------------------------------------------------
              Estimate   Std. Error   t value   Pr(>|t|)         
----------------- ---------- ------------ --------- ---------- -------
 **TreatmentB**     1.951       0.0826      23.62       0       * * * 

 **TreatmentC**     3.007       0.0802      37.49       0       * * * 

     **v1**         0.4555      0.0659      6.915       0       * * * 

     **v2**        -0.6845      0.0659     -10.39       0       * * * 

 **(Intercept)**    0.0864      0.0763      1.131     0.2581          
----------------------------------------------------------------------


(Dispersion parameter for  gaussian  family taken to be  1.083267 )


-------------------- -----------------------------------
   Null deviance:     2803  on 999  degrees of freedom  

 Residual deviance:   1078  on 995  degrees of freedom  
-------------------- -----------------------------------