Some plots to illustrate the differences between predict vs simulate in Lme4

predict_vs_simulate.R

library(lme4)

head(sleepstudy)
summary(sleepstudy)
sapply(sleepstudy, class)

# Fit model using the original data:
d <- sleepstudy
d$Subject <- factor(rep(1:18, each=10))
fm1 <- lmer(Reaction ~ Days + (Days|Subject), d)
summary(fm1)

# We need this to avoid a bug in simualte present in older versions of
# lme4:
fm1.parameters <- list(
    beta = fixef(fm1),
    theta = getME(fm1, "theta"),
    sigma = getME(fm1, "sigma"))

# Add 18 new subjects:
d <- rbind(sleepstudy, sleepstudy)
d$Subject <- factor(rep(1:36, each=10))
d$Reaction <- ifelse(d$Subject %in% 19:36, NA, d$Reaction)

# Predict and simulate data for old and new subjects:
d$predicted <- predict (fm1, newdata=d, allow.new.levels=T)
# Drop Reaction to avoid bug in simulate.merMod and use newparam:
d$simulated <- simulate(~ Days + (Days|Subject), seed=1,
                        newdata=d[-1], family=gaussian,
                        newparams=fm1.parameters,
                        allow.new.levels=T)$sim_1

# Plot subject means, simulated data:

f <- function(x, ...) {
    plot(x, xlab="Subject", ylim=c(0, 500), ...)
    grid()
}
par(mfrow=c(1,3), mar=c(4,4,1,1))
with(d, f(tapply(Reaction,  Subject, mean), main="Original data", ylab="Reaction", xlim=c(1, 36)))
with(d, f(tapply(predicted, Subject, mean), main="Predicted data", ylab="", col=rep(1:2, each=18)))
with(d, f(tapply(simulated, Subject, mean), main="Simulated data", ylab="", col=rep(1:2, each=18)))
legend("bottomright", pch=c(1,1), col=1:2, c("old subjects", "new subjects"), bg="white")

# Plot effect of Day per subject:

x <- with(d, as.matrix(tapply(Reaction,  list(Subject, Days), mean), ncol=9))
y <- with(d, as.matrix(tapply(predicted, list(Subject, Days), mean), ncol=9))
z <- with(d, as.matrix(tapply(simulated, list(Subject, Days), mean), ncol=9))

f <- function(...) {
    plot(c(1, 10), c(0, 500), ylab="Raction", xlab="Day", t="n", ...)
    grid()
}
par(mfrow=c(2,3), mar=c(4,4,1,1))

f(main="Day effect, original data")
for (i in 1:18)
  lines(x[i,])
f(main="Day effect, predicted data old subjects")
for (i in 1:18)
  lines(y[i,])
f(main="Day effect, predicted data new subjects")
for (i in 19:36)
  lines(y[i,], col="red")
f(main="Day effect, original data")
for (i in 1:18)
  lines(x[i,])
f(main="Day effect, simulated data old subjects")
for (i in 1:18)
  lines(z[i,])
f(main="Day effect, simulated data new subjects")
for (i in 19:36)
  lines(z[i,], col="red")
legend("bottomright", lty=c(1,1), col=1:2, c("old subjects", "new subjects"), bg="white")

# Lmer simulated data, old subjects vs new subjects:

fm1 <- lmer(Reaction  ~ Days + (Days | Subject), subset(d, Subject%in%1:18))
fm2 <- lmer(simulated ~ Days + (Days | Subject), subset(d, Subject%in%1:18))
fm3 <- lmer(simulated ~ Days + (Days | Subject), subset(d, Subject%in%19:36))

summary(fm1)
summary(fm2)
summary(fm3)