在多个维度上报告多重回归模型的最佳方法不断发展的模型公式和数据年份

问题描述

情况

我正在拟合一系列不断发展的回归模型。就这个问题而言，我们可以将这些模型视为模型 A、模型 B 和模型 C。所有模型至少共享一个相同的协变量。

我还将这些模型拟合为两个不同年份的数据。同样，对于这个问题，年份将是 2000 年和 2010 年。

为了简化结果报告，我尝试将回归报告合并到一个表格中，该表格具有以下某种格式：


                     2000        2010
Model A

    Coef Ex1
 
Model B

    Coef Ex1

    Coef Ex2

Model C

    Coef Ex1

    Coef Ex2

    Coef Ex3

这个想法是，人们可以跨多个模型和年份快速查看 Coef Ex1。

我尝试了什么

我尝试使用 R stargazer 和 kable 包来实现上表。使用 stargazer，我可以获得多年内单个模型公式的完全格式化表（例如，stargazer(modelA2000,modelA2010)，但我无法弄清楚如何在行上堆叠其他模型公式。

对于 kable，我已经能够堆叠水平模型，但我无法添加额外的年份（例如 coefs <- bind_rows(tidy(modelA2000),tidy(modelB2000),tidy(modelC2000)); coefs %>% kable()）。

问题：我如何使用 stargazer 或 kable 来报告行中不断变化的回归模型（共享相同的协变量）以及列中的横截面年份？我想我可以以某种方式扩展 here 发布的答案，尽管我不确定如何。

可重现的示例


# Load the data
mtcars <- mtcars

# Create example results for models A,B,and C for 2000
modelA2000 <- lm(mpg ~ cyl,data = mtcars)
modelB2000 <- lm(mpg ~ cyl + wt,data = mtcars)
modelC2000 <- lm(mpg ~ cyl + wt + disp,data = mtcars)

# Slightly modify data for second set of results
mtcars$cyl <- mtcars$cyl*runif(1)

# Fit second set of results. Same models,pretending it's a different year. 
modelA2010 <- lm(mpg ~ cyl,data = mtcars)
modelB2010 <- lm(mpg ~ cyl + wt,data = mtcars)
modelC2010 <- lm(mpg ~ cyl + wt + disp,data = mtcars)

解决方法

开始前的两个注意事项：

您想要一个漂亮的“自定义”表，因此几乎不可避免地需要一些手动操作。
我的回答依赖于 modelsummary 的开发版本，您可以像这样安装：

library(remotes)
install_github("vincentarelbundock/modelsummary")

我们需要 4 个概念，其中许多与 broom 包相关：

broom::tidy 一个函数，它采用统计模型并返回一个估计数据框，每个系数一行。
broom::glance 一个采用统计模型并返回具有模型特征（例如，观察次数）的单行数据框的函数
modelsummary_list 一个包含 2 个元素的列表，名为“tidy”和“glance”，类名为“modelsummary_list”。

modelsummary 包允许您绘制回归表。在幕后，它使用 broom::tidy 和 broom::glance 从这些模型中提取信息。用户还可以通过提供我们分配类 modelsummary_list 的列表来提供他们自己的模型信息，如 documented here.

modelsummary_wide 是一个函数，最初设计用于“堆叠”来自具有多组系数的多个模型的结果。这对于多项模型之类的东西很有用，但它也有助于我们在您的情况下，您在多个组中拥有多个模型（此处：年）。

首先，我们加载包、调整数据并估计我们的模型：

library(modelsummary)
library(broom)
library(dplyr)

mtcars2010 <- mtcars
mtcars2010$cyl <- mtcars$cyl * runif(1)

models <- list(
  "A" = list(
    lm(mpg ~ cyl,data = mtcars),lm(mpg ~ cyl,data = mtcars2010)),"B" = list(
    lm(mpg ~ cyl + wt,lm(mpg ~ cyl + wt,"C" = list(
    lm(mpg ~ cyl + wt + disp,lm(mpg ~ cyl + wt + disp,data = mtcars2010)))

请注意，我们将模型保存在三组列表中。

然后，我们定义一个 tidy_model 函数，该函数接受两个模型的列表（每年一个），结合这两个模型的信息，并创建一个 modelsummary_list对象（再次请参阅documentation）。请注意，我们将“年份”信息分配给 tidy 对象中的“组”列。

我们使用 lapply 将此函数应用于三组模型中的每组。

tidy_model <- function(model_list) {
    # tidy estimates
    tidy_2000 <- broom::tidy(model_list[[1]])
    tidy_2010 <- broom::tidy(model_list[[2]])
    # create a "group" column
    tidy_2000$group <- 2000
    tidy_2010$group <- 2010
    ti <- bind_rows(tidy_2000,tidy_2010)
    # glance estimates
    gl <- data.frame("N" = stats::nobs(model_list[[1]]))
    # output
    out <- list(tidy = ti,glance = gl)
    class(out) <- "modelsummary_list"
    return(out)
}

models <- lapply(models,tidy_model)

最后，我们使用 modelsummary_wide 参数调用 stacking="vertical" 来获取此表：

modelsummary_wide(models,stacking = "vertical")

当然，可以使用 modelsummary_wide 函数的其他参数或 kableExtra 或 output 参数支持的其他包来调整表、重命名系数等。

broom kable modelsummary r r stargazer