Multilevel (Hierarchical) Linear Models

Multilevel models (also called hierarchical or mixed-effects models) handle data with nested structures — students within schools, patients within hospitals, repeated measures within individuals.

Why Multilevel Models?

| Issue | Consequence | Solution |

|-------|------------|---------|

| Non-independence | Standard errors too small | Random effects |

| Varying intercepts | Different group baselines | Random intercept model |

| Varying slopes | Different group effects | Random slope model |

| Crossed effects | Multiple grouping factors | Crossed random effects |

Random Intercept Model

Random Slope Model

Allows the effect of X to vary across groups:

Intraclass Correlation (ICC)

The ICC measures the proportion of total variance that lies between groups.

| ICC | Interpretation |

|-----|---------------|

| 0.00 - 0.05 | Negligible clustering — standard regression may suffice |

| 0.05 - 0.15 | Moderate clustering — multilevel modeling recommended |

| 0.15 - 0.25 | Substantial clustering — multilevel modeling essential |

| > 0.25 | Very strong clustering — must model group structure |

Variance Partition Coefficient

The VPC (same as ICC in the random intercept model) answers: How much does the outcome vary between groups?

Estimation Methods

| Method | Description | When to Use |

|--------|-------------|-------------|

| REML (Restricted Maximum Likelihood) | Unbiased variance estimates | Default; preferred |

| ML (Maximum Likelihood) | Biased variance estimates | When comparing models with different fixed effects |

| Bayesian | Full posterior distribution | Complex models; small samples |

Model Comparison

Test the significance of random effects by comparing models with and without the random effect.

Python Implementation

import numpy as np

import pandas as pd

import statsmodels.api as sm

from statsmodels.regression.mixed_linear_model import MixedLM

from scipy import stats

import matplotlib.pyplot as plt



np.random.seed(42)



# Simulate hierarchical data: students in schools

n_schools = 30

n_students_per = 20

n = n_schools * n_students_per



school_id = np.repeat(np.arange(n_schools), n_students_per)

school_effect = np.random.randn(n_schools) * 2  # Random intercepts

u_0 = school_effect[school_id]



X = np.random.randn(n)

Y = 5 + u_0 + 0.8 * X + np.random.randn(n) * 1.5



df = pd.DataFrame({'Y': Y, 'X': X, 'school': school_id})



# ICC calculation

from statsmodels.formula.api import ols

null_model = ols('Y ~ 1', data=df).fit()

total_var = np.var(null_model.resid)



# Between-school variance

school_means = df.groupby('school')['Y'].mean()

between_var = np.var(school_means)



# Within-school variance

within_var = total_var - between_var * n_students_per / (n_students_per - 1)

ICC = between_var / (between_var + within_var)

print(f"ICC: {ICC:.3f}")



# Random intercept model

ri_model = MixedLM.from_formula('Y ~ X', groups='school', data=df)

ri_result = ri_model.fit()

print(f"\nRandom Intercept Model:")

print(ri_result.summary())



# Random slope model

rs_model = MixedLM.from_formula('Y ~ X', groups='school', re_formula='~X', data=df)

rs_result = rs_model.fit()

print(f"\nRandom Slope Model:")

print(rs_result.summary())



# Compare models using AIC

print(f"\nRandom Intercept AIC: {ri_result.aic:.1f}")

print(f"Random Slope AIC: {rs_result.aic:.1f}")

Worked Example

Key Takeaways

Related Topics

• See Panel Data Analysis for longitudinal extensions

• See Mixed Effects Models for more complex structures

• See Regression Diagnostics for checking assumptions