获取每个多面组的置信区间答案

【问题标题】：Get confidence intervals per faceted group获取每个多面组的置信区间
【发布时间】：2022-08-06 20:24:41
【问题描述】：

我正在使用lme4 包来运行线性混合效应模型。我想在 ggplot 中添加每个组级别的拟合线的置信区间。

我的数据： data 是一个数据框，包含： Plot_label：字符变量 // PD_avg：数字变量 // Year：因子 // GS_Prec：数字变量 // Direction：因子

我的代码如下：

#Run the model
mixed.lm <- lmer(PD_avg ~ log(GS_Prec) * Direction + (1|Plot_label) + (1|Year), data = data, REML=TRUE)

#Predict
pred1 <- predict(mixed.lm, newdata = data, re.form = NA) 

#Plot
ggplot(data, aes(log(GS_Prec), PD_avg, colour = Direction)) +
  geom_point(alpha = .2) +
  facet_wrap(~Direction) +
  geom_smooth(aes(y = pred1, colour = Direction), method = \"lm\", size = 1.5, se = T)

我在这里得到的图：

要添加 CI，我设置了 se = T，但它不起作用。所以我试图使用geom_ribbon，但它也不起作用。

我发现一个类似的主题有同样的问题（https://stats.stackexchange.com/questions/552734/r-plotting-lmer-confidence-intervals-per-faceted-group）。我确实关注了这个话题，顺便得到了一个意想不到的结果。

我的代码：

gr <- ref_grid(mixed.lm, cov.keep = c(\"GS_Prec\", \"Direction\"))
emm <- emmeans(gr, spec = c(\"GS_Prec\",\"Direction\"), level = 0.95)
emm

ggplot(data, aes(log(GS_Prec), PD_avg, colour = Direction)) +
  geom_point(alpha = .2) +
  facet_wrap(~Direction) +
  geom_smooth(aes(y = pred1, colour = Direction), method = \"lm\", size = 1.5) +
  geom_ribbon(data = data.frame(emm), aes(ymin = lower.CL, ymax = upper.CL, y = NULL, fill = Direction), alpha = 0.1)+
  geom_smooth(aes(y = pred1, colour = Direction), method = \"lm\", size = 1.5)

我想将置信区间的长度与点的范围联系起来。有谁知道如何正确表示 CI？

这是我的子集数据

data.1 <- data.frame(Plot_label = c(\"BT 1-1-3\", \"BT 1-1-3\", \"BT 1-2-1\", \"BT 1-2-1\",
                                    \"GW 1-1-1\", \"GW 1-1-1\", \"GW 1-5-2\", \"GW 1-5-2\",
                                    \"SP 1-5-2\", \"SP 1-5-2\", \"SP 2-8-2\", \"SP 2-8-2\"),
                     PD_avg = c(\"1196.61\", \"1323.15\", \"1172.17\", \"757.18\",
                                \"1516.02\", \"801.87\", \"1422.93\", \"1062.10\",
                                \"1580.51\", \"1520.30\", \"1326.25\", \"1321.89\"),
                     Year = c(\"2016\", \"2017\", \"2016\", 2017,
                              \"2016\", \"2017\", \"2016\", \"2017\",
                              \"2016\", \"2017\", \"2016\", \"2017\"),
                     Direction = c(\"BT-BT\", \"BT-BT\", \"BT-BT\", \"BT-BT\",
                                   \"GW-BT\", \"GW-BT\", \"GW-BT\", \"GW-BT\",
                                   \"SP-SP\", \"SP-SP\", \"SP-SP\", \"SP-SP\"),
                     GS_Prec = c(\"130.5\", \"190.5\", \"130.5\", \"190.5\",
                                 \"130.5\", \"190.5\", \"130.5\", \"190.5\",
                                 \"593.26\", \"480.29\", \"593.26\", \"593.26\"))

您能否使用dput 分享一些可重现的数据？
@Quinten嗨，我已经用我的子集数据更新了我的问题。你能检查一下吗？
请参阅getting confident interval in mixed effect models 上的@BenBolker 建议。
@AdamQuek 感谢您提供非常有用的链接。我试过了，但它对我不起作用。我得到了每个组的 CI，但是我无法在 ggplot 中手动绘制。他们一直说Error: Aesthetics must be either length 1 or the same as the data (162): ymin and ymax。
所以，我获取 CI CI <- as.data.frame(confint(mixed.lm, method=\"Wald\")) 的代码。然后我排除了NA 值。然后我确实使用了一堆代码来绘制：ggplot(data) + geom_point(aes(x = GS_Prec, y = PD_avg, colour = Direction)) + facet_wrap(~Direction) + geom_ribbon(data = CI, aes( ymin = CI$2.5 %, ymax = CI$97.5 %), alpha = 0.5)

标签： r ggplot2

【解决方案1】：

您还可以使用 ggeffects 包中的 ggpredict 函数，如下所示：

data <- data.frame(Plot_label = c("BT 1-1-3", "BT 1-1-3", "BT 1-2-1", "BT 1-2-1",
                                    "GW 1-1-1", "GW 1-1-1", "GW 1-5-2", "GW 1-5-2",
                                    "SP 1-5-2", "SP 1-5-2", "SP 2-8-2", "SP 2-8-2"),
                     PD_avg = c("1196.61", "1323.15", "1172.17", "757.18",
                                "1516.02", "801.87", "1422.93", "1062.10",
                                "1580.51", "1520.30", "1326.25", "1321.89"),
                     Year = c("2016", "2017", "2016", "2017",
                              "2016", "2017", "2016", "2017",
                              "2016", "2017", "2016", "2017"),
                     Direction = c("BT-BT", "BT-BT", "BT-BT", "BT-BT",
                                   "GW-BT", "GW-BT", "GW-BT", "GW-BT",
                                   "SP-SP", "SP-SP", "SP-SP", "SP-SP"),
                     GS_Prec = c("130.5", "190.5", "130.5", "190.5",
                                 "130.5", "190.5", "130.5", "190.5",
                                 "593.26", "480.29", "593.26", "593.26"))

library(lme4)
library(ggplot2)

# make columns numeric
data$GS_Prec <- as.numeric(data$GS_Prec)
data$PD_avg <- as.numeric(data$PD_avg)

#Run the model
mixed.lm <- lmer(PD_avg ~ log(GS_Prec) * Direction + (1|Plot_label) + (1|Year), data = data, REML=TRUE)
#> boundary (singular) fit: see help('isSingular')

library(ggeffects)

#Predict
pred1 <- ggpredict(mixed.lm, c("GS_Prec", "Direction")) 

#Plot
plot(pred1, add.data = TRUE)

^{由reprex package (v2.0.1) 于 2022 年 7 月 14 日创建}

【讨论】：

谢谢你。我使用了ggpredict 包，我也得到了相同的结果。如您所见，每条拟合线的 CI 都被夸大了。我的意思是 CI 应该限制为 2 个点。你是这么想的吗？

【解决方案2】：

并且使用参数 limit.range 和 facets 允许您将预测限制在实际数据的范围内并为绘图创建事实：

data <- data.frame(Plot_label = c("BT 1-1-3", "BT 1-1-3", "BT 1-2-1", "BT 1-2-1",
                                    "GW 1-1-1", "GW 1-1-1", "GW 1-5-2", "GW 1-5-2",
                                    "SP 1-5-2", "SP 1-5-2", "SP 2-8-2", "SP 2-8-2"),
                     PD_avg = c("1196.61", "1323.15", "1172.17", "757.18",
                                "1516.02", "801.87", "1422.93", "1062.10",
                                "1580.51", "1520.30", "1326.25", "1321.89"),
                     Year = c("2016", "2017", "2016", "2017",
                              "2016", "2017", "2016", "2017",
                              "2016", "2017", "2016", "2017"),
                     Direction = c("BT-BT", "BT-BT", "BT-BT", "BT-BT",
                                   "GW-BT", "GW-BT", "GW-BT", "GW-BT",
                                   "SP-SP", "SP-SP", "SP-SP", "SP-SP"),
                     GS_Prec = c("130.5", "190.5", "130.5", "190.5",
                                 "130.5", "190.5", "130.5", "190.5",
                                 "593.26", "480.29", "593.26", "593.26"))

library(lme4)
#> Loading required package: Matrix
library(ggplot2)

# make columns numeric
data$GS_Prec <- as.numeric(data$GS_Prec)
data$PD_avg <- as.numeric(data$PD_avg)

#Run the model
mixed.lm <- lmer(PD_avg ~ log(GS_Prec) * Direction + (1|Plot_label) + (1|Year), data = data, REML=TRUE)
#> boundary (singular) fit: see help('isSingular')

library(ggeffects)

#Predict
pred1 <- ggpredict(mixed.lm, c("GS_Prec", "Direction")) 

#Plot
plot(pred1, add.data = TRUE, limit.range = TRUE)

# with facets
plot(pred1, add.data = TRUE, limit.range = TRUE, facets = TRUE)

^{由reprex package (v2.0.1) 于 2022 年 8 月 6 日创建}

【讨论】：