[R-sig-ME] Plotting lmer estimates with individual points from residuals

[Andy Fugard]

Dear all,

Are there any good reasons not to graph data points from

   fixed effects estimates + observation-level residuals

in models with random effects at levels above residuals?

Some code below to demonstrate what I mean.

I can see arguments against doing this when there are varying slopes in 
the model, but perhaps a separate graph can be used nearby to show how 
the slopes vary between (e.g.) subjects.

Interesting for the sleepstudy dataset that plotting the varying 
intercept (only) model shows a bow-tie effect, pointing to a need for 
random slopes.

Words of wisdom and pointers to literature very welcome.

Best wishes,

Andy


######################################################################

require(lme4)

# Varying intercept model
M1 = lmer(Reaction ~ Days + (1|Subject), sleepstudy)

# Varying slope model
M2 = lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

# Grab the fixed effect estimates
cs1 = fixef(M1)
cs2 = fixef(M1)

days = sleepstudy$Days

# Save residuals and predictions from fixed effects
sleepstudy$resid1     = resid(M1)
sleepstudy$fixedPred1 = with(sleepstudy, cs1[1] + cs1[2]*Days)
sleepstudy$resid2     = resid(M2)
sleepstudy$fixedPred2 = with(sleepstudy, cs2[1] + cs2[2]*Days)

# Plot
par(mfrow = c(1,3))
   plot(Reaction ~ Days, sleepstudy, main = "Original data")
   plot(I(fixedPred1+resid1) ~ Days, sleepstudy,
       main = "Varying intercept", ylab = "Reaction")
   curve(cs1[1] + cs1[2]*x, min(days), max(days), add=T)
   plot(I(fixedPred2+resid2) ~ Days, sleepstudy,
       main = "Varying slope + intercept", ylab = "Reaction")
   curve(cs2[1] + cs2[2]*x, min(days), max(days), add=T)
par(mfrow = c(1,1))

######################################################################


-- 
Andy Fugard, Post-doc, ESF LogICCC (LcpR) project
Fachbereich Psychologie, Universitaet Salzburg
   Hellbrunnerstr. 34, 5020 Salzburg, Austria
+43 (0)680 2199 346  http://figuraleffect.googlepages.com