0

情况

我正在拟合一系列不断发展的回归模型。出于这个问题的目的，我们可以根据模型 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多个模型和年份。

我试过什么

我尝试使用 Rstargazer和kable包来实现上表。有了stargazer我可以获得跨多年的单个模型公式的完全格式化的表格（例如，stargazer(modelA2000, modelA2010)但我不知道如何在行上堆叠其他模型公式。

因为kable我已经能够堆叠水平模型，但我无法再添加几年（例如， coefs <- bind_rows(tidy(modelA2000), tidy(modelB2000), tidy(modelC2000)); coefs %>% kable()）。

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

可重现的例子


# 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)
4

1 回答 1

0

开始前的两个注意事项：

  1. 您想要一个漂亮的“自定义”表，因此几乎不可避免地需要进行一些手动操作。
  2. 我的答案依赖于 的开发版本modelsummary，您可以像这样安装它：
library(remotes)
install_github("vincentarelbundock/modelsummary")

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

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

该modelsummary软件包允许您绘制回归表。在后台，它使用broom::tidy并broom::glance从这些模型中提取信息。用户还可以通过提供我们分配类的列表来提供他们自己的模型信息modelsummary_list，如此处所述。

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

带有有用讨论的过时示例

该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对象（再次，请参阅文档）。请注意，我们将“年份”信息分配给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。

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

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
一个 （截距） 36.800 36.800
[30.034, 42.403] [30.034, 42.403]
圆柱体 -2.700 -67.944
[-3.465, -1.792] [-87.204，-45.102]
乙 （截距） 38.871 38.871
[30.972, 42.896] [30.972, 42.896]
圆柱体 -1.743 -43.858
[-2.154, -0.535] [-54.215，-13.472]
重量 -2.679 -2.679
[-5.313，-1.53​​1] [-5.313，-1.53​​1]
C （截距） 40.683 40.683
[31.235, 47.507] [31.235, 47.507]
圆柱体 -1.993 -50.162
[-3.137，-1.322] [-78.948, -33.258]
重量 -2.937 -2.937
[-5.443, -1.362] [-5.443, -1.362]
显示 0.003 0.003
[-0.009, 0.035] [-0.009, 0.035]
于 2021-03-23T13:15:24.177 回答