如何构建glm的参数和AICcWt的摘要table

Question

我有一个用于构建广义线性模型的数据集。响应变量是二元的 (absence/presence)，解释变量是分类的。

代码

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)

问题

我想构建一个像下面这样的 table 来总结我的候选集中模型的系数、标准误差和 Akaike 权重。

问题

谁能建议如何使用我的模型列表和 AIC table 最好地构建这个 table？

Answer 1

只是指出：broom 通过将模型输出转换为数据框，然后您可以对其进行整形，从而使您到达想要到达的位置的一半。

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

如果您必须坚持 table 的布局，我想出了以下（可以说是笨拙的）重新排列所有内容的方法：

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

或者，您可能希望自己熟悉用于显示模型输出以供发布的包，例如想到 texreg 或 stargazer：

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

如何构建glm的参数和AICcWt的摘要table

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

r

data-manipulation

glm