Search code examples
rdata-manipulationglm

How to build a summary table of glm's parameters and AICcWt

I have a dataset that I am using to build generalised linear models. The response variable is binary (absence/presence) and the explanatory variables are categorical.

CODE

library(tidyverse)
library(AICcmodavg)

# Data
set.seed(123)
t <- tibble(ID = 1:100,
            A = as.factor(sample(c(0, 1), 100, T)),
            B = as.factor(sample(c("black", "white"), 100, T)),
            C = as.factor(sample(c("pos", "neg", "either"), 100, T)))

# Candidate set of models - Binomial family because response variable
# is binary (0 for absent & 1 for present)
# Global model is A ~ B_black + C_either
m1 <- glm(A ~ 1, binomial, t)
m2 <- glm(A ~ B, binomial, t)
m3 <- glm(A ~ C, binomial, t)
m4 <- glm(A ~ B + C, binomial, t)

# List with all models
ms <- list(null = m1, m_B = m2, m_C = m3, m_BC = m4)

# Summary table
aic_tbl <- aictab(ms)

PROBLEM

I want to build a table like the one below that summarises the coefficients, standard errors, and Akaike weights of the models within my candidate set.

enter image description here

QUESTION

Can anyone suggest how to best build this table using my list of models and AIC table?

Solution

  • Just to point it out: broom gets you half-way to where you want to get by turning the model output into a dataframe, which you can then reshape.

    library(broom)
bind_rows(lapply(ms, tidy), .id="key")
    key        term             estimate std.error            statistic p.value
1  null (Intercept) -0.12014431182649532     0.200 -0.59963969517107030   0.549
2   m_B (Intercept)  0.00000000000000123     0.283  0.00000000000000433   1.000
3   m_B      Bwhite -0.24116205496397874     0.401 -0.60071814968372905   0.548
4   m_C (Intercept) -0.47957308026188367     0.353 -1.35892869678271544   0.174
5   m_C        Cneg  0.80499548069651150     0.507  1.58784953814722285   0.112
6   m_C        Cpos  0.30772282333522433     0.490  0.62856402205887851   0.530
7  m_BC (Intercept) -0.36339654526926718     0.399 -0.90984856337213305   0.363
8  m_BC      Bwhite -0.25083209866475475     0.408 -0.61515191157571303   0.538
9  m_BC        Cneg  0.81144822536950656     0.508  1.59682131202527056   0.110
10 m_BC        Cpos  0.32706970242195277     0.492  0.66527127770403538   0.506

    And if you must insist of the layout of your table, I came up with the following (arguably clumsy) way of rearranging everything:

    out <- bind_rows(lapply(ms, tidy), .id="mod")

t1 <- out %>% select(mod, term, estimate) %>% spread(term, estimate) %>% base::t
t2 <- out %>% select(mod, term, std.error) %>% spread(term, std.error) %>% base::t
rownames(t2) <- paste0(rownames(t2), "_std_e")
tmp <- rbind(t1, t2[-1,])
new_t <- as.data.frame(tmp[-1,])
colnames(new_t) <- tmp[1,]
new_t

    Alternatively, you may want to familiarise yourself with packages that are meant to display model output for publication, e.g. texreg or stargazer come to mind:

    library(texreg)
screenreg(ms)

==================================================
                null     m_B      m_C      m_BC   
--------------------------------------------------
(Intercept)      -0.12     0.00    -0.48    -0.36 
                 (0.20)   (0.28)   (0.35)   (0.40)
Bwhite                    -0.24             -0.25 
                          (0.40)            (0.41)
Cneg                                0.80     0.81 
                                   (0.51)   (0.51)
Cpos                                0.31     0.33 
                                   (0.49)   (0.49)
--------------------------------------------------
AIC             140.27   141.91   141.66   143.28 
BIC             142.87   147.12   149.48   153.70 
Log Likelihood  -69.13   -68.95   -67.83   -67.64 
Deviance        138.27   137.91   135.66   135.28 
Num. obs.       100      100      100      100    
==================================================
*** p < 0.001, ** p < 0.01, * p < 0.05