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

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.

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)

grand = df.test_score.mean()
raw = df.groupby("school")["test_score"].mean() - grand
re = mi.random_effects
shrunk = pd.Series({g: float(re[g].iloc[0]) for g in raw.index})
comp = pd.DataFrame({"raw_dev": raw, "shrunk_dev": shrunk}).round(2)
print(comp.head(8))
print(f"\nmean |raw| = {raw.abs().mean():.2f}  ->  mean |shrunk| = {shrunk.abs().mean():.2f}")

        raw_dev  shrunk_dev
SCH100    -1.80       -0.59
SCH101     2.35        2.50
SCH102     4.22        3.43
SCH103    -4.51       -4.28
SCH104     6.08        4.75
SCH105     2.12        1.29
SCH106     6.04        4.16
SCH107     2.73        2.59

mean |raw| = 4.83  ->  mean |shrunk| = 4.24

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.

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.

Challenge 1 · Pooled vs grouped¶

Challenge 2 · Compute the ICC¶

Challenge 3 · See the shrinkage¶

Challenge 4 · Add a random slope¶

Challenge 5 · A latent factor¶