""" Purpose: - Run the report spending classsize effects workflow for this repository. Data files required: - `SchoolDataWithAddressesAndCensusSES.csv` - `SchoolDataMathWithAddressesAndCensusSES.csv` - `SchoolDataScienceWithAddressesAndCensusSES.csv` Invocation peculiarities: - Run from repository root: `python3 scripts/report_spending_classsize_effects.py`. - Inputs and outputs are typically controlled by top-level constants in this file. Statistical techniques used: - Correlation and descriptive association analysis. """ import csv from datetime import datetime from pathlib import Path import numpy as np import pandas as pd BASE_DIR = Path(__file__).resolve().parents[1] DATA_DIR = BASE_DIR / "data" / "processed" DOCS_DIR = BASE_DIR / "docs" DATASETS = [ ("English", DATA_DIR / "SchoolDataWithAddressesAndCensusSES.csv"), ("Math", DATA_DIR / "SchoolDataMathWithAddressesAndCensusSES.csv"), ("Science", DATA_DIR / "SchoolDataScienceWithAddressesAndCensusSES.csv"), ] OUT_PATH = DOCS_DIR / "spending_class_size_effects_report.txt" TARGET = "Percent Proficient" WEIGHT = "Scored Performance Denominator" DISTRICT = "District" GROUP = "Student Group" GRADE = "Grade Level" SCHOOL_ID = "School ID" PREDICTORS = [ ("Per Student Spending", "Overall spending per student"), ("Classroom Spending per Student", "Classroom spending per student"), ("Median Class Size For School", "Median class size"), ] CONTROLS = [ ("ACS_Per_Capita_Income", "Per-capita income"), ("ACS_Ed_BA_or_Higher_Rate", "Adult BA+ rate"), ("Attendance_Percent_Regular_Attenders", "Regular attendance rate"), ] def weighted_mean(x, w): wsum = np.sum(w) if wsum <= 0: return np.nan return np.sum(x * w) / wsum def weighted_var(x, w): mu = weighted_mean(x, w) if not np.isfinite(mu): return np.nan wsum = np.sum(w) if wsum <= 0: return np.nan return np.sum(w * (x - mu) ** 2) / wsum def weighted_std(x, w): v = weighted_var(x, w) if not np.isfinite(v) or v <= 0: return np.nan return float(np.sqrt(v)) def weighted_corr(x, y, w): sx = weighted_std(x, w) sy = weighted_std(y, w) if not np.isfinite(sx) or not np.isfinite(sy) or sx <= 0 or sy <= 0: return np.nan mx = weighted_mean(x, w) my = weighted_mean(y, w) cov = np.sum(w * (x - mx) * (y - my)) / np.sum(w) return cov / (sx * sy) def weighted_regression(df, y_col, x_cols, w_col): y = df[y_col].to_numpy(dtype=float) w = df[w_col].to_numpy(dtype=float) X = df[x_cols].to_numpy(dtype=float) X = np.column_stack([np.ones(len(X)), X]) sw = np.sqrt(w) Xw = X * sw[:, None] yw = y * sw beta, _, _, _ = np.linalg.lstsq(Xw, yw, rcond=None) yhat = X @ beta ybar = weighted_mean(y, w) sse = np.sum(w * (y - yhat) ** 2) sst = np.sum(w * (y - ybar) ** 2) r2 = np.nan if sst <= 0 else 1 - sse / sst result = { "intercept": float(beta[0]), "coef": {x_cols[i]: float(beta[i + 1]) for i in range(len(x_cols))}, "r2": float(r2) if np.isfinite(r2) else np.nan, } return result def standardized_betas(df, y_col, x_cols, w_col): w = df[w_col].to_numpy(dtype=float) out = pd.DataFrame(index=df.index) out[w_col] = w for col in [y_col] + x_cols: x = df[col].to_numpy(dtype=float) mu = weighted_mean(x, w) sd = weighted_std(x, w) if not np.isfinite(sd) or sd <= 0: return None out[col] = (x - mu) / sd fit = weighted_regression(out, y_col, x_cols, w_col) return fit["coef"] def weighted_demean_by_group(df, cols, group_col, w_col): out = df.copy() for col in cols: means = {} for g, sub in df.groupby(group_col): x = sub[col].to_numpy(dtype=float) w = sub[w_col].to_numpy(dtype=float) means[g] = weighted_mean(x, w) out[col] = out[col] - out[group_col].map(means) return out def select_total_population_rows(path): with path.open(newline="", encoding="utf-8") as f: rows = list(csv.DictReader(f)) df = pd.DataFrame(rows) if df.empty: return df df = df[df[GROUP].astype(str).str.strip() == "Total Population (All Students)"].copy() if df.empty: return df # Keep school-grade rows and drop redundant "All Grades" rows where school has explicit grades. school_has_non_all = ( df.assign(_grade=df[GRADE].astype(str).str.strip().str.lower()) .query("_grade != 'all grades'") .groupby(SCHOOL_ID) .size() ) has_non_all = school_has_non_all.index.astype(str).tolist() mask_all = df[GRADE].astype(str).str.strip().str.lower() == "all grades" mask_school = df[SCHOOL_ID].astype(str).isin(has_non_all) df = df[~(mask_all & mask_school)].copy() return df def prepare_analysis_frame(df): numeric_cols = [TARGET, WEIGHT, DISTRICT] + [p[0] for p in PREDICTORS] + [c[0] for c in CONTROLS] for col in numeric_cols: if col == DISTRICT: continue df[col] = pd.to_numeric(df[col], errors="coerce") keep_cols = [DISTRICT, TARGET, WEIGHT] + [p[0] for p in PREDICTORS] + [c[0] for c in CONTROLS] frame = df[keep_cols].copy() frame = frame.replace([np.inf, -np.inf], np.nan) frame = frame.dropna(subset=[TARGET, WEIGHT]) frame = frame[frame[WEIGHT] > 0].copy() return frame def fmt(v, d=3): if v is None or (isinstance(v, float) and not np.isfinite(v)): return "n/a" return f"{v:.{d}f}" def run_for_dataset(label, path): base = select_total_population_rows(path) if base.empty: return {"label": label, "error": "No rows after total-population filtering."} frame = prepare_analysis_frame(base) if frame.empty: return {"label": label, "error": "No numeric rows after preparation."} result = { "label": label, "rows": int(len(frame)), "scored_students": float(frame[WEIGHT].sum()), "bivariate": [], "context_assoc": [], "controlled": [], "joint": None, "district_fe": None, } y = frame[TARGET].to_numpy(dtype=float) w = frame[WEIGHT].to_numpy(dtype=float) # Bivariate signals (weighted correlation + simple slope). for col, pretty in PREDICTORS: sub = frame.dropna(subset=[col]).copy() if sub.empty: result["bivariate"].append((pretty, np.nan, np.nan)) continue corr = weighted_corr(sub[col].to_numpy(dtype=float), sub[TARGET].to_numpy(dtype=float), sub[WEIGHT].to_numpy(dtype=float)) fit = weighted_regression(sub, TARGET, [col], WEIGHT) slope = fit["coef"][col] result["bivariate"].append((pretty, corr, slope)) # Predictor association with SES/attendance context vars. for col, pretty in PREDICTORS: row_vals = [pretty] sub = frame.dropna(subset=[col] + [c[0] for c in CONTROLS]).copy() if sub.empty: row_vals.extend([np.nan, np.nan, np.nan]) else: x = sub[col].to_numpy(dtype=float) w_sub = sub[WEIGHT].to_numpy(dtype=float) for c, _ in CONTROLS: row_vals.append(weighted_corr(x, sub[c].to_numpy(dtype=float), w_sub)) result["context_assoc"].append(tuple(row_vals)) # Controlled models: each spending/class-size predictor with SES+attendance controls. control_cols = [c[0] for c in CONTROLS] for col, pretty in PREDICTORS: model_cols = [col] + control_cols sub = frame.dropna(subset=model_cols + [TARGET, WEIGHT]).copy() if sub.empty: result["controlled"].append((pretty, np.nan, np.nan, np.nan)) continue fit = weighted_regression(sub, TARGET, model_cols, WEIGHT) betas = standardized_betas(sub, TARGET, model_cols, WEIGHT) beta_std = np.nan if betas is None else betas.get(col, np.nan) result["controlled"].append((pretty, fit["coef"][col], beta_std, fit["r2"])) # Joint model with all three predictors + controls. joint_cols = [p[0] for p in PREDICTORS] + control_cols sub_joint = frame.dropna(subset=joint_cols + [TARGET, WEIGHT]).copy() if not sub_joint.empty: joint_fit = weighted_regression(sub_joint, TARGET, joint_cols, WEIGHT) joint_betas = standardized_betas(sub_joint, TARGET, joint_cols, WEIGHT) or {} result["joint"] = { "n": int(len(sub_joint)), "r2": joint_fit["r2"], "coef": joint_fit["coef"], "beta": joint_betas, } # Within-district fixed-effects style check via weighted demeaning. if not sub_joint.empty: dm_cols = [TARGET] + joint_cols dm = weighted_demean_by_group(sub_joint[[DISTRICT, WEIGHT] + dm_cols].copy(), dm_cols, DISTRICT, WEIGHT) dm_fit = weighted_regression(dm, TARGET, joint_cols, WEIGHT) dm_betas = standardized_betas(dm, TARGET, joint_cols, WEIGHT) or {} result["district_fe"] = { "n": int(len(dm)), "r2": dm_fit["r2"], "coef": dm_fit["coef"], "beta": dm_betas, } # Decile slope hints for visibility of weak but monotonic patterns. deciles = {} for col, pretty in PREDICTORS: sub = frame.dropna(subset=[col, TARGET, WEIGHT]).copy() if len(sub) < 20: continue sub["_q"] = pd.qcut(sub[col], q=10, labels=False, duplicates="drop") grouped = sub.groupby("_q") points = [] for q, g in grouped: xv = weighted_mean(g[col].to_numpy(float), g[WEIGHT].to_numpy(float)) yv = weighted_mean(g[TARGET].to_numpy(float), g[WEIGHT].to_numpy(float)) points.append((float(xv), float(yv))) if len(points) >= 2: deciles[pretty] = (points[0], points[-1], points[-1][1] - points[0][1]) result["deciles"] = deciles return result def write_report(results): lines = [] lines.append("Spending/Class-Size vs Performance: Numerical Signal Check") lines.append(f"Generated: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}") lines.append("") lines.append("Design") lines.append("- Outcome: Percent Proficient") lines.append("- Rows: Total Population school rows; \"All Grades\" rows dropped when school has explicit grade rows") lines.append("- Weighting: students with reported Level 1-4 results") lines.append("- Predictors tested: Overall spending per student, Classroom spending per student, Median class size") lines.append("- Controls: Per-capita income, Adult BA+ rate, Regular attendance rate") lines.append("- Extra checks: joint model and within-district de-meaned (fixed-effects style) model") lines.append("") for r in results: lines.append(f"{r['label']} dataset") if "error" in r: lines.append(f"- ERROR: {r['error']}") lines.append("") continue lines.append(f"- Rows analyzed: {r['rows']:,}") lines.append(f"- Students with reported Level 1-4 results: {r['scored_students']:,.0f}") lines.append("- Bivariate weighted signals") for pretty, corr, slope in r["bivariate"]: lines.append(f" - {pretty}: corr={fmt(corr)}, slope={fmt(slope,4)}") lines.append("- Controlled models (predictor + SES/attendance controls)") for pretty, coef, beta_std, r2 in r["controlled"]: lines.append( f" - {pretty}: coef={fmt(coef,4)}, std_beta={fmt(beta_std,4)}, model_R2={fmt(r2,4)}" ) lines.append("- Predictor association with context variables (weighted corr)") lines.append(" - Columns: with Per-capita income, Adult BA+ rate, Regular attendance") for pretty, c1, c2, c3 in r["context_assoc"]: lines.append(f" - {pretty}: {fmt(c1,3)}, {fmt(c2,3)}, {fmt(c3,3)}") if r["joint"]: j = r["joint"] lines.append(f"- Joint model (all three predictors + controls), n={j['n']:,}, R2={fmt(j['r2'],4)}") for col, pretty in PREDICTORS: lines.append( f" - {pretty}: coef={fmt(j['coef'].get(col),4)}, std_beta={fmt(j['beta'].get(col),4)}" ) if r["district_fe"]: d = r["district_fe"] lines.append(f"- Within-district check (demeaned), n={d['n']:,}, R2={fmt(d['r2'],4)}") for col, pretty in PREDICTORS: lines.append( f" - {pretty}: coef={fmt(d['coef'].get(col),4)}, std_beta={fmt(d['beta'].get(col),4)}" ) if r.get("deciles"): lines.append("- Decile edge-to-edge change in Percent Proficient (Q10 minus Q1)") for pretty, (p0, p9, dy) in r["deciles"].items(): lines.append( f" - {pretty}: y@Q1={fmt(p0[1],2)} -> y@Q10={fmt(p9[1],2)} (delta={fmt(dy,2)} pts)" ) lines.append("") # High-level readout: lines.append("Interpretation guide") lines.append("- If bivariate effects are visible but controlled std_beta is near zero, gross correlation is likely confounded by SES/attendance.") lines.append("- If within-district std_beta is near zero, differences are mostly between districts rather than within-district resource variation.") lines.append("- Joint-model std_beta allows apples-to-apples comparison of relative effect sizes.") lines.append("") OUT_PATH.write_text("\n".join(lines), encoding="utf-8") def main(): results = [run_for_dataset(label, path) for label, path in DATASETS] write_report(results) print(f"Wrote {OUT_PATH}") if __name__ == "__main__": main()