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)

bayesglm isnt in the list `?stargazer:::\`stargazer models\`` — rawr, Mar 25 '16 at 17:15
good lord, stargazer is [_a single function_](https://github.com/cran/stargazer/blob/master/R/stargazer-internal.R#L8), perhaps that's why no one has made an extension for bayesglm — rawr, Mar 25 '16 at 17:16

score 2 · Accepted Answer · answered Mar 25 '16 at 17:09

2

It looks like this does the trick:

library(texreg)
htmlreg(fit)

and for text:

screenreg(list(fit))

answered Mar 25 '16 at 17:09

Ignacio

7,646
16
60
113

Ben Bolker · Answer 2 · 2016-03-25T17:45:27.617

1

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.

edited Mar 25 '16 at 17:45

answered Mar 25 '16 at 17:29

Ben Bolker

211,554
25
370
453

1

that's what I thought of doing but it didn't work--I didn't have the second line. I needed `quote(lm())` for this to work, though – rawr Mar 25 '16 at 17:43
yup. I ended up hacking through `.stargazer.wrap` until I found that stargazer is (ugh ugh ugh) using a special `model.identify()` function that looks at the call object as well as the class ... – Ben Bolker Mar 25 '16 at 17:45
Yeah, that's why ``texreg`` uses a generic function called ``extract``, so users can write and register their own methods. ``extract`` is built on a similar principle as ``broom``, though ``broom`` was developed later. The resulting ``texreg`` objects are S4 objects which store the relevant data like coefficients etc., and then the ``texreg``, ``htmlreg``, ``screenreg`` etc. functions can use those data to construct tables. You can either let the ``texreg`` function call ``extract`` internally, or you can save the intermediate results and manipulate them before you pass them on to ``texreg``. – Philip Leifeld Aug 03 '16 at 12:51

score 1 · Answer 3 · answered Mar 25 '16 at 19:47

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

> pander(fit, round = 4)

--------------------------------------------------------------
     &nbsp;        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)

----------------------------------------------------------------------
     &nbsp;        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  
-------------------- -----------------------------------

Regression table for bayesglm?

3 Answers3