森林图面网格比较多个模型的回归模型系数

Posted 2023-02-16

【中文标题】森林图面网格比较多个模型的回归模型系数

【英文标题】：Forest plot facet grid comparing regression model coefficients from multiple models

【发布时间】：2021-12-24 02:47:59

【问题描述】：

我目前正在处理 30 个具有相同列名但数值数据不同的数据集。我需要将线性混合模型和广义线性模型应用于数据集的每个实例，并将生成的固定效应系数绘制在森林图上。

数据目前的结构如下（为简单起见，对每个列表元素使用相同的数据集）：

library(lme4)

data_list <- list()

# There's definitely a better way of doing this through lapply(), I just can't figure out how
for (i in 1:30)
     data_list[[i]] <- tibble::as_tibble(mtcars) # this would originally load different data at every instance


compute_model_lmm <- function(data)
     lmer("mpg ~ hp + disp + drat + (1|cyl)", data = data)


result_list_lmm <- lapply(data_list, compute_model_lmm)

我现在正在做的是

library(modelsummary)

modelplot(result_list_lmm)+
     facet_wrap(~model) #modelplot() takes arguments/functions from ggplot2

这需要很长时间，但它确实有效。

现在，我想在同一个地块上比较另一个模型，如

compute_model_glm <- function(data)
     glm("mpg ~ hp + disp + drat + cyl", data = data)


result_list_glm <- lapply(data_list, compute_model_glm)

modelplot(list(result_list_lmm[[1]], result_list_glm[[1]]))

但对于情节的每个实例。

如何将其指定为modelplot()？

提前致谢！

【问题讨论】：

您能否澄清一下“我想在同一个地块上比较另一个模型但具有网格结构”的意思？我不清楚最终的情节应该是什么样子。

对不起，如果不清楚，希望现在更好

还是不确定。你想要一个像我下面答案中的情节吗？

是的，我就是这个意思

好的，太好了！如果它满足您的需要，请接受该答案，或者解释为什么它仍然不充分。

【参考方案1】：

modelplot 函数为您提供了一些绘制系数和区间的基本方法（例如，检查 facet 参数）。

然而，函数的真正威力来自于使用draw=FALSE 参数。在这种情况下，modelplot 会在一个方便的数据框中为您提供估计值，其中包含 modelplot 函数的所有重命名、稳健的标准错误和其他实用程序。然后，您可以使用该数据框通过ggplot2 自己进行绘图，以实现无限自定义。

library(modelsummary)
library(ggplot2)

results_lm <- lapply(1:10, function(x) lm(hp ~ mpg, data = mtcars)) |>
    modelplot(draw = FALSE) |>
    transform("Function" = "lm()")
results_glm <- lapply(1:10, function(x) glm(hp ~ mpg, data = mtcars)) |>
    modelplot(draw = FALSE) |>
    transform("Function" = "glm()")
results <- rbind(results_lm, results_glm) 

head(results)
          term   model estimate std.error conf.low conf.high Function
1  (Intercept) Model 1 324.0823   27.4333  268.056  380.1086     lm()
3  (Intercept) Model 2 324.0823   27.4333  268.056  380.1086     lm()
5  (Intercept) Model 3 324.0823   27.4333  268.056  380.1086     lm()
7  (Intercept) Model 4 324.0823   27.4333  268.056  380.1086     lm()
9  (Intercept) Model 5 324.0823   27.4333  268.056  380.1086     lm()
11 (Intercept) Model 6 324.0823   27.4333  268.056  380.1086     lm()

ggplot(results, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high)) +
    geom_pointrange(aes(color = Function), position = position_dodge(width = .5)) +
    facet_wrap(~model)

【讨论】：

请注意，我在示例中使用的是lmer()，而不是lm()（我认为它除了公式之外没有任何变化）