Как построить сводную таблицу параметров glm и AICcWt - PullRequest
Новая
       52

Как построить сводную таблицу параметров glm и AICcWt

0 голосов
/ 31 декабря 2018

У меня есть набор данных, который я использую для построения обобщенных линейных моделей.Переменная ответа является двоичной (отсутствие / присутствие), а пояснительные переменные являются категориальными.

КОД 

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)

ПРОБЛЕМА

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

enter image description here

ВОПРОС

Кто-нибудь может подсказать, как лучше построить эту таблицу, используя мой список моделей и таблицу AIC?

1 Ответ

0 голосов
/ 31 декабря 2018

Просто чтобы указать на это: 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

А если выдолжен настаивать на макете вашей таблицы, я придумал следующий (возможно неуклюжий) способ перестановки всего: 

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
...