R: Extract and plot confidence intervals from a lmer object using ggplot [closed]

Question

I am trying to use lmer function from lme4 package to estimate differences between two response curves from a control and treatment responses over time, leaving Subjects as random effect. Here the data:

> df
   Day Subject    Levels   Response
1   10    A001   Control 0.19672131
2   10    A002 Treatment 0.16830515
3   10    A003   Control 0.21355398
4   10    A004   Control 0.18644068
5   10    A005 Treatment 0.17231538
6   10    A007 Treatment 0.18448729
7   11    A001   Control 0.23774081
8   11    A002 Treatment 0.25000000
9   11    A003   Control 0.17288616
10  11    A004   Control 0.25843209
11  11    A005 Treatment 0.29505507
12  11    A007 Treatment 0.27315358
13  12    A001   Control 0.37851189
14  12    A002 Treatment 0.39753941
15  12    A003   Control 0.30925738
16  12    A004   Control 0.45247148
17  12    A005 Treatment 0.37485050
18  12    A007 Treatment 0.41668477
19  13    A001   Control 0.47589286
20  13    A002 Treatment 0.48965316
21  13    A003   Control 0.46696617
22  13    A004   Control 0.50611299
23  13    A005 Treatment 0.41968785
24  13    A007 Treatment 0.51708049
25  14    A001   Control 0.58793970
26  14    A002 Treatment 0.45247189
27  14    A003   Control 0.43121189
28  14    A004   Control 0.56663276
29  14    A005 Treatment 0.37929057
30  14    A007 Treatment 0.46441606
31  15    A001   Control 0.44310684
32  15    A002 Treatment 0.38066676
33  15    A003   Control 0.32576304
34  15    A004   Control 0.39422772
35  15    A005 Treatment 0.28628568
36  15    A007 Treatment 0.34023209
37  16    A001   Control 0.25967359
38  16    A002 Treatment 0.20789686
39  16    A003   Control 0.23629368
40  16    A004   Control 0.22833444
41  16    A005 Treatment 0.24163539
42  16    A007 Treatment 0.21100646
43  17    A001   Control 0.17009653
44  17    A002 Treatment 0.13781610
45  17    A003   Control 0.19149637
46  17    A004   Control 0.21317316
47  17    A005 Treatment 0.17746651
48  17    A007 Treatment 0.15096285
49  18    A001   Control 0.15408115
50  18    A002 Treatment 0.16038546
51  18    A003   Control 0.18361628
52  18    A004   Control 0.18867523
53  18    A005 Treatment 0.20131984
54  18    A007 Treatment 0.19504027
55  19    A001   Control 0.21285064
56  19    A002 Treatment 0.19435679
57  19    A003   Control 0.23979739
58  19    A004   Control 0.24010952
59  19    A005 Treatment 0.20209201
60  19    A007 Treatment 0.25806452
61  20    A001   Control 0.23613019
62  20    A002 Treatment 0.20014232
63  20    A003   Control 0.26122983
64  20    A004   Control 0.26375544
65  20    A005 Treatment 0.17656201
66  20    A007 Treatment 0.22391777
67  21    A001   Control 0.20523904
68  21    A002 Treatment 0.18967355
69  21    A003   Control 0.22878808
70  21    A004   Control 0.26186233
71  21    A005 Treatment 0.18644467
72  21    A007 Treatment 0.18347698
73  22    A001   Control 0.19849361
74  22    A002 Treatment 0.16430202
75  22    A003   Control 0.23331322
76  22    A004   Control 0.25791045
77  22    A005 Treatment 0.18159936
78  22    A007 Treatment 0.17076203
79  23    A001   Control 0.17558492
80  23    A002 Treatment 0.12551814
81  23    A003   Control 0.21406131
82  23    A004   Control 0.22028128
83  23    A005 Treatment 0.17529323
84  23    A007 Treatment 0.14576150
85  24    A001   Control 0.15733775
86  24    A002 Treatment 0.12099877
87  24    A003   Control 0.22833499
88  24    A004   Control 0.15324628
89  24    A005 Treatment 0.15217124
90  24    A007 Treatment 0.09604689

Now I try to fit a 6th order polynomial with a base model with no categorical variables, one to assess the intercept and one to assess the interaction between terms

library(lme4)

model.base=lmer(Response ~ poly(Day, 6, raw=FALSE)+(Day | Subject), df)
model.1=lmer(Response ~ poly(Day, 6, raw=FALSE)+Levels+(Day | Subject), df)
model.2=lmer(Response ~ poly(Day, 6, raw=FALSE)*Levels+(Day | Subject), df)

Then I use anova function to assess the model improvement

> anova(model.base,model.1,model.2)
refitting model(s) with ML (instead of REML)
Data: df
Models:
model.base: Response ~ poly(Day, 6, raw = FALSE) + (Day | Subject)
model.1: Response ~ poly(Day, 6, raw = FALSE) + Levels + (Day | Subject)
model.2: Response ~ poly(Day, 6, raw = FALSE) * Levels + (Day | Subject)
           Df     AIC     BIC logLik deviance   Chisq Chi Df Pr(>Chisq)   
model.base 11 -302.85 -275.35 162.42  -324.85                             
model.1    12 -309.60 -279.61 166.80  -333.60  8.7579      1   0.003083 **
model.2    18 -313.00 -268.00 174.50  -349.00 15.3978      6   0.017378 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

and now my question is how can I plot the fitted data from the model and their confidence interval around the fitted lines similar to this example in ggplot

ggplot(df, aes(Day, Response, color = Levels)) +
  geom_point()+
  scale_x_continuous(breaks = c(seq(10,26,2)), limits = c(9.5,26.5))+
  stat_smooth(method="lm", se=TRUE, 
              formula=y ~ poly(x, 6, raw=FALSE))

So far I have tried confint, effects and lsmeans packages to extract the confidence intervals, being unsuccessful.

Do you have any idea how this could be done?

Answer 1

The plot shows the prediction interval when only varying poly(x, 6, raw = FALSE). If you want to replicate this using lme4 you can use the following (I use the sleepstudy data for reproducabilaty).

library(lme4)
library(merTools)
library(ggplot2)
library(dplyr)

df_new <- sleepstudy %>% filter(Subject %in% c(309, 332))
fit <- lmer(Reaction ~ Days + ( 1 | Subject), sleepstudy)
pred <- cbind(df_new, predictInterval(fit, df_new))

ggplot(pred) + 
  geom_line(aes(Days, fit, fill = Subject)) +
  geom_ribbon(aes(Days, ymin = lwr, ymax = upr, fill = Subject), alpha = .2) +
  geom_point(aes(Days, y = Reaction, color = Subject))

Note that these datapoints are the same as that was trained on. You could create new datapoints by altering only the Days (or in your case Day) variable.