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

Question

情况

我正在拟合一系列不断发展的回归模型。出于这个问题的目的，我们可以从模型 A、模型 B 和模型 C 的角度来考虑这些模型。所有模型共享至少一个相同的协变量。

我也在为这些模型拟合两年的数据。同样，为了这个问题的目的，年份将是 2000 年和 2010 年。

为了简化结果报告，我试图将回归报告合并为一个 table，格式如下：


                     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 包来实现上述 table。使用 stargazer，我可以获得多年来单个模型公式的完整格式 table（例如，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)

Answer 1

开始前的两个注意事项：

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

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

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

broom::tidy 一个采用统计模型的函数和 returns 一个 data.frame 估计值，每个系数一行。
broom::glance 一个采用统计模型的函数和 returns 具有模型特征（例如，观察数）的单行 data.frame
modelsummary_list 包含名为“tidy”和“glance”的 2 个元素的列表，class 名称为“modelsummary_list”。

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

编辑：在 modelsummary 中推荐的方法是使用 group 参数。滚动到此 post 的末尾以获取说明性代码。

带有有用讨论的过时示例

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, data = mtcars),
    lm(mpg ~ cyl + wt, data = mtcars2010)),
  "C" = list(
    lm(mpg ~ cyl + wt + disp, data = mtcars),
    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)

最后我们用stacking="vertical"参数调用modelsummary_wide得到这个table:

modelsummary_wide(models, stacking = "vertical")

当然，table 可以调整，系数重命名等，使用 modelsummary_wide 函数的其他参数或 kableExtra 或 output 参数。

没有详细解释的更现代的例子

library("modelsummary")
library("broom")
library("quantreg")


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

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

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)
    # output
    out <- list(tidy = ti, glance = data.frame("nobs 2010" = length(model_list[[1]]$fitted.values)))
    class(out) <- "modelsummary_list"
    return(out)
}

models <- lapply(models, tidy_model)

modelsummary(models, 
             group = model + term ~ group, 
             statistic = "conf.int")

		2000	2010
A	(Intercept)	36.800	36.800
		[30.034, 42.403]	[30.034, 42.403]
	cyl	-2.700	-67.944
		[-3.465, -1.792]	[-87.204, -45.102]
B	(Intercept)	38.871	38.871
		[30.972, 42.896]	[30.972, 42.896]
	cyl	-1.743	-43.858
		[-2.154, -0.535]	[-54.215, -13.472]
	wt	-2.679	-2.679
		[-5.313, -1.531]	[-5.313, -1.531]
C	(Intercept)	40.683	40.683
		[31.235, 47.507]	[31.235, 47.507]
	cyl	-1.993	-50.162
		[-3.137, -1.322]	[-78.948, -33.258]
	wt	-2.937	-2.937
		[-5.443, -1.362]	[-5.443, -1.362]
	disp	0.003	0.003
		[-0.009, 0.035]	[-0.009, 0.035]

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

Best way to report multiple regression models on several dimensions (evolving model formulation and year of data)

r

kable

stargazer

broom

modelsummary

带有有用讨论的过时示例

没有详细解释的更现代的例子