Including random effects in prediction with Linear Mixed Model

I have a question regarding Linear Mixed Modeling using statsmodels.

The first picture below shows the mixed model I fitted. My dummy dataset only contains one variable, and multiple groups. I would like to predict using not only the fixed intercept and coefficient (see results Intercept and A), but also the random effects of the groups. It seems like it is possible in R, but not using Python (https://www.rdocumentation.org/packages/lme4/versions/1.1-23/topics/predict.merMod)

The second picture shows the predict functionality and some dummy data. As you can see, the prediction is a result of the a combination of only the fixed effects.

Or would it be weird to fit a linear mixed model and also use the random effects in prediction?

Thanks in advance!

3 Answers

There are situations where it would make sense to include the predicted random effects (BLUPs) in a prediction. If you have observations on a group (e.g. a person) and you want to make predictions about additional responses from the same group, then you would generally want to include the BLUP.

This is possible in statsmodels, as illustrated below. We should add an option to make this easier. For now, you have to do it manually, as shown below.

import pandas as pd
import statsmodels.api as sm
import numpy as np

# Simulate data for illustration
group_size = 5
n_groups = 100
x1 = np.random.normal(size=group_size*n_groups)
x2 = np.random.normal(size=group_size*n_groups)
u = np.kron(np.random.normal(size=n_groups), np.ones(group_size))
g = np.kron(np.arange(n_groups), np.ones(group_size))
e = np.random.normal(size=group_size*n_groups)
y = x1 - x2 + u + e

# Fit a multilevel model
df = pd.DataFrame({"y":y, "x1":x1, "x2":x2, "g":g})
model = sm.MixedLM.from_formula("y ~ x1 + x2", groups="g", data=df)
result = model.fit()

# The BLUPs
re = result.random_effects

# Multiply each BLUP by the random effects design matrix for one group
rex = [np.dot(model.exog_re_li[j], re[k]) for (j, k) in enumerate(model.group_labels)]

# Add the fixed and random terms to get the overall prediction
rex = np.concatenate(rex)
yp = result.fittedvalues + rex

• Thank you! We have implemented this approach. However, we cannot figure out what exog_re_li exactly does. Do you have any source on this? Also, we noticed that the current approach shows in sample predictions. Would it be possible to use this for out of sample prediction also, in case you predict for a group that was available during training? And would it still make sense to include this? May 27, 2020 at 10:33
• exog_re_li is a list containing the design matrices for the random effects. The responses for group i have an additive term of the form dot(exog_re_li[i], u), where u is distributed according to the random effects distribution. I think exog_re_li[i] is often called Z_i in the literature. May 29, 2020 at 3:40
• The only reason you can predict the random effect for one block is because you have partial information about that block. If you don't have any data at all from a block, then you should predict its values using the fixed effects only. A typical example would be if you were predicting current blood pressure based on each subject's age, and also blood pressure values recorded in the past. If you want to use age and past blood pressure values for prediction, you should use the fixed effects and the BLUP. If you only want to use age, you should predict based on fixed effects only. May 29, 2020 at 3:49
• It is true that marginally, the random effects only relate to the variance and covariance of your responses, not to the mean values. However if you have access to some data on each block, then using the covariance implied by the random effects allows you to do a better job by predicting conditionally on the partially-observed data, compared to what you would by predicting using the fixed effects alone. If you don't want to use this partial information when making predictions, then most of the advantages of using mixed models are lost. May 29, 2020 at 3:56
• How would I get the prediction for the random effect but in R if may I ask ? Mar 14 at 17:05

To make predictions purely on fixed effects, you can do

md.predict(mdf.fe_params, exog=random_df)


To make predictions on random effects, you can just change the parameters with specifying the particular group name (e.g. "1.5")

md.predict(mdf.random_effects["1.5"], exog=random_df).


I assume the order of features in the test data should follow the same order as what you give as the model's parameters. You can add them together to get the prediction specifically for a group.

I think the simple answer is "yes it would be weird to use random effects for prediction." Random effects in a linear model only effect the variance of your model, not the coefficients. So random effects will change the size of your confidence intervals etc, but not predictions. So if you include random effects in your model, it may change the significance of certain variables, among other things, but it won't change the prediction values.