Таблица регрессии для Bayesglm?

stargazer — отличный инструмент для создания таблицы регрессии, если не используя bayesglm. Например, предположим, что у меня есть следующие данные:

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)

Я могу запустить lm, а затем создать вывод html (или текст) для моей уценки r:

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

Есть ли способ сделать что-то подобное для bayesglm? В этом случае pointEstimate имеет коэффициенты, а standardError — стандартные ошибки.

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)

r stargazer pander bayesglm
person Ignacio    schedule 25.03.2016    source источник
comment
Bayesglm нет в списке ?stargazer:::`stargazer models`   -  person rawr    schedule 25.03.2016
comment
господи, звездочет — это одна функция< /i>, возможно поэтому никто не сделал расширения для байесглм   -  person rawr    schedule 25.03.2016

Ответы (3)


arrow_upward
2
arrow_downward

Похоже, это помогает:

library(texreg)
htmlreg(fit)

и для текста:

screenreg(list(fit))
person Ignacio    schedule 25.03.2016

arrow_upward
1
arrow_downward

Как отмечает @rawr в комментариях, stargazer хоть и полезен, но, к сожалению, написан в довольно монолитном, трудно расширяемом стиле. Пакет broom представляет собой (IMO) прекрасно спроектированную модульную/объектно-ориентированную структуру для извлечения сводки. информацию из большого диапазона типов моделей, но она не ориентирована на создание текстовых/табличных сводок. Для тех, кто любит такие вещи, было бы здорово, если бы stargazer переднюю часть можно было привить к broom задней части, но это потребовало бы много работы. В то же время, за исключением взлома внутренней функции stargazer:::.stargazer.wrap, этот метод (обмануть stargazer, заставив его поверить, что fit на самом деле является моделью lm()) вроде как работает:

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

В нем представлены коэффициенты и стандартные ошибки, встроенные в объект fit, а не выходные данные ваших апостериорных рисунков, но в этом примере, по крайней мере, эти ответы очень похожи.

person Ben Bolker    schedule 25.03.2016
comment
вот что я думал сделать, но это не сработало - у меня не было второй строки. Однако мне нужно было quote(lm()), чтобы это сработало. - person rawr; 25.03.2016
comment
ага. Я закончил тем, что взломал .stargazer.wrap, пока не обнаружил, что звездочет (тьфу тьфу тьфу) использует специальную функцию model.identify(), которая смотрит на объект вызова, а также на класс ... - person Ben Bolker; 25.03.2016
comment
Да, именно поэтому texreg использует общую функцию extract, чтобы пользователи могли писать и регистрировать свои собственные методы. extract построен по тому же принципу, что и broom, хотя broom был разработан позже. Результирующие объекты texreg являются объектами S4, в которых хранятся соответствующие данные, такие как коэффициенты и т. д., а затем функции texreg, htmlreg, screenreg и т. д. могут использовать эти данные для построения таблиц. Вы можете либо разрешить функции texreg вызывать extract внутри, либо сохранить промежуточные результаты и манипулировать ими, прежде чем передать их texreg. - person Philip Leifeld; 03.08.2016

arrow_upward
1
arrow_downward

Если вас устраивает уценка, то общий метод pander S3 обычно делает довольно хорошая работа:

> 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

Обратите внимание, что в этом примере я принудительно округлил числа, так как некоторые значения P были чрезвычайно низкими, поэтому по умолчанию digits и другие глобальные параметры привели бы к очень широкой таблице. Но есть некоторые другие параметры, которые вы можете использовать, например:

> 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  
-------------------- -----------------------------------
person daroczig    schedule 25.03.2016