How can I better predict with (g)lmer with missing values? Suppose I'm building a mixed model in R, and I want to use that model to predict new data for which I might not know the value of all the features. Or in some cases, it might not be so much that I don't know the value of all the features but more that I want a prediction for a "neutral" value of those features in some meaningful sense. Is there a "proper" way to do this?
Here's a silly example.
install.packages("lme4")
require(lme4)
mtcars2 = mtcars
mtcars2$cyl = paste0(mtcars2$cyl, " cylinders")
# I want to make sure we don't treat 'cyl' numerically. This
# probably isn't the best approach, but it should work.

model = lmer(mpg ~ (1|cyl) + hp, mtcars2)

With my silly model now built, I can use it to make a prediction about a new row of data.
testExample = data.frame(c("4 cylinders"), c(100))
names(testExample) = c("cyl", "hp")

testPrediction = predict(model, testExample)
print(testPrediction)

This yields a numerical result, as expected of course. But suppose I wanted to know what the model would predict if I only knew the hp value and didn't know the cyl value. 
testExample$cyl = c("New value not seen before")
testPrediction = predict(model, testExample, allow.new.levels=TRUE)
print(testPrediction)

I can add allow.new.levels=TRUE, and this at least lets me get a number (21.66101). I'm not convinced it's the right number, but we'll get back to that in a minute.
It seems natural to me that I should be able to apply the same sort of logic if instead I had a legitimate value for cyl but hp were missing. The natural way for me would seem to be to put NA in place of the numeric value for hp as follows.
testExample$cyl = c("4 cylinders") # restoring the original value
testExample$hp = c(NA)
testPrediction = predict(model, testExample, allow.new.levels=TRUE)
print(testPrediction)

This approach does not yield a number. I'm not sure if the difference between the two examples lies more in the difference between fixed and random effects or in one being numeric and the other not.
So is there some direct way of handling this in R? I think (but I'm not sure) that what I might want is for the mean to be imputed in my numeric case of hp. That doesn't seem as straight-forward in my cyl case, as I can't take the mean of some strings. It seems that what allow.new.levels effectively does in my case is ignore the part of the model for that feature entirely as though it produced a value of 0, and so it just computes the fixed effect intercept plus 100 times the fixed effect slope. Are there standard ways for handling these sorts of things?
 A: A couple of points:


*

*Lets start with a general definition of the linear mixed model, namely, let $y_i$ be the outcome for the $i$-th subject ($i = 1, \ldots, n$), $X_i$ the design matrix for the fixed effects $\beta$, and $Z_i$ the design matrix for the random effects $b_i$, then we have:
$$\left\{
\begin{array}{l}
y_i = X_i \beta + Z_i b_i + \varepsilon_i,\\\\
b_i \sim \mathcal N(0, D), \quad \varepsilon_i \sim \mathcal N(0, \sigma^2 \mbox{I}),
\end{array}
\right.$$
where $\varepsilon_i$ are the error terms. When you fit this model you get estimates for the fixed effects $\hat \beta$. You can also get estimates of the random effects $\hat b_i$ which are defined as the modes of the conditional distribution:
$$
p(b_i \mid y_i; \hat \theta) = \frac{p(y_i \mid b_i; \hat \theta) \, p(b_i; \hat 
\theta)}{p(y_i; \hat \theta)},$$
where $\theta = (\beta, \mbox{vech}(D), \sigma^2)$ denotes all model parameters.

*From this model you can get two types of predictions, namely population predictions defined as $$\hat y_i^{pop} = X_i \hat \beta,$$ and subject-specific predictions, defined as $$\hat y_i^{subj} = X_i \hat \beta + Z_i \hat b_i.$$

*Say now that you want to calculate prediction for a new subject $j$. If you have no outcome data for this subject, then the only prediction you can do is the population prediciton, namely $$\hat y_j^{pop} = X_j \hat \beta,$$ which will be the average for his group. If however you also have some data $y_j^{new}$ for this subject, then you could estimate his random effect $\hat b_j$ as the mode of $[b_j \mid y_j^{new}; \hat \theta]$, and then use this into the above equation of the subject-specific prediction, i.e., $$\hat y_j^{subj} = X_j \hat \beta + Z_j \hat b_j.$$
You can find additional information regarding this in Chapter 3, Section 3.4 of my Repeated Measurements course notes.
