The Statsmodels imputation of linear mixed models (MixedLM) closely follows the approach outlined in Lindstrom and Bates (JASA 1988). This is also the approach followed in the R package LME4. Other packages such as Stata, SAS, etc. should also be consistent with this approach, as the basic techniques in this area are mostly mature.

Here we show how linear mixed models can be fit using the MixedLM procedure in Statsmodels. Results from R (LME4) are included for comparison.

Here are our import statements:

In [11]:

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

These are longitudinal data from a factorial experiment. The outcome variable is the weight of each pig, and the only predictor variable we will use here is "time". First we fit a model that expresses the mean weight as a linear function of time, with a random intercept for each pig. The model is specified using formulas. Since the random effects structure is not specified, the default random effects structure (a random intercept for each group) is automatically used.

In [12]:

```
data = pd.read_csv("dietox.csv")
model = sm.MixedLM.from_formula("Weight ~ Time", data, groups=data["Pig"])
result = model.fit()
print result.summary()
```

Here is the same model fit in R using LMER:

```
Linear mixed model fit by REML
Formula: Weight ~ Time + (1 | Pig)
Data: data
AIC BIC logLik deviance REMLdev
4818 4837 -2405 4806 4810
Random effects:
Groups Name Variance Std.Dev.
Pig (Intercept) 40.394 6.3556
Residual 11.367 3.3715
Number of obs: 861, groups: Pig, 72
Fixed effects:
Estimate Std. Error t value
(Intercept) 15.72352 0.78805 19.95
Time 6.94251 0.03339 207.94
Correlation of Fixed Effects:
(Intr)
Time -0.275
```

Note that in the Statsmodels summary of results, the fixed effects and random effects parameter estimates are shown in a single table. The random effects are labeled "RE" in the Statmodels output. In the LME4 output, this effect is the pig intercept under the random effects section.

There has been a lot of debate about whether the standard errors for random effect variance and covariance parameters are useful. In LME4, these standard errors are not displayed, because the authors of the package believe they are not very informative. While there is good reason to question their utility, we elected to include the standard errors in the summary table, but do not show the corresponding Wald confidence intervals.

Next we fit a model with two random effects for each animal: a random intercept, and a random slope (with respect to time). This means that each pig may have a different baseline weight, as well as growing at a different rate. The formula specifies that "Time" is a covariate with a random coefficient. By default, formulas always include an intercept (which could be suppressed here using "0 + Time" as the formula).

In [13]:

```
model = sm.MixedLM.from_formula("Weight ~ Time", data, re_formula="Time", groups=data["Pig"])
result = model.fit()
print result.summary()
```

Here is the same model fit using LMER in R:

```
Linear mixed model fit by REML
Formula: Weight ~ Time + (1 + Time | Pig)
Data: data
AIC BIC logLik deviance REMLdev
4446 4475 -2217 4432 4434
Random effects:
Groups Name Variance Std.Dev. Corr
Pig (Intercept) 19.49346 4.41514
Time 0.41606 0.64503 0.103
Residual 6.03745 2.45712
Number of obs: 861, groups: Pig, 72
Fixed effects:
Estimate Std. Error t value
(Intercept) 15.73865 0.55013 28.61
Time 6.93901 0.07982 86.93
Correlation of Fixed Effects:
(Intr)
Time 0.005
```

In [14]:

```
from statsmodels.regression.mixed_linear_model import MixedLMParams
model = sm.MixedLM.from_formula("Weight ~ Time", data, re_formula="Time", groups=data["Pig"])
free = MixedLMParams(2, 2)
free.set_fe_params(np.ones(2))
free.set_cov_re(np.eye(2))
result = model.fit(free=free)
print result.summary()
```

The likelihood drops by 0.3 when we fix the correlation parameter to 0. Comparing 2 x 0.3 = 0.6 to the chi^2 1 df reference distribution suggests that the data are very consistent with a model in which this parameter is equal to 0.

Here is the same model fit using LMER in R (note that here R is reporting the REML criterion instead of the likelihood, where the REML criterion is twice the log likeihood):

```
Linear mixed model fit by REML ['lmerMod']
Formula: Weight ~ Time + (1 | Pig) + (0 + Time | Pig)
Data: data
REML criterion at convergence: 4434.7
Scaled residuals:
Min 1Q Median 3Q Max
-6.4281 -0.5527 -0.0405 0.4840 3.5661
Random effects:
Groups Name Variance Std.Dev.
Pig (Intercept) 19.8409 4.4543
Pig.1 Time 0.4234 0.6507
Residual 6.0282 2.4552
Number of obs: 861, groups: Pig, 72
Fixed effects:
Estimate Std. Error t value
(Intercept) 15.73875 0.55445 28.39
Time 6.93899 0.08045 86.25
Correlation of Fixed Effects:
(Intr)
Time -0.086
```

This is one of the example data sets provided in the LMER R library. The outcome variable is the size of the tree, and the covariate used here is a time value. The data are grouped by tree.

In [15]:

```
data = pd.read_csv("Sitka.csv")
endog = data["size"]
exog = np.ones((data.shape[0], 2), dtype=np.float64)
exog[:,1] = data["Time"]
exog = pd.DataFrame(data=exog)
exog.columns = ["Intercept", "Time"]
```

In [17]:

```
model = sm.MixedLM(endog, exog, data["tree"], exog.iloc[:,0])
result = model.fit()
print result.summary()
```

Here is the same model fit in R using LMER:

```
Linear mixed model fit by REML
Formula: size ~ Time + (1 | tree)
Data: Sitka
AIC BIC logLik deviance REMLdev
172.8 188.7 -82.39 146.6 164.8
Random effects:
Groups Name Variance Std.Dev.
tree (Intercept) 0.374512 0.61197
Residual 0.039206 0.19800
Number of obs: 395, groups: tree, 79
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.2732443 0.0878948 25.86
Time 0.0126855 0.0002654 47.80
Correlation of Fixed Effects:
(Intr)
Time -0.611
```

```
Linear mixed model fit by REML
Formula: size ~ Time + (1 + Time | tree)
Data: Sitka
AIC BIC logLik deviance REMLdev
176.6 200.4 -82.28 146.4 164.6
Random effects:
Groups Name Variance Std.Dev. Corr
tree (Intercept) 3.4418e-01 0.58666576
Time 1.5640e-08 0.00012506 1.000
Residual 3.9179e-02 0.19793693
Number of obs: 395, groups: tree, 79
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.2732443 0.0856713 26.53
Time 0.0126855 0.0002657 47.75
Correlation of Fixed Effects:
(Intr)
Time -0.585
```

In [18]:

```
exog_re = np.asarray(exog.copy())
model = sm.MixedLM(endog, exog, data["tree"], exog_re)
result = model.fit()
print result.summary()
```

In [ ]:

```
import warnings
warnings.filterwarnings("ignore")
re = result.cov_re.iloc[0, 0]
likev = result.profile_re(0, dist_low=0.1, dist_high=0.1)
```

In [21]:

```
plt.figure(figsize=(10,8))
plt.plot(likev[:,0], 2*likev[:,1])
plt.xlabel("Variance of random slope", size=17)
plt.ylabel("-2 times profile log likelihood", size=17)
```

Out[21]:

In [ ]:

```
re = result.cov_re.iloc[1, 1]
likev = result.profile_re(1, dist_low=0.8*re, dist_high=0.8*re)
plt.figure(figsize=(10, 8))
plt.plot(likev[:,0], 2*likev[:,1])
plt.xlabel("Variance of random slope", size=17)
plt.ylabel("-2 times profile log likelihood", size=17)
```

Out[ ]: