pipeline_miami/analysis/residual_z.py
Served verbatim from the project repository. Internal working document conventions apply: documents may reference file paths, branch names, and findings-ledger anchors from the repo.
"""Residual-z outlier engine for Miami-Dade (docs/design/nyc-school-outcomes-spec.md §3).
For each outcome × year: weighted least squares of the school value on the
same-year covariates (pipeline_miami/config/covariates.yaml), weights = the
cell's sampling N (outcomes.yaml n_rule). The residual is standardised by
the between-school residual spread and the cell's own sampling error:
z = resid / sqrt(tau^2 + se_i^2),
tau^2 = max(0, s^2 - mean_w(se_i^2)), s^2 = sum(w r^2) / sum(w) * n / (n - k),
se_i = binomial (percent) or Poisson (rate) SE, or 0 where the outcome
is not a share of independent students (outcomes.yaml sampling_error:
none — attendance is a share of student-days, review 2026-09-16b N11).
s^2 already contains the sampling variance, so it is removed once rather
than added twice, and the n/(n-k) factor is the degrees-of-freedom
correction (review 2026-09-16, M3). Cells with N below min_n are ineligible.
|z| >= z_threshold flags; the sign is oriented so "positive" always means
doing well. Persistence: >= min_flags of the same sign within the last
`window` calendar years inside one comparability regime (excluded years
leave a gap; flagged years are fitted and marked). The fit universe is every
school with the outcome that year, closed schools included.
Everything below is descriptive, not causal. Reads the served DB read-only;
writes data/cities/miami/reference/residual_z.csv (per cell; loaded into the
metrics table by scripts/loaders/miami/residual_z.ts) and
data/analysis/miami-residual-z.json (fit statistics, the pre-registered
checks, the covariate ablation and the flip-diff against peer percentiles).
Run: python -m pipeline_miami.analysis.residual_z
"""
from __future__ import annotations
import datetime as dt
import hashlib
import json
import math
import re
import numpy as np
import pandas as pd
import yaml
from scipy import stats
from pipeline_miami.verify import common
CONFIG_DIR = common.REPO_ROOT / "pipeline_miami" / "config"
OUT_CSV = common.REPO_ROOT / "data" / "cities" / "miami" / "reference" / "residual_z.csv"
OUT_JSON = common.REPO_ROOT / "data" / "analysis" / "miami-residual-z.json"
GRADES = ["PK", "KG", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"]
# Pre-registered checks, fixed before the first fit (2026-09-15). The residual
# is meant to be the part of achievement not explained by who enrols. Florida's
# learning gains are an independent growth measure computed by the state from
# each student's own prior score, so:
# P1 Spearman rho(z_achievement, learning_gains) > 0 with p < 0.001 in every fitted year;
# P2 rho(z, gains) >= rho(raw achievement, gains) in every fitted year
# (adjusting for demographics should not weaken the link to growth).
# A demographic gradient does not make failure mathematically inevitable.
# These are Spearman correlations of z, not a Pearson decomposition of
# raw = fitted + residual. The fitted/gains and lowest-quartile comparisons
# are descriptive sensitivity checks, not proofs of a mechanism.
# P3 |rho(z, pct_econ_disadv)| < 0.10 (a descriptive diagnostic, not proof of independence);
# P4 fitted values outside [0, 100] for fewer than 2% of eligible cells (bounded-rate guardrail);
# P5 flagged share (|z| >= 2) between 2% and 20% of eligible cells.
PREREGISTERED = {
"P1": {"rule": "rho(z, gains) > 0 and p < 0.001, every fitted year", "applies_to": ["ela_achievement", "math_achievement"]},
"P2": {"rule": "rho(z, gains) >= rho(raw, gains), every fitted year", "applies_to": ["ela_achievement", "math_achievement"],
"caveat": "empirical rank-correlation comparison; a demographic gradient does not predetermine its result, and lowest-quartile gains use the same gain rules"},
"P3": {"rule": "|rho(z, pct_econ_disadv)| < 0.10, every fitted year", "applies_to": "all"},
"P4": {"rule": "share of fitted outside [0,100] < 2%", "applies_to": "percent-scale outcomes"},
"P5": {"rule": "flagged share in [2%, 20%]", "applies_to": "all"},
}
def _clean(o):
"""NaN/inf → null so the packet is valid JSON for the app's static import."""
if isinstance(o, dict):
return {k: _clean(v) for k, v in o.items()}
if isinstance(o, (list, tuple)):
return [_clean(v) for v in o]
if isinstance(o, (float, np.floating)):
return None if not math.isfinite(float(o)) else float(o)
if isinstance(o, (np.integer,)):
return int(o)
if isinstance(o, (np.bool_,)):
return bool(o)
return o
def spring_year(y: str) -> int:
return int(y[:4]) + 1
def load_config():
cov = yaml.safe_load(open(CONFIG_DIR / "covariates.yaml"))
out = yaml.safe_load(open(CONFIG_DIR / "outcomes.yaml"))
comp = yaml.safe_load(open(CONFIG_DIR / "comparability.yaml"))
return cov, out, comp
def regime_windows(comp: dict, regime: str) -> tuple[dict[str, str], set[str], set[str]]:
"""year → window id within a regime, from boundary years only; excluded years dropped; flagged years fitted and marked."""
r = comp["regimes"][regime]
bounds = sorted(str(b["year"]) for b in r.get("boundaries", []))
excluded = set(str(y) for y in r.get("excluded_years", []))
flagged = set(str(y) for y in r.get("flagged_years", []))
def window_of(year: str) -> str:
w = "pre"
for b in bounds:
if year >= b:
w = b
return f"{regime}:{w}"
return window_of, excluded, flagged # type: ignore[return-value]
def parse_n_rule(rule: str):
"""'denominator' | 'enrollment' | 'grades:3-10' | 'grades:5,8,9-12;share:0.25'."""
if rule in ("denominator", "enrollment"):
return {"kind": rule}
m = re.match(r"grades:([0-9A-Z,\-]+)(?:;share:([0-9.]+))?$", rule)
if not m:
raise SystemExit(f"bad n_rule {rule}")
grades: list[str] = []
for part in m.group(1).split(","):
if "-" in part:
a, b = part.split("-")
grades += [str(g) for g in range(int(a), int(b) + 1)]
else:
grades.append(part)
return {"kind": "grades", "grades": grades, "share": float(m.group(2)) if m.group(2) else 1.0}
def fetch_frame(conn, cov_cfg, outcome_keys: list[str]) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
num_keys = [c["key"] for c in cov_cfg["numeric"]]
keys = sorted(set(num_keys + outcome_keys + ["enrollment_oct", "pct_econ_disadv", "pct_ell_upper_bound", "learning_gains_ela", "learning_gains_math", "learning_gains_ela_l25", "learning_gains_math_l25"] + [f"enrollment_oct_g{g}" for g in GRADES]))
q = """
SELECT m.school_number, m.year, m.metric_key, m.value::float8, m.denominator::float8, m.suppressed,
m.comparison_group_percentile::float8, m.citywide_percentile::float8
FROM miami_school_year_metrics m
WHERE m.subgroup = 'ALL' AND m.metric_key = ANY(%s)
"""
rows = conn.execute(q, (keys,)).fetchall()
df = pd.DataFrame(rows, columns=["school", "year", "key", "value", "den", "suppressed", "peer_pct", "city_pct"])
sch = pd.DataFrame(conn.execute("SELECT school_number, school_level, is_charter, magnet_status, activity_code, closed_at, include_in_default_comparisons, name FROM miami_schools").fetchall(),
columns=["school", "level", "is_charter", "magnet_status", "activity_code", "closed_at", "include", "name"])
sch["sector"] = np.where(sch.is_charter, "charter", "district")
sch["magnet"] = np.where(sch.magnet_status == "S", "magnet", "other") # S = whole-school magnet (MSID code table)
vals = df[~df.suppressed & df.value.notna()].pivot_table(index=["school", "year"], columns="key", values="value", aggfunc="first")
dens = df[df.key.isin(outcome_keys)].pivot_table(index=["school", "year"], columns="key", values="den", aggfunc="first")
dens.columns = [f"den__{c}" for c in dens.columns]
pcts = df[df.key.isin(outcome_keys)].pivot_table(index=["school", "year"], columns="key", values="peer_pct", aggfunc="first")
pcts.columns = [f"peer__{c}" for c in pcts.columns]
frame = vals.join(dens, how="left").join(pcts, how="left").reset_index().merge(sch, on="school", how="left")
return frame, sch, df
def cell_n(frame: pd.DataFrame, key: str, rule: dict) -> pd.Series:
if rule["kind"] == "denominator":
return frame.get(f"den__{key}", pd.Series(np.nan, index=frame.index))
if rule["kind"] == "enrollment":
return frame.get("enrollment_oct", pd.Series(np.nan, index=frame.index))
cols = [f"enrollment_oct_g{g}" for g in rule["grades"] if f"enrollment_oct_g{g}" in frame.columns]
n = frame[cols].sum(axis=1, min_count=1) * rule["share"]
return n
def design(df: pd.DataFrame, cov_cfg, drop: set[str] = frozenset()) -> tuple[np.ndarray, list[str]]:
cols = [np.ones((len(df), 1))]
names = ["intercept"]
for c in cov_cfg["numeric"]:
if c["key"] in drop:
continue
x = df[c["key"]].to_numpy(dtype=float)
if c.get("transform") == "log":
x = np.log(np.maximum(x, 1.0))
names.append(f"log_{c['key']}")
else:
names.append(c["key"])
cols.append(x[:, None])
for c in cov_cfg["categorical"]:
if c["key"] in drop:
continue
d = pd.get_dummies(df[c["key"]].astype("category"), prefix=c["key"], drop_first=True)
d = d.reindex(sorted(d.columns), axis=1)
if d.shape[1]:
cols.append(d.to_numpy(dtype=float))
names += list(d.columns)
return np.hstack(cols), names
def fit_wls(X: np.ndarray, y: np.ndarray, w: np.ndarray) -> dict:
sw = np.sqrt(w)
beta, *_ = np.linalg.lstsq(X * sw[:, None], y * sw, rcond=None)
fitted = X @ beta
resid = y - fitted
n, k = X.shape
s = float(np.sqrt(np.sum(w * resid ** 2) / np.sum(w) * n / max(1, n - k))) # weighted residual RMSE, dof-corrected
ybar = float(np.sum(w * y) / np.sum(w))
r2 = 1.0 - float(np.sum(w * resid ** 2) / np.sum(w * (y - ybar) ** 2)) if np.sum(w * (y - ybar) ** 2) > 0 else float("nan")
return {"beta": beta, "fitted": fitted, "resid": resid, "s": s, "r2": r2, "k": k}
def sampling_se(value: np.ndarray, n: np.ndarray, scale: str) -> np.ndarray:
"""Per-cell sampling SE. scale = the outcome's sampling model: percent (binomial on N students),
rate_per_100 (Poisson on the count), none (no per-cell noise term: z = resid / tau)."""
n = np.where(np.isfinite(n) & (n > 0), n, np.nan)
if scale == "none":
return np.zeros_like(value, dtype=float)
if scale == "percent":
p = np.clip(value / 100.0, 0, 1)
p = np.clip(p, 0.5 / n, 1 - 0.5 / n) # never zero SE at 0% or 100%
return 100.0 * np.sqrt(p * (1 - p) / n)
if scale == "rate_per_100":
lam = np.maximum(value, 0.5 * 100.0 / n) # events per 100 students; Poisson on the count
return np.sqrt(100.0 * lam / n)
raise SystemExit(f"unknown scale {scale}")
def between_variance(s: float, se: np.ndarray, w: np.ndarray) -> float:
"""tau^2: the residual variance net of sampling noise (weighted mean se_i^2 removed once)."""
m = np.isfinite(se)
mean_se2 = float(np.sum(w[m] * se[m] ** 2) / np.sum(w[m])) if m.any() else 0.0
return max(0.0, s ** 2 - mean_se2)
def standardize(resid: np.ndarray, tau2: float, se: np.ndarray) -> np.ndarray:
denom = np.sqrt(tau2 + se ** 2)
return resid / np.where(denom > 0, denom, np.nan)
def orient(z: np.ndarray, direction: str) -> np.ndarray:
return -z if direction == "LOWER_BETTER" else z
def spearman(a: pd.Series, b: pd.Series) -> tuple[float, float, int]:
m = a.notna() & b.notna()
if m.sum() < 10:
return float("nan"), float("nan"), int(m.sum())
r, p = stats.spearmanr(a[m], b[m])
return float(r), float(p), int(m.sum())
def ell_endpoint_sensitivity(e: pd.DataFrame, key: str, cov_cfg: dict, fit: dict,
z: np.ndarray, flag: np.ndarray, se: np.ndarray,
direction: str, z_thr: float) -> dict:
"""Refit the same eligible schools at the ELL upper endpoint.
Neither endpoint is an exact share for a school with hidden LY counts.
This comparison is a sensitivity check, not an interval that contains all
possible fits (schools' hidden counts can vary independently).
"""
lower, upper = "pct_ell_lower_bound", "pct_ell_upper_bound"
if lower not in {c["key"] for c in cov_cfg["numeric"]}:
return {"skipped": "ELL lower-bound covariate is not configured"}
if upper not in e or e[upper].isna().any():
return {"skipped": "ELL upper bounds are missing for eligible schools"}
alternate = e.copy()
alternate[lower] = e[upper]
X, names = design(alternate, cov_cfg)
w = e.n.to_numpy(dtype=float)
upper_fit = fit_wls(X, e[key].to_numpy(dtype=float), w)
upper_tau2 = between_variance(upper_fit["s"], se, w)
upper_z = orient(standardize(upper_fit["resid"], upper_tau2, se), direction)
upper_flags = np.where(upper_z >= z_thr, "positive", np.where(upper_z <= -z_thr, "negative", "none"))
ell_idx = names.index(lower)
return {
"baseline_covariate": lower, "alternate_covariate": upper,
"n_cells": len(e), "n_uncertain": int((e[upper] > e[lower]).sum()),
"baseline_ell_coefficient": float(fit["beta"][ell_idx]),
"upper_ell_coefficient": float(upper_fit["beta"][ell_idx]),
"upper_r2_weighted": float(upper_fit["r2"]), "upper_tau_between": math.sqrt(upper_tau2),
"flag_changes": int((upper_flags != flag).sum()),
"max_abs_z_change": float(np.max(np.abs(upper_z - z))),
"rho_z_baseline_vs_upper": float(stats.spearmanr(z, upper_z)[0]),
"note": "same eligible schools; all lower endpoints versus all upper endpoints, not a bound on every possible hidden-count allocation",
}
def run_outcome(key: str, ocfg: dict, defaults: dict, frame: pd.DataFrame, cov_cfg, comp) -> tuple[pd.DataFrame, dict]:
rule = parse_n_rule(ocfg["n_rule"])
min_n = ocfg.get("min_n", defaults["min_n"])
z_thr = defaults["z_threshold"]
direction = ocfg["direction"]
scale = ocfg["scale"]
se_model = ocfg.get("sampling_error", scale) # outcomes.yaml: the noise model may differ from the value's scale
window_of, excluded, flagged = regime_windows(comp, ocfg["comparability_regime"])
num_keys = [c["key"] for c in cov_cfg["numeric"]]
cat_keys = [c["key"] for c in cov_cfg["categorical"]]
# Fit universe: every school with the outcome that year — closed and inactive schools included (C9).
f = frame[frame[key].notna()].copy()
f["n"] = cell_n(f, key, rule)
f["covariates_complete"] = f[num_keys + cat_keys].notna().all(axis=1)
f["eligible"] = f.covariates_complete & (f.n >= min_n) & ~f.year.isin(excluded)
f["regime"] = f.year.map(window_of)
out_rows = []
per_year = {}
ablation = []
ref_key = ocfg.get("validation_reference")
for year, g in f.groupby("year", sort=True):
e = g[g.eligible].copy()
if len(e) < 30:
per_year[year] = {"skipped": f"only {len(e)} eligible cells"}
continue
X, names = design(e, cov_cfg)
y = e[key].to_numpy(dtype=float)
w = e.n.to_numpy(dtype=float)
fit = fit_wls(X, y, w)
se = sampling_se(y, w, se_model)
tau2 = between_variance(fit["s"], se, w)
z_raw = standardize(fit["resid"], tau2, se)
z = orient(z_raw, direction)
flag = np.where(z >= z_thr, "positive", np.where(z <= -z_thr, "negative", "none"))
e["fitted"], e["resid"], e["se_i"], e["z"], e["flag"] = fit["fitted"], fit["resid"], se, z, flag
out_rows.append(e)
stat = {
"n_cells": int(len(g)), "n_eligible": int(len(e)), "ineligible_low_n": int((((g.n < min_n) | g.n.isna()) & g.covariates_complete & ~g.year.isin(excluded)).sum()),
"ineligible_missing_covariates": int((~g.covariates_complete).sum()), "r2_weighted": round(fit["r2"], 4), "s_weighted_rmse": round(fit["s"], 3),
"tau_between": round(math.sqrt(tau2), 3), "mean_sampling_se": round(float(np.sqrt(np.nansum(w * se ** 2) / np.sum(w[np.isfinite(se)]))) if np.isfinite(se).any() else 0.0, 3),
"sampling_error_model": se_model, "k_parameters": int(fit["k"]),
"n_closed_or_inactive_included": int(((e.activity_code != "A") | e.closed_at.notna()).sum()), "flagged_year": year in flagged,
"beta": {nm: round(float(b), 4) for nm, b in zip(names, fit["beta"])},
"fitted_below_0": int((fit["fitted"] < 0).sum()), "fitted_above_100": int((fit["fitted"] > 100).sum()) if scale == "percent" else 0,
"flagged_positive": int((flag == "positive").sum()), "flagged_negative": int((flag == "negative").sum()),
"z_p10": round(float(np.nanpercentile(z, 10)), 3), "z_p50": round(float(np.nanpercentile(z, 50)), 3), "z_p90": round(float(np.nanpercentile(z, 90)), 3),
"median_n": float(np.median(w)), "median_se_i": round(float(np.nanmedian(se)), 3),
}
# Pre-registered checks
checks = {}
if ref_key and ref_key in e.columns and e[ref_key].notna().sum() >= 30:
r_z, p_z, n_pair = spearman(pd.Series(z, index=e.index), e[ref_key])
r_raw, _, _ = spearman(e[key], e[ref_key])
r_econ_gains, _, _ = spearman(e["pct_econ_disadv"], e[ref_key])
r_fit_gains, _, _ = spearman(pd.Series(fit["fitted"], index=e.index), e[ref_key])
checks["P1"] = {"rho_z_gains": round(r_z, 4), "p": p_z, "n": n_pair, "pass": bool(r_z > 0 and p_z < 0.001)}
checks["P2"] = {"rho_z_gains": round(r_z, 4), "rho_raw_gains": round(r_raw, 4), "rho_econ_gains": round(r_econ_gains, 4), "pass": bool(r_z >= r_raw),
"rho_fitted_gains": round(r_fit_gains, 4), "sd_fitted": round(float(np.std(fit["fitted"])), 3), "sd_resid": round(float(np.std(fit["resid"])), 3)}
l25 = f"{ref_key}_l25"
if l25 in e.columns and e[l25].notna().sum() >= 30:
r_z25, _, n25 = spearman(pd.Series(z, index=e.index), e[l25])
r_raw25, _, _ = spearman(e[key], e[l25])
r_econ25, _, _ = spearman(e["pct_econ_disadv"], e[l25])
r_fit25, _, _ = spearman(pd.Series(fit["fitted"], index=e.index), e[l25])
checks["P2_lowest_quartile_gains"] = {"reference": l25, "n": n25, "rho_z_gains": round(r_z25, 4), "rho_raw_gains": round(r_raw25, 4), "rho_econ_gains": round(r_econ25, 4),
"rho_fitted_gains": round(r_fit25, 4), "z_beats_raw": bool(r_z25 >= r_raw25),
"note": "not pre-registered; the lowest quartile is relative to the school's prior-year assessment groups and uses the same gain rules, including applicable level-holding rules"}
elif ref_key:
checks["P1"] = {"pass": None, "note": f"{ref_key} not published for {year}"}
checks["P2"] = {"pass": None, "note": f"{ref_key} not published for {year}"}
r_ed, _, _ = spearman(pd.Series(z, index=e.index), e["pct_econ_disadv"])
checks["P3"] = {"rho_z_econ": round(r_ed, 4), "pass": bool(abs(r_ed) < 0.10)}
if scale == "percent":
share_out = (stat["fitted_below_0"] + stat["fitted_above_100"]) / len(e)
checks["P4"] = {"share_fitted_outside_0_100": round(share_out, 4), "pass": bool(share_out < 0.02)}
share_flag = (stat["flagged_positive"] + stat["flagged_negative"]) / len(e)
checks["P5"] = {"flagged_share": round(share_flag, 4), "pass": bool(0.02 <= share_flag <= 0.20)}
stat["checks"] = checks
stat["ell_bounds_sensitivity"] = ell_endpoint_sensitivity(e, key, cov_cfg, fit, z, flag, se, direction, z_thr)
per_year[year] = stat
# Ablation: drop one covariate at a time; no covariates (raw standardised value); alternative n rules (weights and se_i).
variants: list[tuple[str, np.ndarray | None, np.ndarray]] = [("ALL (intercept only)", None, w)] + [(k_, None, w) for k_ in num_keys + cat_keys]
for alt in ocfg.get("n_rule_alternatives", []):
variants.append((f"n_rule={alt}", None, cell_n(e, key, parse_n_rule(alt)).to_numpy(dtype=float)))
for label, _, wa in variants:
if label == "ALL (intercept only)":
Xa = np.ones((len(e), 1))
elif label.startswith("n_rule="):
Xa = X
else:
Xa, _ = design(e, cov_cfg, drop={label})
ok = np.isfinite(wa) & (wa > 0)
if not ok.all():
wa = np.where(ok, wa, np.nanmedian(wa[ok]) if ok.any() else 1.0)
fa = fit_wls(Xa, y, wa)
sea = sampling_se(y, wa, se_model)
za = orient(standardize(fa["resid"], between_variance(fa["s"], sea, wa), sea), direction)
fla = np.where(za >= z_thr, "positive", np.where(za <= -z_thr, "negative", "none"))
row = {"year": year, "dropped": label, "r2_weighted": round(fa["r2"], 4),
"rho_z_full_vs_ablated": round(float(stats.spearmanr(z, za)[0]), 4), "flag_changes": int((fla != flag).sum()), "flag_change_share": round(float((fla != flag).mean()), 4)}
if ref_key and ref_key in e.columns:
row["rho_z_gains"] = round(spearman(pd.Series(za, index=e.index), e[ref_key])[0], 4)
ablation.append(row)
cells = pd.concat(out_rows) if out_rows else pd.DataFrame()
if len(cells):
cells = add_persistence(cells, defaults["persistence"])
return cells, {"config": {**ocfg, "min_n": min_n, "z_threshold": z_thr}, "years": per_year, "ablation": ablation}
def add_persistence(cells: pd.DataFrame, pcfg: dict) -> pd.DataFrame:
"""Persistent = the same flag in >= min_flags of the fitted years within the last `window` CALENDAR
years (spring years) inside one regime window, so a missing middle year is a gap, not bridged."""
window, min_flags = int(pcfg["window"]), int(pcfg["min_flags"])
cells = cells.sort_values(["school", "year"]).copy()
cells["persistent_dir"] = "none"
cells["persistent_n"] = 0
for (school, regime), g in cells.groupby(["school", "regime"]):
springs = [spring_year(y) for y in g.year]
flags = list(g.flag)
idx = list(g.index)
for i in range(len(springs)):
win = [flags[j] for j in range(len(springs)) if springs[i] - window < springs[j] <= springs[i]]
for d in ("positive", "negative"):
c = win.count(d)
if c >= min_flags and flags[i] == d:
cells.loc[idx[i], "persistent_dir"] = d
cells.loc[idx[i], "persistent_n"] = c
return cells
def flip_diff(cells: pd.DataFrame, key: str) -> dict:
"""Residual-z flags vs the served peer-group percentile (top/bottom decile)."""
e = cells[cells.eligible & cells[f"peer__{key}"].notna()] if f"peer__{key}" in cells.columns else pd.DataFrame()
if not len(e):
return {"n": 0}
top_pct = e[f"peer__{key}"] >= 0.9
bot_pct = e[f"peer__{key}"] <= 0.1
top_z = e.flag == "positive"
bot_z = e.flag == "negative"
rho = float(stats.spearmanr(e.z, e[f"peer__{key}"])[0])
return {
"n": int(len(e)), "rho_z_vs_peer_percentile": round(rho, 4),
"top": {"percentile_only": int((top_pct & ~top_z).sum()), "z_only": int((top_z & ~top_pct).sum()), "both": int((top_pct & top_z).sum())},
"bottom": {"percentile_only": int((bot_pct & ~bot_z).sum()), "z_only": int((bot_z & ~bot_pct).sum()), "both": int((bot_pct & bot_z).sum())},
}
def main() -> None:
import psycopg
cov_cfg, out_cfg, comp = load_config()
outcome_keys = list(out_cfg["outcomes"].keys())
with psycopg.connect(common.served_db_dsn()) as conn:
conn.read_only = True
frame, sch, raw = fetch_frame(conn, cov_cfg, outcome_keys)
all_cells = []
report = {}
for key, ocfg in out_cfg["outcomes"].items():
if key not in frame.columns:
report[key] = {"skipped": "not served"}
continue
cells, rep = run_outcome(key, ocfg, out_cfg["defaults"], frame, cov_cfg, comp)
if len(cells):
rep["flip_diff_vs_peer_percentile"] = flip_diff(cells, key)
cells = cells.assign(metric_key=key)
all_cells.append(cells[["school", "year", "metric_key", key, "n", "fitted", "resid", "se_i", "z", "flag", "eligible", "persistent_dir", "persistent_n", "regime"]].rename(columns={key: "value", "school": "school_number"}))
report[key] = rep
fitted_years = [y for y, v in rep["years"].items() if "skipped" not in v]
print(f"[residual_z] {key:28s} years {','.join(fitted_years)} cells {len(cells)}")
cells = pd.concat(all_cells) if all_cells else pd.DataFrame()
OUT_CSV.parent.mkdir(parents=True, exist_ok=True)
cells.to_csv(OUT_CSV, index=False, float_format="%.12g")
cfg_hash = hashlib.sha256((open(CONFIG_DIR / "covariates.yaml").read() + open(CONFIG_DIR / "outcomes.yaml").read() + open(CONFIG_DIR / "comparability.yaml").read()).encode()).hexdigest()[:12]
packet = {
"meta": {"generated_at": dt.datetime.now(dt.timezone.utc).isoformat(), "config_sha256_12": cfg_hash, "n_cells": int(len(cells)),
"n_eligible": int(cells.eligible.sum()) if len(cells) else 0, "outcomes": outcome_keys,
"covariates": {"numeric": [c["key"] for c in cov_cfg["numeric"]], "categorical": [c["key"] for c in cov_cfg["categorical"]]},
"method": "WLS per outcome x year on every school with the outcome (closed schools included); z = resid / sqrt(tau^2 + se_i^2) with tau^2 = max(0, s^2 - mean_w se_i^2) and s^2 dof-corrected; flags at |z| >= 2; persistence 2 of the last 3 calendar years within regime. Descriptive, not causal."},
"preregistered": PREREGISTERED,
"outcomes": report,
}
OUT_JSON.write_text(json.dumps(_clean(packet), indent=1, default=str, allow_nan=False))
# Console summary of the pre-registered checks
for key, rep in report.items():
if "years" not in rep:
continue
for y, st in rep["years"].items():
if "checks" in st:
ch = st["checks"]
extra = f" rho(z,gains)={ch['P2']['rho_z_gains']:+.3f} rho(raw,gains)={ch['P2']['rho_raw_gains']:+.3f} rho(econ,gains)={ch['P2']['rho_econ_gains']:+.3f}" if ch.get("P2", {}).get("pass") is not None else ""
print(f" {key:26s} {y}: n={st['n_eligible']:3d} R2={st['r2_weighted']:.3f} flags +{st['flagged_positive']} -{st['flagged_negative']} " + " ".join(f"{k}={'n/a' if v['pass'] is None else 'ok' if v['pass'] else 'FAIL'}" for k, v in ch.items() if "pass" in v) + extra)
print(f"[residual_z] wrote {OUT_CSV} ({len(cells)} cells) and {OUT_JSON}")
if __name__ == "__main__":
main()