Chapter 145 · Advanced & Applied Topics
SEM & Mixed Models · Challenge Solutions
Worked solutions: pooled vs grouped, the ICC, partial-pooling shrinkage, a random slope, and a latent factor from survey items.
In [1]:
import numpy as np, pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import warnings
from statsmodels.tools.sm_exceptions import ConvergenceWarning
# A mixed model with a variance near zero (like our random slope) prints a benign
# "MLE may be on the boundary" convergence note. Silence it for clean output; it does
# not affect the estimates. Remove this line if you want to see the diagnostics.
warnings.simplefilter("ignore", ConvergenceWarning)
plt.rcParams.update({"figure.dpi":110,"axes.grid":True,"grid.alpha":0.25,"font.size":11})
FU, VI, GR = "#a21caf", "#7c3aed", "#16a34a"
BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "structural-equation-modeling-and-mixed-models--student-scores.xlsx"
try:
df = pd.read_excel("../../data/" + fn)
except FileNotFoundError:
df = pd.read_excel(BASE + fn)
print(df.shape); print(df.head())
(1016, 3) school hours_studied test_score 0 SCH100 4.6 57.6 1 SCH100 6.3 74.5 2 SCH100 4.2 73.7 3 SCH100 5.9 65.1 4 SCH100 2.7 73.2
Challenge 1 · Pooled vs grouped¶
OLS ignores schools; the mixed model models them. The slope is similar, but only the mixed model exposes the school structure.
In [2]:
ols = smf.ols("test_score ~ hours_studied", data=df).fit()
mi = smf.mixedlm("test_score ~ hours_studied", df, groups=df["school"]).fit()
print(f"OLS slope = {ols.params['hours_studied']:.3f}")
print(f"Mixed slope = {mi.fe_params['hours_studied']:.3f}")
print("OLS reports one line for everyone; the mixed model adds a per-school intercept.")
OLS slope = 2.369 Mixed slope = 2.310 OLS reports one line for everyone; the mixed model adds a per-school intercept.
Challenge 2 · Compute the ICC¶
ICC = between-school variance / (between + within).
In [3]:
gv, rv = float(mi.cov_re.iloc[0,0]), float(mi.scale)
print(f"between-school var = {gv:.1f}, within var = {rv:.1f}")
print(f"ICC = {gv/(gv+rv):.3f} ({gv/(gv+rv):.0%} of leftover variance is between schools)")
between-school var = 29.0, within var = 66.9 ICC = 0.303 (30% of leftover variance is between schools)
Challenge 3 · See the shrinkage¶
Compare each school's own deviation (its mean residual after the fixed hours effect) with the model's shrunk deviation. Partial pooling pulls every school toward zero.
In [4]:
fixed = mi.fe_params["Intercept"] + mi.fe_params["hours_studied"] * df.hours_studied
raw = (df.test_score - fixed).groupby(df.school).mean() # own deviation, after the hours effect
re = mi.random_effects
shrunk = pd.Series({g: float(re[g].iloc[0]) for g in raw.index})
comp = pd.DataFrame({"own_dev": raw, "shrunk_dev": shrunk}).round(2)
print(comp.head(8))
print(f"\nmean |own| = {raw.abs().mean():.2f} -> mean |shrunk| = {shrunk.abs().mean():.2f}")
print(f"schools pulled toward zero: {(shrunk.abs() < raw.abs()).sum()} of {len(raw)}")
own_dev shrunk_dev SCH100 -0.65 -0.59 SCH101 2.77 2.50 SCH102 3.72 3.43 SCH103 -4.60 -4.28 SCH104 5.09 4.75 SCH105 1.43 1.29 SCH106 4.50 4.16 SCH107 2.79 2.59 mean |own| = 4.63 -> mean |shrunk| = 4.24 schools pulled toward zero: 40 of 40
Challenge 4 · Add a random slope¶
Refit with a random slope and inspect its variance.
In [5]:
ms = smf.mixedlm("test_score ~ hours_studied", df, groups=df["school"], re_formula="~hours_studied").fit()
print(f"slope variance across schools = {ms.cov_re.iloc[1,1]:.3f}")
print("Near zero: the study effect is essentially constant across schools; only baselines differ.")
slope variance across schools = 0.008 Near zero: the study effect is essentially constant across schools; only baselines differ.
Challenge 5 · A latent factor¶
Extract one latent factor from six correlated items and read the loadings.
In [6]:
from sklearn.decomposition import FactorAnalysis
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(7); N = 600
mot = rng.normal(0,1,N)
items = pd.DataFrame({f"item{k+1}": L*mot + np.sqrt(1-L**2)*rng.normal(0,1,N)
for k,L in enumerate([0.8,0.75,0.7,0.65,0.6,0.55])})
fa = FactorAnalysis(n_components=1, random_state=0).fit(StandardScaler().fit_transform(items))
load = np.round(fa.components_[0], 2)
for name, L in zip(items.columns, load): print(f"{name}: loading {L}")
print("\nItems 1-3 load highest, so they measure the trait best.")
item1: loading 0.78 item2: loading 0.7 item3: loading 0.7 item4: loading 0.61 item5: loading 0.58 item6: loading 0.49 Items 1-3 load highest, so they measure the trait best.