Predict with new random effects

Question

When fitting a mixed effects model (or any other model), it's often useful to predict "counterfactual" fitted values, using new data or new model parameters. In my situation, I would like to predict fitted values with new random effects.

The lme4 package allows one to insert new parameters to be used in place of the beta, theta and sigma slots of the merMod class. This would allow one to predict fitted values under counterfactual fixed effects coefficients. How can I do the same for the random effects estimates?

My first thought was to directly modify the u slot of the merMod object, but that doesn't seem to do anything. What can I do?

Example code:

library(lme4)

# use sleepstudy example
fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)

# estimate predictions
fm1Predictions <- predict(fm1)

# estimate predictions with new fixed effects (arbitrarily set to 10)
cfPredictions <- predict(fm1, newparams=list('theta'=10)) # different than fm1Predictions

# estimate predictions with new random effects
fm2 <- fm1
fm2@u <- rep(10,length(fm2@u))
fm2Predictions <- predict(fm2) # same as fm1Predictions

Answer 1

This is not the most elegant solution, but it might work in the meantime.

Using merTools you can return the components of a predicted value from a multilevel model in a data.frame.

library(merTools)
fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
predInt2 <- predictInterval(fm1, which = "all", seed = 8231,
                        include.resid.var = TRUE)

Specifying which = "all" means that the components of the predicted value are returned in addition to the full predicted value.

predInt2[predInt2$obs == 3,]

Returns:

      effect       fit       upr         lwr obs
3   combined 313.16116 353.33704 268.5411738   3
183  Subject  40.02862  83.56716  -0.8552207   3
363    fixed 272.02953 313.85544 230.8799770   3

This data.frame shows you the predicted value for observation 3 from the sleepstudy data broken down by the random component for Subject, the fixed component, and the combined.

You could recompute the combined effect by summing the fixed effect and some counterfactual value for Subject. This is not directly supported by a function in R, but with a little data manipulation you can do:

counter_fact <- predInt2[predInt2$effect == "fixed", ]
counter_fact$est <- counter_fact$fit + 20
head(counter_fact)

And now you have an estimated value if the random component was constant with an effect of 20 for all observations.