Plotting proper mixed models regression slope

Question

I have a data set to which I am fitting a mixed models regression with lme4.

dat <- structure(list(dv280 = c(41L, 68L, 0L, 6L, 20L, 30L, 8L, 1L, 
15L, NA, 59L, 5L, 21L, 41L, 11L, 14L, -2L, 20L, 25L, 33L, 32L, 
30L, 68L, 16L, 11L, -1L, 8L, 0L), v0 = c(55L, 90L, 30L, 23L, 
74L, 48L, 25L, 25L, 46L, NA, 60L, 69L, 55L, 41L, 34L, 41L, 53L, 
76L, 72L, 64L, 34L, 37L, 75L, 21L, 26L, 14L, 24L, 19L), treatment = structure(c(2L, 
1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L), .Label = c("hc", 
"nhc"), class = "factor"), cse = structure(c(2, 2, 6, 6, -4, 
-4, 5, 5, NA, NA, -4, -4, -3, -3, -2, -2, 3, 3, 2, 2, -4, -4, 
-7, -7, 4, 4, 2, 2), .Dim = c(28L, 1L)), pp = structure(c(1L, 
1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 
9L, 10L, 10L, 11L, 11L, 12L, 12L, 13L, 13L, 14L, 14L), .Label = c("j.1", 
"j.3", "j.6", "j.11", "j.13", "j.16", "j.17", "j.18", "j.19", 
"j.22", "j.24", "j.30", "j.32", "j.36"), class = "factor")), .Names = c("dv280", 
"v0", "treatment", "cse", "pp"), row.names = c(NA, 28L), class = "data.frame")
head(dat)
require(lme4)
m <- lmer(dv280 ~ 1 + v0:treatment + cse + (0 + v0 | pp), data=dat, REML=TRUE)
summary(m)

Next I would like to plot the resulting fit, so using the example code from http://glmm.wikidot.com/faq I generated predictions for the original values.

newdat <- data.frame(
  v0=dat$v0,
  treatment=dat$treatment, 
  cse=dat$cse,
  dv280=0)
newdat <- newdat[-c(9,10),]
mm <- model.matrix(terms(m), newdat)
newdat$dv280 <- mm %*% fixef(m)
pvar1 <- diag(mm %*% tcrossprod(vcov(m), mm))
tvar1 <- pvar1 + VarCorr(m)$pp[1]
newdat <- data.frame(newdat, plo=newdat$dv280 - 2 * sqrt(pvar1), 
                     phi=newdat$dv280 + 2 * sqrt(pvar1), 
                     tlo=newdat$dv280 - 2 * sqrt(tvar1), 
                     thi=newdat$dv280 + 2 * sqrt(tvar1)) 

This is easy to plot with ggplot:

p <- ggplot(data=newdat, mapping=aes(x=v0, y=dv280, colour=treatment)) +
  geom_point() +
  geom_smooth(method='lm', se=TRUE) +
  scale_colour_discrete(guide=guide_legend(title.position='left', title.hjust=1))
p + .mytheme + coord_cartesian(xlim=c(-20,100)) +     
  geom_hline(yintercept=0, colour='gray35', linetype='dashed') +
  geom_vline(xintercept=0, colour='gray35', linetype='dashed') 

But as the image shows, the intercepts of the regression lines plotted by ggplot are off compared to the estimated intercept of the regression; at least this is the case for the nhc line which clearly goes to a negative intercept whereas the common estimated intercept is 0.54.

My mistake is easy (I presume): using geom_smooth results in a different fit than what the regression said given that geom_smooth takes the data separately per treatment and also at face value without the context of the fitted model. But I am at a loss how to get the lines properly plotted.

Answer 1

You need to consider the effect of cse. If I set cse=0 in newdat, the lines go through the intercept:

p <- ggplot(data=newdat, mapping=aes(x=v0, y=dv280, colour=treatment)) +
  geom_point() +
  geom_smooth(method='lm', se=TRUE, fullrange=TRUE) +
  scale_colour_discrete(guide=guide_legend(title.position='left', title.hjust=1))
p + theme_bw() + coord_cartesian(xlim=c(-1,1), ylim=c(0.4,0.65)) +     
  geom_hline(yintercept=0, colour='gray35', linetype='dashed') +
  geom_vline(xintercept=0, colour='gray35', linetype='dashed') 

Note that it is normally not advisable to include an interaction in the model without the corresponding main effects.

Edit:

According to your comments you might want something like this (with your original newdat):

p <- ggplot(data=newdat, mapping=aes(x=v0, y=dv280, ymin=tlo, ymax=thi, colour=treatment, fill=treatment)) +
  geom_ribbon(alpha=0.2, aes(colour=NULL)) +
  geom_point() +
  geom_line() +
  scale_colour_discrete(guide=guide_legend(title.position='left', title.hjust=1))
p + theme_bw() +   
  geom_hline(yintercept=0, colour='gray35', linetype='dashed') +
  geom_vline(xintercept=0, colour='gray35', linetype='dashed') 

However, this plot doesn't seem very useful (in particular the lines between points) as it only depicts the projection on the dv280-v0-plane and you can't even see the linear relationship this way.