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

Question

我目前正在处理 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()？

提前致谢！

Answer 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)