Test how simulate.merMod deals with new factor levels

test_simulate.R

library(lme4)

head(sleepstudy)
summary(sleepstudy)

# Relabel subjects:
d <- sleepstudy
d$Subject <- factor(rep(1:18, each=10))

# Fit model:
fm1 <- lmer(Reaction ~ Days + (Days | Subject), d)
summary(fm1)

# Add new subjects:
d <- data.frame(rbind(sleepstudy, sleepstudy))
d$Subject <- factor(rep(1:36, each=10))

# Simulate new data for old and new subjects:
d$simulated <- simulate(fm1, nsim=1, newdata=d, allow.new.levels=T)$sim_1
# No measurements for new subjects:
d$Reaction  <- ifelse(d$Subject %in% 19:36, NA, d$Reaction)

# Plot subject means:

f <- function(x, ...) {
    plot(x, xlab="Subject", ylim=c(0, 500), ...)
    grid()
}
par(mfrow=c(1,2), 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(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 <- as.matrix(with(d, tapply(Reaction, list(Subject, Days), mean)), ncol=9)
y <- as.matrix(with(d, tapply(rsim,     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(1,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, simulated data old subjects")
for (i in 1:18)
  lines(y[i,])

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