Get confidence intervals per faceted group

Question

I am using lme4 package to run linear mixed-effect model. I would like to add the confidence interval of the fitting line per group level in a ggplot.

My data: data is a data frame containing: Plot_label: charactor variable // PD_avg: numeric variable // Year: Factor // GS_Prec: Numeric variable // Direction: Factor

My code as follows:

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

The figure I got here:

To add the CI, I was setting se = T, but it did not work. So I was trying to use geom_ribbon, but it did not work also.

I found one similar topic having the same problem (https://stats.stackexchange.com/questions/552734/r-plotting-lmer-confidence-intervals-per-faceted-group). I did follow the topic, incidentally I got a unexpected result.

My code:

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)

I would like to have the length of the confidence interval should be linked to the range of points. Does anyone know how to represent the CI properly?

This is my subset data

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

@Quinten Hi, I already updated my question with my subset data. Can you check it out? — Anh, Jul 13 '22 at 19:48
See @BenBolker suggestion on [getting confident interval in mixed effect models](https://stackoverflow.com/questions/11072544/how-to-get-coefficients-and-their-confidence-intervals-in-mixed-effects-models). — Adam Quek, Jul 13 '22 at 22:24
@AdamQuek Thank you for very useful link. I tried that but it didnt work for me. I got the CIs for each group, but with that I could not plot manually in `ggplot`. They kept saying that `Error: Aesthetics must be either length 1 or the same as the data (162): ymin and ymax`. — Anh, Jul 14 '22 at 13:54
So, my code to get CI `CI <- as.data.frame(confint(mixed.lm, method="Wald"))`. Then I excluded the `NA` values. Then I did use the bunch of codes to plot: `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)` — Anh, Jul 14 '22 at 13:58

score 1 · Answer 1 · answered Jul 14 '22 at 09:07

You could also use the ggpredict function from the ggeffects package like this:

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)

^{Created on 2022-07-14 by the reprex package (v2.0.1)}

thank you for that. I used the `ggpredict` package and I got the same result too. As you can see the CIs of each fitted line were exaggerated. I mean the CIs should be limitted by 2 poits. Did you think so? — Anh, Jul 14 '22 at 12:43

score 0 · Answer 2 · answered Aug 06 '22 at 12:24

And using the arguments limit.range and facets allows you to limit the predictions to the range of the actual data and creates facts for the plot:

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)

^{Created on 2022-08-06 by the reprex package (v2.0.1)}

Get confidence intervals per faceted group

2 Answers2