Mixed effects and multilevel models in Python: random intercepts and slopes

Data often come in groups: pupils in schools, patients in hospitals, repeated measurements on the same person. Observations in the same group share influences the data do not record, so they are not independent, and the effect of a regressor may itself differ from group to group. A linear mixed model describes both. Fixed effects give the average relationship; random effects describe how each group departs from it.

What a mixed model estimates

With a random intercept and a random slope on zᵢⱼ, observation i in group j follows

yᵢⱼ = xᵢⱼ′β + u₀ⱼ + u₁ⱼ zᵢⱼ + eᵢⱼ
(u₀ⱼ, u₁ⱼ) ~ N(0, G),   eᵢⱼ ~ N(0, σₑ²)
ICC = var(u₀) / (var(u₀) + σₑ²)   in a random-intercept model

Laird and Ware (1982) set out this general form for longitudinal data. It has two kinds of parameters.

• Fixed effects β describe the average relationship, for example the average gain in test score from one unit of socioeconomic status. Regressors that are constant within a group, such as whether a school is private, can be estimated, which is not possible with fixed group effects.
• Variance components describe the groups. var(u₀) is the variance of the group intercepts, var(u₁) the variance of the group slopes, and their covariance says whether groups with high intercepts also have steep slopes. In a random-intercept model, the intraclass correlation (ICC) is the share of the variance that lies between groups, which is also the correlation between two observations from the same group.

ML or REML. Maximum likelihood estimates of the variance components ignore the degrees of freedom used up by β and are biased downward in small samples. Restricted maximum likelihood (REML) maximizes the likelihood of residual contrasts that do not depend on β, which takes those degrees of freedom into account (Harville 1977). REML likelihoods of models with different fixed effects are not comparable, so likelihood-ratio tests of fixed effects use ML.

Predicted group effects. The best linear unbiased predictors (BLUPs) of uⱼ combine a group's own data with the estimated distribution of the effects. They are shrunk towards zero, more so for small groups and for effects the data pin down poorly (Robinson 1991).

Testing a variance. The hypothesis that a variance is zero lies on the boundary of the parameter space, so the usual chi-square reference distribution does not apply. With one variance component, the likelihood-ratio statistic follows an equal mixture of χ²(0) and χ²(1) (Self and Liang 1987). With several, a chi-square test with as many degrees of freedom as variance parameters is conservative.

A random-intercept model assumes that the group effects are uncorrelated with the regressors, the assumption of random effects in panel regressions. If that is doubtful, the within estimator is the safer choice for regressors that vary within groups.

In OpenEconometrics

oe.mixed fits linear mixed models with one grouping level, or two nested ones. The model in the example corresponds to the Stata command mixed score ses private || school: ses, covariance(unstructured) reml:

oe.mixed(data=df, y="score", x=["ses", "private"], group="school", random=["ses"], covstructure="unstructured", method="reml")

• group="school", or group=["school", "class"] with the top level first, which corresponds to || school: || class:. Class labels are read within their school.
• random= lists the columns with random slopes; random_intercept=False drops the random intercept. covstructure is "independent" (the default), "unstructured", "exchangeable" or "identity".
• method="ml" (the default) or "reml".
• covariance="nonrobust" (the default) uses the model-based covariance; "robust" clusters on the top-level groups, and "cluster" on a column that nests them. Both need method="ml".
• Fixed effects have z tests. Variance components are reported on the variance scale, as terms such as /var(_cons[school]); the intervals of variances are formed on the log scale. tests["lr_vs_linear"] compares the model with a linear regression, and metrics["icc"] reports the ICC for random-intercept models.
• oe.mixed_predict(fit, df, kind="reffects") returns the BLUPs; kind="xb" and kind="fitted" give the fixed-part and the group-conditional predictions.

Related functions cover other needs: oe.mixed_contrast gives Kenward–Roger (REML fits) and Satterthwaite t and F tests of fixed effects, oe.mixedresid adds AR(1) or exchangeable residual correlation, and oe.melogit, oe.meprobit and oe.mepoisson fit binary and count outcomes. The mixed-model guide gives the likelihoods, the parameterization and the covariance rules.

Supported options

Options of oe.mixed

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 120 schools with 10 to 40 pupils each. A pupil's score is 50 + 4·ses + 3·private, plus a school intercept effect (standard deviation 3), a school-specific ses slope (standard deviation 1.5, correlation 0.4 with the intercept effect) and a pupil error (standard deviation 6). The true variance components are therefore var(_cons) = 9, var(ses) = 2.25, cov = 1.8 and var(Residual) = 36.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(schools=120, seed=21):
    """Test scores of pupils in schools, with a random intercept and a random ses slope per school."""
    rng = np.random.default_rng(seed)
    sd, corr = np.array([3.0, 1.5]), 0.4                      # sd of intercept and slope effects
    cov = np.outer(sd, sd) * np.array([[1, corr], [corr, 1]])
    effects = rng.multivariate_normal([0, 0], cov, size=schools)
    private = (rng.random(schools) < 0.3).astype(float)      # a school-level regressor
    size = rng.integers(10, 41, size=schools)                # 10 to 40 pupils per school
    school = np.repeat(np.arange(schools), size)
    ses = rng.normal(size=school.size)
    score = (50 + 4 * ses + 3 * private[school] + effects[school, 0] + effects[school, 1] * ses
             + rng.normal(scale=6, size=school.size))
    df = pd.DataFrame({"school": school, "ses": ses, "private": private[school], "score": score})
    truth = pd.DataFrame(effects, columns=["_cons", "ses"])
    return df, truth


if __name__ == "__main__":
    df, truth = simulate()
    intercept_only = oe.mixed(data=df, y="score", x=["ses", "private"], group="school")
    print(f"Intraclass correlation, random-intercept model: {intercept_only.metrics['icc']:.3f}")
    print()
    fit = oe.mixed(data=df, y="score", x=["ses", "private"], group="school",
                   random=["ses"], covstructure="unstructured", method="reml")
    print(fit.summary())
    print()
    blups = oe.mixed_predict(fit, df, kind="reffects")
    blups = blups.pivot(index="group", columns="effect", values="blup").sort_index()
    for effect in ["_cons", "ses"]:
        r = np.corrcoef(blups[effect], truth[effect])[0, 1]
        print(f"BLUP of {effect:<5}: sd {blups[effect].std():.3f} (true effects {truth[effect].std():.3f}), "
              f"correlation with the true effects {r:.3f}")

With openecon 0.3.20a1 from PyPI the script prints:

Intraclass correlation, random-intercept model: 0.182

Mixed-effects REML regression — score
Observations: 2840  |  Covariance: nonrobust  |  Confidence: 95%

Term                     Estimate  Std. error        z      P>|stat|  CI lower  CI upper
-----------------------  --------  ----------  -------  ------------  --------  --------
[score]                                                                                 
Intercept                 50.4978    0.347505  145.315             0   49.8167   51.1789
ses                       3.99692    0.168533  23.7159  2.47024e-124    3.6666   4.32724
private                   2.08364    0.637242  3.26977    0.00107633  0.834666   3.33261
[school]                                                                                
/var(ses[school])         1.63393    0.418008  3.90884   9.27409e-05  0.989621   2.69771
/var(_cons[school])       8.67751     1.35585  6.40003   1.55344e-10   6.38845   11.7868
/cov(ses,_cons[school])    1.3049     0.55925   2.3333     0.0196323  0.208789   2.40101
[Residual]                                                                              
/var(Residual)            35.7894    0.990513  36.1322  7.08623e-286   33.8997   37.7844

log_likelihood: -9257.89  |  aic: 18529.8  |  bic: 18571.4  |  n_groups: 120
group_size_min: 10  |  group_size_avg: 23.6667  |  group_size_max: 40
Wald chi2 test of the fixed slopes: chi2(2) = 571.88, p = 6.574e-125
LR test vs. linear model (conservative: variance parameters on the boundary): chi2(3) = 422.988, p = 2.32e-91

BLUP of _cons: sd 2.674 (true effects 2.818), correlation with the true effects 0.914
BLUP of ses  : sd 0.911 (true effects 1.295), correlation with the true effects 0.671

Reading the output

• The fixed effects are recovered. The ses slope is 3.997 (95% interval 3.667 to 4.327) and the intercept 50.50 (49.82 to 51.18). The private-school effect, estimated from only 120 schools, is 2.084 with a wide interval from 0.835 to 3.333 that still contains 3.
• So are the variance components. REML gives var(_cons) = 8.678, var(ses) = 1.634, cov = 1.305 and var(Residual) = 35.79. Each interval contains the true value. The slope variance and the covariance are the least precise, with intervals from 0.990 to 2.698 and from 0.209 to 2.401.
• The groups matter. In the random-intercept model, 18.2% of the variance lies between schools (ICC 0.182). The likelihood-ratio test against a linear regression gives χ²(3) = 423.0; with three variance parameters the test is conservative, and it still rejects decisively.
• The BLUPs are shrunk. The predicted school intercepts have a standard deviation of 2.674 against 2.818 for the true effects, and correlate 0.914 with them. The predicted slopes are shrunk much more (0.911 against 1.295) and correlate only 0.671 with the truth, because 10 to 40 pupils say little about a school's own slope.

Limitations

• Cluster-robust covariance for oe.mixed allows at most 1 clustering dimension.
• At most two nested levels. Random slopes are available at the lowest level only; the top level of a two-level model has a random intercept only.
• oe.mixed has no residual-correlation structures and no small-sample degrees of freedom; use oe.mixedresid and oe.mixed_contrast for those.
• Only frequency weights are accepted; level-specific sampling weights are not offered.
• robust and cluster covariances need method="ml", and no likelihood-ratio test is reported with them.
• A variance component estimated at zero is refused with boundary_solution rather than reported with meaningless standard errors; refit without that component.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358), 320–338. doi:10.1080/01621459.1977.10480998
• Laird, N. M. and Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963–974. doi:10.2307/2529876
• Robinson, G. K. (1991). That BLUP is a good thing: The estimation of random effects. Statistical Science, 6(1), 15–32. doi:10.1214/ss/1177011926
• Self, S. G. and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. doi:10.1080/01621459.1987.10478472