"""
Build an exploratory read of Andrew Ho's proficiency-threshold critique
against Oregon 2024-25 school achievement data.

This script intentionally works only with the achievement-level counts that
are available in the public ODE extracts. It does not infer scale-score means
or percentiles.
"""

from __future__ import annotations

import os
from pathlib import Path
import textwrap

os.environ.setdefault("MPLCONFIGDIR", "/tmp/orschoolanalysis_matplotlib")

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


ROOT = Path(__file__).resolve().parents[3]
OUT_DIR = Path(__file__).resolve().parent
DATE = "2026-04-26"

SUBJECT_INPUTS = {
    "English Language Arts": ROOT / "data/processed/SchoolDataWithAddressesAndCensusSES.csv",
    "Mathematics": ROOT / "data/processed/SchoolDataMathWithAddressesAndCensusSES.csv",
    "Science": ROOT / "data/processed/SchoolDataScienceWithAddressesAndCensusSES.csv",
}

LEVEL_COUNT_COLS = [f"Number Level {level}" for level in range(1, 5)]
CONTEXT_FIELDS = [
    "ACS_Ed_BA_or_Higher_Rate",
    "Students Experiencing Poverty",
    "Attendance_Percent_Regular_Attenders",
]
PREDICTOR_LABELS = {
    "ACS_Ed_BA_or_Higher_Rate": "BA+ adults",
    "Students Experiencing Poverty": "Students in poverty",
    "Attendance_Percent_Regular_Attenders": "Regular attendance",
}
OUTCOME_LABELS = {
    "pct_proficient": "Percent proficient",
    "level_index": "4-level index",
    "pct_l4": "Percent Level 4",
    "pct_l1": "Percent Level 1",
    "near_cut_pct": "Level 2+3 share",
}
STRESS_SCHEMES = {
    "Percent proficient": (0, 0, 1, 1),
    "Equal 1-4": (1, 2, 3, 4),
    "Level 4 premium": (1, 2, 3, 5),
    "Strong Level 4 premium": (1, 2, 3, 6),
    "Level 1 penalty": (0, 2, 3, 4),
    "Compressed middle": (1, 1.7, 2.3, 4),
    "Near-pass friendly": (1, 2.5, 3, 4),
    "Level 4 only": (0, 0, 0, 1),
    "Avoid Level 1": (0, 1, 1, 1),
}
CORE_STRESS_SCHEMES = [
    "Percent proficient",
    "Equal 1-4",
    "Level 4 premium",
    "Strong Level 4 premium",
    "Level 1 penalty",
    "Compressed middle",
    "Near-pass friendly",
]

COLORS = {
    "Level 1": "#b23a48",
    "Level 2": "#e0a526",
    "Level 3": "#2a9d8f",
    "Level 4": "#31588a",
}
SUBJECT_COLORS = {
    "English Language Arts": "#2a9d8f",
    "Mathematics": "#31588a",
    "Science": "#b23a48",
}


def is_yes(series: pd.Series) -> pd.Series:
    return series.astype(str).str.strip().str.lower().isin({"yes", "true", "1"})


def clean_number(series: pd.Series) -> pd.Series:
    return pd.to_numeric(series, errors="coerce")


def weighted_average(values: pd.Series, weights: pd.Series) -> float:
    mask = values.notna() & weights.notna() & (weights > 0)
    if not mask.any():
        return np.nan
    return float(np.average(values[mask], weights=weights[mask]))


def weighted_corr(x: pd.Series, y: pd.Series, w: pd.Series) -> float:
    x_arr = pd.to_numeric(x, errors="coerce").to_numpy(dtype=float)
    y_arr = pd.to_numeric(y, errors="coerce").to_numpy(dtype=float)
    w_arr = pd.to_numeric(w, errors="coerce").to_numpy(dtype=float)
    mask = np.isfinite(x_arr) & np.isfinite(y_arr) & np.isfinite(w_arr) & (w_arr > 0)
    if mask.sum() < 3:
        return np.nan
    x_arr = x_arr[mask]
    y_arr = y_arr[mask]
    w_arr = w_arr[mask]
    x_mean = np.average(x_arr, weights=w_arr)
    y_mean = np.average(y_arr, weights=w_arr)
    cov = np.average((x_arr - x_mean) * (y_arr - y_mean), weights=w_arr)
    x_var = np.average((x_arr - x_mean) ** 2, weights=w_arr)
    y_var = np.average((y_arr - y_mean) ** 2, weights=w_arr)
    if x_var <= 0 or y_var <= 0:
        return np.nan
    return float(cov / np.sqrt(x_var * y_var))


def build_school_subject_dataset() -> pd.DataFrame:
    frames: list[pd.DataFrame] = []

    for subject, path in SUBJECT_INPUTS.items():
        df = pd.read_csv(path, low_memory=False)
        df = df[df["Student Group"].eq("Total Population (All Students)")].copy()
        df = df[df["Grade Level"].ne("All Grades")]
        ordinary_mask = ~is_yes(df["Is Charter School"]) & ~is_yes(df["Is Virtual School"])
        if "Is Special Enrollment School" in df.columns:
            ordinary_mask &= ~is_yes(df["Is Special Enrollment School"])
        df = df[ordinary_mask]

        numeric_cols = LEVEL_COUNT_COLS + [
            "Number of Participants",
            "Scored Performance Denominator",
            "Participation Rate",
        ] + CONTEXT_FIELDS
        for col in numeric_cols:
            df[col] = clean_number(df[col])
        df = df[df["Scored Performance Denominator"] > 0]

        aggs: dict[str, str] = {
            col: "sum"
            for col in LEVEL_COUNT_COLS
            + ["Number of Participants", "Scored Performance Denominator"]
        }
        for field in CONTEXT_FIELDS:
            df[f"{field}__wx"] = df[field] * df["Scored Performance Denominator"]
            df[f"{field}__w"] = np.where(
                df[field].notna(), df["Scored Performance Denominator"], 0
            )
            aggs[f"{field}__wx"] = "sum"
            aggs[f"{field}__w"] = "sum"

        grouped = (
            df.groupby(["School ID", "School", "District", "Subject"], dropna=False)
            .agg(aggs)
            .reset_index()
        )
        grouped["Subject"] = subject
        grouped["level_count_sum"] = grouped[LEVEL_COUNT_COLS].sum(axis=1)
        grouped["level_count_gap"] = grouped["Number of Participants"] - grouped["level_count_sum"]

        n = grouped["Scored Performance Denominator"]
        for level in range(1, 5):
            grouped[f"pct_l{level}"] = 100 * grouped[f"Number Level {level}"] / n
        grouped["pct_not_proficient"] = grouped["pct_l1"] + grouped["pct_l2"]
        grouped["pct_proficient"] = grouped["pct_l3"] + grouped["pct_l4"]
        grouped["level_index"] = (
            grouped["Number Level 1"]
            + 2 * grouped["Number Level 2"]
            + 3 * grouped["Number Level 3"]
            + 4 * grouped["Number Level 4"]
        ) / n
        grouped["near_cut_pct"] = grouped["pct_l2"] + grouped["pct_l3"]
        grouped["level4_of_proficient_pct"] = np.where(
            grouped["pct_proficient"] > 0,
            100 * grouped["pct_l4"] / grouped["pct_proficient"],
            np.nan,
        )

        for field in CONTEXT_FIELDS:
            grouped[field] = grouped[f"{field}__wx"] / grouped[f"{field}__w"]

        frames.append(grouped)

    school_subject = pd.concat(frames, ignore_index=True)
    school_subject = school_subject[
        school_subject["Scored Performance Denominator"] >= 30
    ].copy()
    school_subject["pp_nearest_5"] = (school_subject["pct_proficient"] / 5).round() * 5
    return school_subject


def subject_summary(school_subject: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for subject, sub in school_subject.groupby("Subject", sort=False):
        weights = sub["Scored Performance Denominator"]
        rows.append(
            {
                "subject": subject,
                "school_subject_rows": len(sub),
                "scored_students": int(weights.sum()),
                "percent_proficient": weighted_average(sub["pct_proficient"], weights),
                "four_level_index": weighted_average(sub["level_index"], weights),
                "percent_level_1": weighted_average(sub["pct_l1"], weights),
                "percent_level_2": weighted_average(sub["pct_l2"], weights),
                "percent_level_3": weighted_average(sub["pct_l3"], weights),
                "percent_level_4": weighted_average(sub["pct_l4"], weights),
                "level_2_3_share": weighted_average(sub["near_cut_pct"], weights),
                "median_school_level_2_3_share": float(sub["near_cut_pct"].median()),
                "share_school_subjects_level_2_3_ge_50": float((sub["near_cut_pct"] >= 50).mean() * 100),
                "weighted_corr_pp_index": weighted_corr(sub["pct_proficient"], sub["level_index"], weights),
            }
        )
    return pd.DataFrame(rows)


def pair_examples(school_subject: pd.DataFrame) -> pd.DataFrame:
    examples = []
    candidates = school_subject[
        school_subject["Scored Performance Denominator"] >= 80
    ].copy()

    for subject, sub in candidates.groupby("Subject", sort=False):
        sub = sub.sort_values("pct_proficient").reset_index(drop=True)
        best = None
        for idx, row in sub.iterrows():
            close = sub[
                (sub["pct_proficient"] >= row["pct_proficient"])
                & (sub["pct_proficient"] <= row["pct_proficient"] + 1.0)
            ]
            if len(close) < 2:
                continue
            low = close.loc[close["level_index"].idxmin()]
            high = close.loc[close["level_index"].idxmax()]
            diff = high["level_index"] - low["level_index"]
            if best is None or diff > best[0]:
                best = (diff, low, high)

        if best is None:
            continue

        diff, low, high = best
        for label, rec in [("lower_index", low), ("higher_index", high)]:
            examples.append(
                {
                    "subject": subject,
                    "comparison_role": label,
                    "school": rec["School"],
                    "district": rec["District"],
                    "scored_students": rec["Scored Performance Denominator"],
                    "percent_proficient": rec["pct_proficient"],
                    "four_level_index": rec["level_index"],
                    "percent_level_1": rec["pct_l1"],
                    "percent_level_2": rec["pct_l2"],
                    "percent_level_3": rec["pct_l3"],
                    "percent_level_4": rec["pct_l4"],
                    "index_difference_with_pair": diff,
                }
            )

    return pd.DataFrame(examples)


def correlation_sensitivity(school_subject: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for subject, sub in school_subject.groupby("Subject", sort=False):
        weights = sub["Scored Performance Denominator"]
        for outcome, outcome_label in OUTCOME_LABELS.items():
            for predictor, predictor_label in PREDICTOR_LABELS.items():
                rows.append(
                    {
                        "subject": subject,
                        "outcome": outcome,
                        "outcome_label": outcome_label,
                        "predictor": predictor,
                        "predictor_label": predictor_label,
                        "weighted_r": weighted_corr(sub[predictor], sub[outcome], weights),
                        "n_school_subject_rows": int((sub[predictor].notna() & sub[outcome].notna()).sum()),
                        "scored_student_weight": int(
                            weights[sub[predictor].notna() & sub[outcome].notna()].sum()
                        ),
                    }
                )
    return pd.DataFrame(rows)


def stress_test_level_spacing(school_subject: pd.DataFrame) -> pd.DataFrame:
    rows = []
    equal_score = (
        school_subject["Number Level 1"]
        + 2 * school_subject["Number Level 2"]
        + 3 * school_subject["Number Level 3"]
        + 4 * school_subject["Number Level 4"]
    ) / school_subject["Scored Performance Denominator"]

    for subject, sub in school_subject.groupby("Subject", sort=False):
        subject_equal = equal_score.loc[sub.index]
        for scheme, level_scores in STRESS_SCHEMES.items():
            score = (
                level_scores[0] * sub["Number Level 1"]
                + level_scores[1] * sub["Number Level 2"]
                + level_scores[2] * sub["Number Level 3"]
                + level_scores[3] * sub["Number Level 4"]
            ) / sub["Scored Performance Denominator"]
            row = {
                "subject": subject,
                "scheme": scheme,
                "level_1_score": level_scores[0],
                "level_2_score": level_scores[1],
                "level_3_score": level_scores[2],
                "level_4_score": level_scores[3],
                "is_core_stress_scheme": scheme in CORE_STRESS_SCHEMES,
                "weighted_mean": weighted_average(
                    score, sub["Scored Performance Denominator"]
                ),
                "rank_corr_vs_equal_index": weighted_corr(
                    subject_equal.rank(method="average"),
                    score.rank(method="average"),
                    pd.Series(1, index=sub.index),
                ),
            }
            for predictor, predictor_label in PREDICTOR_LABELS.items():
                row[f"weighted_r__{predictor}"] = weighted_corr(
                    sub[predictor], score, sub["Scored Performance Denominator"]
                )
                row[f"weighted_r_label__{predictor}"] = predictor_label
            rows.append(row)
    return pd.DataFrame(rows)


def plot_statewide_levels(summary: pd.DataFrame) -> Path:
    fig, ax = plt.subplots(figsize=(10, 5.8))
    subjects = summary["subject"].tolist()
    y = np.arange(len(subjects))
    left = np.zeros(len(subjects))

    level_fields = [
        ("Level 1", "percent_level_1"),
        ("Level 2", "percent_level_2"),
        ("Level 3", "percent_level_3"),
        ("Level 4", "percent_level_4"),
    ]
    for label, field in level_fields:
        values = summary[field].to_numpy()
        ax.barh(y, values, left=left, label=label, color=COLORS[label], height=0.62)
        for i, value in enumerate(values):
            if value >= 6:
                ax.text(left[i] + value / 2, y[i], f"{value:.0f}%", ha="center", va="center", color="white", fontsize=9)
        left += values

    for i, row in summary.iterrows():
        ax.text(101, i, f"Proficient: {row['percent_proficient']:.1f}%", va="center", fontsize=10)

    ax.set_yticks(y, subjects)
    ax.invert_yaxis()
    ax.set_xlim(0, 118)
    ax.set_xlabel("Student-weighted share of tested students")
    ax.set_title("Oregon 2024-25 achievement levels: proficiency is Level 3 + Level 4")
    ax.legend(ncol=4, loc="lower center", bbox_to_anchor=(0.5, -0.22), frameon=False)
    ax.spines[["top", "right", "left"]].set_visible(False)
    ax.grid(axis="x", color="#d9d9d9", linewidth=0.8, alpha=0.7)
    fig.tight_layout()
    out = OUT_DIR / f"ho_proficiency_thresholds_statewide_levels_{DATE}.png"
    fig.savefig(out, dpi=180)
    plt.close(fig)
    return out


def plot_pp_vs_index(school_subject: pd.DataFrame) -> Path:
    fig, ax = plt.subplots(figsize=(9.5, 6.2))
    for subject, sub in school_subject.groupby("Subject", sort=False):
        sizes = np.clip(np.sqrt(sub["Scored Performance Denominator"]) * 5, 18, 95)
        ax.scatter(
            sub["pct_proficient"],
            sub["level_index"],
            s=sizes,
            alpha=0.34,
            color=SUBJECT_COLORS[subject],
            label=subject,
            linewidths=0,
        )

    ax.set_xlim(-2, 102)
    ax.set_ylim(1.0, 4.0)
    ax.set_xlabel("Percent proficient")
    ax.set_ylabel("Four-level achievement index (Level 1=1, Level 4=4)")
    ax.set_title("Same threshold result can mask different achievement distributions")
    ax.text(
        2,
        3.78,
        "The vertical spread is the information lost when Level 1-4 counts collapse to one cut score.",
        fontsize=9.5,
        color="#333333",
        bbox={"facecolor": "white", "edgecolor": "#cccccc", "boxstyle": "round,pad=0.35"},
    )
    ax.legend(frameon=False, loc="lower right")
    ax.grid(color="#d9d9d9", linewidth=0.8, alpha=0.7)
    ax.spines[["top", "right"]].set_visible(False)
    fig.tight_layout()
    out = OUT_DIR / f"ho_proficiency_thresholds_pp_vs_index_{DATE}.png"
    fig.savefig(out, dpi=180)
    plt.close(fig)
    return out


def plot_same_pp_pairs(pairs: pd.DataFrame) -> Path:
    fig, ax = plt.subplots(figsize=(10.8, 6.8))
    pair_order = []
    labels = []
    for subject in SUBJECT_INPUTS:
        sub = pairs[pairs["subject"].eq(subject)]
        for _, row in sub.iterrows():
            pair_order.append(row)
            school_label = row["school"]
            if len(school_label) > 28:
                school_label = school_label[:25] + "..."
            labels.append(
                f"{subject.replace('English Language Arts', 'ELA')} - {school_label}\n"
                f"PP {row['percent_proficient']:.1f}%, index {row['four_level_index']:.2f}"
            )

    y = np.arange(len(pair_order))
    left = np.zeros(len(pair_order))
    for level, field in [
        ("Level 1", "percent_level_1"),
        ("Level 2", "percent_level_2"),
        ("Level 3", "percent_level_3"),
        ("Level 4", "percent_level_4"),
    ]:
        values = np.array([row[field] for row in pair_order])
        ax.barh(y, values, left=left, color=COLORS[level], label=level, height=0.68)
        left += values

    ax.set_yticks(y, labels, fontsize=8.5)
    ax.invert_yaxis()
    ax.set_xlim(0, 100)
    ax.set_xlabel("Share of students by achievement level")
    ax.set_title("Matched examples: nearly identical percent proficient, different level mix")
    ax.legend(ncol=4, loc="lower center", bbox_to_anchor=(0.5, -0.2), frameon=False)
    ax.grid(axis="x", color="#d9d9d9", linewidth=0.8, alpha=0.7)
    ax.spines[["top", "right", "left"]].set_visible(False)
    fig.tight_layout()
    out = OUT_DIR / f"ho_proficiency_thresholds_same_pp_pairs_{DATE}.png"
    fig.savefig(out, dpi=180)
    plt.close(fig)
    return out


def plot_correlation_sensitivity(corrs: pd.DataFrame) -> Path:
    subjects = list(SUBJECT_INPUTS.keys())
    outcomes = list(OUTCOME_LABELS.values())
    predictors = list(PREDICTOR_LABELS.values())

    fig, axes = plt.subplots(1, 3, figsize=(14.2, 5.6), sharey=True)
    for ax, subject in zip(axes, subjects):
        sub = corrs[corrs["subject"].eq(subject)]
        matrix = (
            sub.pivot(index="outcome_label", columns="predictor_label", values="weighted_r")
            .reindex(index=outcomes, columns=predictors)
        )
        image = ax.imshow(matrix.to_numpy(), vmin=-0.8, vmax=0.8, cmap="RdBu")
        ax.set_xticks(np.arange(len(predictors)), predictors, rotation=35, ha="right", fontsize=8.5)
        ax.set_yticks(np.arange(len(outcomes)), outcomes, fontsize=9)
        ax.set_title(subject.replace("English Language Arts", "ELA"), fontsize=11)
        for i in range(matrix.shape[0]):
            for j in range(matrix.shape[1]):
                value = matrix.iat[i, j]
                color = "white" if abs(value) > 0.48 else "#222222"
                ax.text(j, i, f"{value:.2f}", ha="center", va="center", fontsize=8.5, color=color)
        ax.tick_params(length=0)

    fig.subplots_adjust(left=0.18, right=0.90, top=0.82, bottom=0.26, wspace=0.34)
    cax = fig.add_axes([0.925, 0.28, 0.015, 0.50])
    cbar = fig.colorbar(image, cax=cax)
    cbar.set_label("Weighted Pearson r")
    fig.suptitle("Context correlations change when the outcome moves away from the proficiency cut", y=0.98)
    out = OUT_DIR / f"ho_proficiency_thresholds_context_correlation_sensitivity_{DATE}.png"
    fig.savefig(out, dpi=180)
    plt.close(fig)
    return out


def write_memo(
    summary: pd.DataFrame,
    pairs: pd.DataFrame,
    corrs: pd.DataFrame,
    stress: pd.DataFrame,
    figure_paths: list[Path],
) -> Path:
    def subject_line(subject: str) -> str:
        row = summary[summary["subject"].eq(subject)].iloc[0]
        return (
            f"- {subject}: {row['percent_proficient']:.1f}% proficient, "
            f"{row['percent_level_1']:.1f}% Level 1, {row['percent_level_2']:.1f}% Level 2, "
            f"{row['percent_level_3']:.1f}% Level 3, {row['percent_level_4']:.1f}% Level 4; "
            f"Level 2+3 adjacent-band share is {row['level_2_3_share']:.1f}%."
        )

    pair_lines = []
    for subject, sub in pairs.groupby("subject", sort=False):
        low = sub[sub["comparison_role"].eq("lower_index")].iloc[0]
        high = sub[sub["comparison_role"].eq("higher_index")].iloc[0]
        pair_lines.append(
            f"- {subject}: {low['school']} and {high['school']} are both around "
            f"{low['percent_proficient']:.1f}-{high['percent_proficient']:.1f}% proficient, "
            f"but their four-level indexes are {low['four_level_index']:.2f} and "
            f"{high['four_level_index']:.2f}."
        )

    corr_lines = []
    for subject in SUBJECT_INPUTS:
        sub = corrs[(corrs["subject"].eq(subject)) & (corrs["outcome"].isin(["pct_proficient", "level_index"]))]
        pp = sub[sub["outcome"].eq("pct_proficient")]
        idx = sub[sub["outcome"].eq("level_index")]
        corr_lines.append(
            f"- {subject}: Percent proficient correlations are "
            f"BA+ {pp[pp['predictor_label'].eq('BA+ adults')]['weighted_r'].iloc[0]:+.2f}, "
            f"poverty {pp[pp['predictor_label'].eq('Students in poverty')]['weighted_r'].iloc[0]:+.2f}, "
            f"attendance {pp[pp['predictor_label'].eq('Regular attendance')]['weighted_r'].iloc[0]:+.2f}; "
            f"the four-level index is similar at "
            f"{idx[idx['predictor_label'].eq('BA+ adults')]['weighted_r'].iloc[0]:+.2f}, "
            f"{idx[idx['predictor_label'].eq('Students in poverty')]['weighted_r'].iloc[0]:+.2f}, "
            f"{idx[idx['predictor_label'].eq('Regular attendance')]['weighted_r'].iloc[0]:+.2f}."
        )

    stress_lines = []
    core_stress = stress[stress["is_core_stress_scheme"]].copy()
    for subject, sub in core_stress.groupby("subject", sort=False):
        ba = sub["weighted_r__ACS_Ed_BA_or_Higher_Rate"]
        poverty = sub["weighted_r__Students Experiencing Poverty"]
        attendance = sub["weighted_r__Attendance_Percent_Regular_Attenders"]
        stress_lines.append(
            f"- {subject}: across the core scoring alternatives, BA+ ranges "
            f"{ba.min():+.2f} to {ba.max():+.2f}, poverty {poverty.min():+.2f} to "
            f"{poverty.max():+.2f}, and attendance {attendance.min():+.2f} to "
            f"{attendance.max():+.2f}."
        )

    extreme_lines = []
    extreme_stress = stress[~stress["is_core_stress_scheme"]].copy()
    for subject, sub in extreme_stress.groupby("subject", sort=False):
        l4 = sub[sub["scheme"].eq("Level 4 only")].iloc[0]
        avoid_l1 = sub[sub["scheme"].eq("Avoid Level 1")].iloc[0]
        extreme_lines.append(
            f"- {subject}: the edge-case `Level 4 only` view has BA+ "
            f"{l4['weighted_r__ACS_Ed_BA_or_Higher_Rate']:+.2f}, poverty "
            f"{l4['weighted_r__Students Experiencing Poverty']:+.2f}, attendance "
            f"{l4['weighted_r__Attendance_Percent_Regular_Attenders']:+.2f}; "
            f"the edge-case `Avoid Level 1` view has BA+ "
            f"{avoid_l1['weighted_r__ACS_Ed_BA_or_Higher_Rate']:+.2f}, poverty "
            f"{avoid_l1['weighted_r__Students Experiencing Poverty']:+.2f}, attendance "
            f"{avoid_l1['weighted_r__Attendance_Percent_Regular_Attenders']:+.2f}."
        )

    figure_lines = "\n".join(f"- `{path.name}`" for path in figure_paths)
    memo = f"""# Ho Proficiency Thresholds: Oregon 2024-25 Exploratory Read

Generated: {DATE}

## Question

Andrew Ho's 2008 paper argues that Percentage of Proficient Students statistics can be technically correct yet misleading as summaries of score distributions. This note asks what that concern means for our Oregon 2024-25 school achievement analyses, which often use `Percent Proficient` as the main outcome.

## Scope

- Data: Oregon 2024-25 school achievement rows in `data/processed/*WithAddressesAndCensusSES.csv`.
- Rows: `Total Population (All Students)`, grade-specific rows only. `All Grades` science rows are excluded to avoid double-counting.
- Filters: charter, virtual, and curated special-enrollment schools excluded.
- Aggregation: school-subject totals across tested grades, using achievement-level counts.
- Inclusion threshold: at least 30 students with scored Level 1-4 results after aggregation.
- Unit: school-subject rows. Scored-student counts are subject-counts, not unique students across subjects.

## Ho's Concern, Applied Here

Ho's central statistical point is that proficiency percentages depend on the selected cut score and on how many students are near that cut (Ho, 2008). A trend, gap, or school comparison can look large, small, or even directionally different depending on which threshold is used. He recommends distribution-wide views such as averages, effect sizes, percentiles, and multi-point graphical summaries.

Oregon's public school-level data do not give us scale-score means, student percentiles, or full score distributions. They do give us counts in Levels 1-4. That means we cannot fully solve Ho's problem, but we can reduce it by checking whether `Percent Proficient` tells the same story as Level 1, Level 4, and a simple four-level index.

## Statewide Level Mix

{subject_line("English Language Arts")}
{subject_line("Mathematics")}
{subject_line("Science")}

Interpretation: the proficiency cut is not summarizing the same distributional shape in each subject. Science, for example, has a large Level 2+3 adjacent-band share but very little Level 4, so a single `Level 3+4` rate can blur a different story than it does in ELA.

## Same Percent Proficient, Different Distribution

{chr(10).join(pair_lines)}

These are not claims about the named schools themselves; they are diagnostic examples. The point is that the official threshold metric can be almost identical while the underlying Level 1-4 mix is materially different.

## Effect on Our Core Context Analyses

{chr(10).join(corr_lines)}

This is reassuring but not a free pass. The big context story does not disappear when we replace `Percent Proficient` with the crude four-level index: poverty, attendance, and adult BA+ context still point in similar directions. But the sensitivity plot shows that Level 4 and Level 1 are not interchangeable with proficiency. Top-end performance, severe low-end performance, and near-cut concentration each emphasize different parts of the distribution.

## Stress Test: Does Equal Spacing Drive the Finding?

Question answered: if the four-level index is only a rough proxy because Level 1, Level 2, Level 3, and Level 4 may not be evenly spaced, do the SES associations disappear when the level scores are changed?

I tested several alternative scoring rules: the original proficient/not-proficient threshold, equal `1-2-3-4` spacing, extra premium for Level 4, stronger premium for Level 4, stronger penalty for Level 1, a compressed middle, and a near-pass-friendly rule that gives Level 2 more credit. These are not official scales; they are robustness checks.

{chr(10).join(stress_lines)}

Interpretation: the exact index value is not sacred, but the main SES story is not driven by the equal-spacing assumption. Adult BA+, poverty, and attendance remain strongly associated with the achievement distribution under the core alternate scoring rules.

I also included two edge-case lenses, not as preferred indexes but as reminders that "success" can mean different parts of the distribution:

{chr(10).join(extreme_lines)}

Interpretation: Level 4-only and avoiding Level 1 are different questions. Level 4-only puts more weight on excellence; avoiding Level 1 puts more weight on severe low-end performance. These variants change some magnitudes, especially in science, which reinforces the caution about precise school rankings and tail-specific claims.

## Practical Read

1. `Percent Proficient` is acceptable as a descriptive snapshot, but it should not carry the whole evidentiary load.
2. For school comparisons, pair it with stacked Level 1-4 distributions or at least Level 1 and Level 4 rates.
3. For "overperforming" or "underperforming" claims, rerun the model with Level 4, Level 1, the four-level index, and at least one alternate level-spacing rule as sensitivity checks.
4. For trend claims, be especially cautious. Without scale-score means or percentiles, Oregon school-level data cannot tell whether students are moving across the whole distribution or mostly around the Level 2/3 threshold.
5. The best long-run fix would be ODE release of school-level scale-score summaries or percentile points, not just achievement levels.

## Figures

{figure_lines}

## Machine-Readable Outputs

- `ho_proficiency_thresholds_school_subject_{DATE}.csv`
- `ho_proficiency_thresholds_summary_{DATE}.csv`
- `ho_proficiency_thresholds_same_pp_examples_{DATE}.csv`
- `ho_proficiency_thresholds_correlation_sensitivity_{DATE}.csv`
- `ho_proficiency_thresholds_level_spacing_stress_test_{DATE}.csv`

## Caveats

The four-level index treats Level 1, Level 2, Level 3, and Level 4 as evenly spaced. That is a useful sensitivity check, not an official scale score. The Level 2+3 share is an adjacent-band proxy for threshold sensitivity, not a count of students literally near the cut score. The underlying public data remain aggregate, not student-level.

## Reference

Ho, A. D. (2008). The problem with "proficiency": Limitations of statistics and policy under No Child Left Behind. *Educational Researcher, 37*(6), 351-360. https://doi.org/10.3102/0013189X08323842
"""
    memo = textwrap.dedent(memo).strip() + "\n"
    out = OUT_DIR / f"ho_proficiency_thresholds_memo_{DATE}.md"
    out.write_text(memo, encoding="utf-8")
    return out


def main() -> None:
    plt.rcParams.update(
        {
            "font.family": "DejaVu Sans",
            "axes.titlesize": 13,
            "axes.labelsize": 10.5,
            "figure.facecolor": "white",
            "axes.facecolor": "white",
        }
    )

    school_subject = build_school_subject_dataset()
    summary = subject_summary(school_subject)
    pairs = pair_examples(school_subject)
    corrs = correlation_sensitivity(school_subject)
    stress = stress_test_level_spacing(school_subject)

    school_subject.to_csv(OUT_DIR / f"ho_proficiency_thresholds_school_subject_{DATE}.csv", index=False)
    summary.to_csv(OUT_DIR / f"ho_proficiency_thresholds_summary_{DATE}.csv", index=False)
    pairs.to_csv(OUT_DIR / f"ho_proficiency_thresholds_same_pp_examples_{DATE}.csv", index=False)
    corrs.to_csv(OUT_DIR / f"ho_proficiency_thresholds_correlation_sensitivity_{DATE}.csv", index=False)
    stress.to_csv(OUT_DIR / f"ho_proficiency_thresholds_level_spacing_stress_test_{DATE}.csv", index=False)

    figures = [
        plot_statewide_levels(summary),
        plot_pp_vs_index(school_subject),
        plot_same_pp_pairs(pairs),
        plot_correlation_sensitivity(corrs),
    ]
    memo = write_memo(summary, pairs, corrs, stress, figures)

    print(f"Wrote {len(school_subject):,} school-subject rows")
    print(f"Wrote memo: {memo}")
    for figure in figures:
        print(f"Wrote figure: {figure}")


if __name__ == "__main__":
    main()
