"""
Statistics coursework sample: hypothesis tests and regression on real data
Python script behind https://assignmenthelper.org/statistics-coursework-sample-r-excel/

What this is: the script that computed every figure on that page. Its output
is published beside it as statistics-coursework-student-por.out.txt, and the
R route (statistics-coursework-student-por.R, with its own output file) was
run afterwards and agrees at every figure the page prints.

Dataset : Student Performance, UCI Machine Learning Repository, id 320.
          Cortez, P. (2008). https://doi.org/10.24432/C5TG7T
          https://archive.ics.uci.edu/dataset/320/student+performance
Download: https://archive.ics.uci.edu/static/public/320/student+performance.zip
          (40,735 bytes; contains student.zip, which contains student-por.csv,
          student-mat.csv, student.txt and student-merge.R)
Licence : Creative Commons Attribution 4.0 International (CC BY 4.0).
File    : student-por.csv (649 rows, 33 columns, semicolon-separated).
Run date: 2026-09-25 (first run 2026-09-24; sections 6b and the extra checks
          in 2, 3 and 5 were added on the 25th).
Versions: Python 3.13.15, pandas 3.0.6, scipy 1.18.1, statsmodels 0.15.0,
          numpy 2.5.3.

Usage
    uv run --with pandas --with scipy --with statsmodels \
        python statistics-coursework-student-por.py /path/to/student-por.csv

The csv is not bundled (attribution belongs at the source and the download
URL is stable); download the zip above, unzip both layers and pass the path.
Nothing is random, so there is no seed.

Sections match the page's H2s:
  1. Descriptives (G3, absences, studytime bands, failures, the G3 zero spike)
  2. Welch two-sample t-test, G3 by internet: Levene first, Welch t, mean
     difference with 95% CI, Cohen's d with an approximate 95% CI, then the
     within-group Shapiro-Wilk and a Mann-Whitney check
  3. One-way ANOVA, G3 by studytime band: Levene, ANOVA table, eta squared,
     Tukey HSD, Welch ANOVA as the fallback if Levene rejects, Kruskal-Wallis
  4. Chi-square test of independence, school x higher: observed, row
     percentages, expected counts, Pearson chi-square with and without the
     Yates correction (R applies it on 2x2 by default; Excel does not),
     Cramer's V
  5. Multiple linear regression G3 ~ studytime + failures + absences + higher
     + school + sex + Medu, with 95% CIs, R2, adjusted R2, F, residual SE;
     the unadjusted slope of G3 on studytime; internet added to the model
  6. Diagnostics: residual summary, Shapiro-Wilk, Breusch-Pagan, Ramsey RESET,
     VIF, Cook's distance, the zero-grade profile; then re-run 1 (HC3 standard
     errors) and re-run 2 (without the zero-grade rows) with the checks
     repeated on re-run 2
  6b. Locating the curvature RESET found: study time and failures each freed
     into categories, the nested F-tests, RESET on each, failures squared,
     the raw means by failures, and re-run 3 (failures as categories)
  7. Leakage: the same model with G1 and G2 added
"""
from __future__ import annotations

import sys

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats
from statsmodels.stats.anova import anova_lm
from statsmodels.stats.diagnostic import het_breuschpagan, linear_reset
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.oneway import anova_oneway
from statsmodels.stats.outliers_influence import variance_inflation_factor

pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 30)
pd.set_option("display.float_format", lambda v: f"{v:,.3f}")

FORMULA = "G3 ~ studytime + failures + absences + higher_yes + school_MS + sex_M + Medu"
LEAK_FORMULA = FORMULA + " + G1 + G2"
PREDICTORS = ["studytime", "failures", "absences", "higher_yes", "school_MS", "sex_M", "Medu"]
BAND = {1: "under 2 h", 2: "2 to 5 h", 3: "5 to 10 h", 4: "over 10 h"}


def section(title: str) -> None:
    print("\n" + "=" * 78)
    print(title)
    print("=" * 78)


def p_str(p: float) -> str:
    return "p < 0.001" if p < 0.001 else f"p = {p:.3f}"


def cohens_d(a: np.ndarray, b: np.ndarray) -> tuple[float, float, float]:
    """Cohen's d with the pooled SD; approximate 95% CI from the large-sample
    standard error sqrt((n1+n2)/(n1 n2) + d^2 / (2 (n1+n2)))."""
    n1, n2 = len(a), len(b)
    s_pooled = np.sqrt(((n1 - 1) * a.var(ddof=1) + (n2 - 1) * b.var(ddof=1)) / (n1 + n2 - 2))
    d = (a.mean() - b.mean()) / s_pooled
    se = np.sqrt((n1 + n2) / (n1 * n2) + d**2 / (2 * (n1 + n2)))
    return d, d - 1.96 * se, d + 1.96 * se


def coef_table(res, label: str) -> None:
    ci = res.conf_int()
    out = pd.DataFrame(
        {
            "coef": res.params,
            "SE": res.bse,
            "t": res.tvalues,
            "p": res.pvalues,
            "CI low": ci[0],
            "CI high": ci[1],
        }
    )
    print(f"\n{label}")
    print(out.to_string())
    print(
        f"n = {int(res.nobs)}; R2 = {res.rsquared:.3f}; adjusted R2 = {res.rsquared_adj:.3f}; "
        f"F({int(res.df_model)}, {int(res.df_resid)}) = {res.fvalue:.2f}, {p_str(res.f_pvalue)}; "
        f"residual SE = {np.sqrt(res.mse_resid):.3f}"
    )


def main(path: str) -> None:
    d = pd.read_csv(path, sep=";")
    d["higher_yes"] = (d.higher == "yes").astype(int)
    d["school_MS"] = (d.school == "MS").astype(int)
    d["sex_M"] = (d.sex == "M").astype(int)

    # ---------------------------------------------------------- 1. describe
    section("1. Descriptives (student-por.csv)")
    print(f"rows: {len(d)}; columns: {d.shape[1] - 3} in the file (+3 dummies built here)")
    print(f"missing values in the file: {int(d.isna().sum().sum())}")
    desc = d[["G3", "G1", "G2", "absences", "age"]].describe(percentiles=[0.5]).T
    desc = desc.rename(columns={"50%": "median"})
    print("\nNumeric summary:")
    print(desc[["count", "mean", "std", "median", "min", "max"]].to_string())
    print(f"\nG3 = 0: {int((d.G3 == 0).sum())} students ({(d.G3 == 0).mean():.1%})")
    print("\nG3 count table (0 to 20):")
    print(d.G3.value_counts().sort_index().to_string())
    print("\nstudytime bands (1: under 2 h, 2: 2 to 5 h, 3: 5 to 10 h, 4: over 10 h):")
    print(d.studytime.value_counts().sort_index().to_string())
    print("\nfailures (past class failures):")
    print(d.failures.value_counts().sort_index().to_string())
    for col in ["internet", "higher", "school", "sex"]:
        print(f"\n{col}:")
        print(d[col].value_counts().to_string())
    print("\nMedu (mother's education, 0 none to 4 higher):")
    print(d.Medu.value_counts().sort_index().to_string())
    print(f"\nabsences: median {d.absences.median():.0f}, share with 0 absences {(d.absences == 0).mean():.1%}, max {d.absences.max()}")

    # -------------------------------------------- 2. Welch t-test, internet
    section("2. Two-sample t-test: G3 by internet access at home")
    yes = d.loc[d.internet == "yes", "G3"].to_numpy(dtype=float)
    no = d.loc[d.internet == "no", "G3"].to_numpy(dtype=float)
    grp = d.groupby("internet").G3.agg(n="size", mean="mean", sd="std", median="median")
    print(grp.to_string())
    lev = stats.levene(yes, no, center="median")
    print(f"\nLevene (Brown-Forsythe, centre = median, as car::leveneTest default): W = {lev.statistic:.3f}, {p_str(lev.pvalue)}")
    welch = stats.ttest_ind(yes, no, equal_var=False)
    ci = welch.confidence_interval(0.95)
    print(f"Welch t = {welch.statistic:.3f}, df = {welch.df:.1f}, {p_str(welch.pvalue)} (two-sided)")
    print(f"mean difference (yes minus no) = {yes.mean() - no.mean():.3f}, 95% CI [{ci.low:.3f}, {ci.high:.3f}]")
    student = stats.ttest_ind(yes, no, equal_var=True)
    print(f"Student (pooled) t = {student.statistic:.3f}, df = {student.df:.0f}, {p_str(student.pvalue)} (for the Excel equal-variance route)")
    dd, dlo, dhi = cohens_d(yes, no)
    print(f"Cohen's d = {dd:.3f}, approximate 95% CI [{dlo:.3f}, {dhi:.3f}]")
    sw_yes, sw_no = stats.shapiro(yes), stats.shapiro(no)
    print(f"Shapiro-Wilk within each group (the zero grades make both non-normal): yes W = {sw_yes.statistic:.3f}, {p_str(sw_yes.pvalue)}; no W = {sw_no.statistic:.3f}, {p_str(sw_no.pvalue)}")
    mw = stats.mannwhitneyu(yes, no, alternative="two-sided")
    print(f"Mann-Whitney U (rank-based check, no normality assumption) = {mw.statistic:,.0f}, {p_str(mw.pvalue)}; R's wilcox.test prints the other group's U, {len(yes) * len(no) - mw.statistic:,.0f}")

    # --------------------------------------------- 3. one-way ANOVA, studytime
    section("3. One-way ANOVA: G3 by weekly study-time band")
    bands = d.groupby("studytime").G3.agg(n="size", mean="mean", sd="std", median="median")
    bands.index = [f"{k} ({BAND[k]})" for k in bands.index]
    print(bands.to_string())
    groups = [g.G3.to_numpy(dtype=float) for _, g in d.groupby("studytime")]
    lev = stats.levene(*groups, center="median")
    print(f"\nLevene (centre = median): W = {lev.statistic:.3f}, {p_str(lev.pvalue)}")
    ols_band = smf.ols("G3 ~ C(studytime)", data=d).fit()
    table = anova_lm(ols_band, typ=1)
    table.index = ["between (studytime)", "within (residual)"]
    print("\nANOVA table (type I, one factor):")
    print(table.to_string())
    ss_b, ss_w = table.loc["between (studytime)", "sum_sq"], table.loc["within (residual)", "sum_sq"]
    f_val, f_p = table.loc["between (studytime)", "F"], table.loc["between (studytime)", "PR(>F)"]
    print(f"F({int(table.df.iloc[0])}, {int(table.df.iloc[1])}) = {f_val:.2f}, {p_str(f_p)}")
    print(f"mean squares to two decimals: between {table.loc['between (studytime)', 'mean_sq']:.2f}, within {table.loc['within (residual)', 'mean_sq']:.2f}; total SS {ss_b + ss_w:.2f} on {int(table.df.sum())} df")
    print(f"eta squared = SS_between / SS_total = {ss_b:.2f} / {ss_b + ss_w:.2f} = {ss_b / (ss_b + ss_w):.3f}")
    welch_a = anova_oneway(d.G3, d.studytime, use_var="unequal", welch_correction=True)
    print(f"Welch ANOVA (fallback if Levene rejects): F({welch_a.df_num:.0f}, {welch_a.df_denom:.1f}) = {welch_a.statistic:.2f}, {p_str(welch_a.pvalue)}")
    kw = stats.kruskal(*groups)
    print(f"Kruskal-Wallis H(3) (rank-based check, no normality assumption) = {kw.statistic:.2f}, {p_str(kw.pvalue)}")
    tukey = pairwise_tukeyhsd(d.G3, d.studytime.map(lambda k: f"band {k}"), alpha=0.05)
    print("\nTukey HSD (alpha = 0.05):")
    print(tukey.summary())

    # --------------------------------------------- 4. chi-square, school x higher
    section("4. Chi-square test of independence: school x higher")
    obs = pd.crosstab(d.school, d.higher)
    print("Observed:")
    print(obs.to_string())
    print("\nRow percentages (share of each school):")
    print((obs.div(obs.sum(axis=1), axis=0) * 100).round(1).to_string())
    chi_y, p_y, dof, expected = stats.chi2_contingency(obs, correction=True)
    chi_n, p_n, _, _ = stats.chi2_contingency(obs, correction=False)
    print("\nExpected counts under independence (row total x column total / N):")
    print(pd.DataFrame(expected, index=obs.index, columns=obs.columns).round(2).to_string())
    print(f"\nPearson chi-square with Yates continuity correction (R chisq.test default on 2x2): chi2({dof}) = {chi_y:.3f}, {p_str(p_y)}")
    print(f"Pearson chi-square without correction (Excel CHISQ.TEST): chi2({dof}) = {chi_n:.3f}, {p_str(p_n)}")
    n = obs.to_numpy().sum()
    v = np.sqrt(chi_n / (n * (min(obs.shape) - 1)))
    print(f"Cramer's V (from the uncorrected statistic) = {v:.3f}")
    print(f"minimum expected count = {expected.min():.2f} (rule: every expected count at least 5)")

    # --------------------------------------------- 5. multiple regression
    section("5. Multiple linear regression (model fixed before looking at the output)")
    print(FORMULA)
    print("dummies: higher_yes = 1 if higher == yes; school_MS = 1 if school == MS; sex_M = 1 if sex == M")
    res = smf.ols(FORMULA, data=d).fit()
    coef_table(res, "OLS coefficients with conventional standard errors:")
    res_raw = smf.ols("G3 ~ studytime", data=d).fit()
    ci_raw = res_raw.conf_int().loc["studytime"]
    print(f"\nUnadjusted slope of G3 on studytime alone (for the adjusted-versus-raw comparison): {res_raw.params['studytime']:.3f} per band, 95% CI [{ci_raw[0]:.3f}, {ci_raw[1]:.3f}], {p_str(res_raw.pvalues['studytime'])}, R2 = {res_raw.rsquared:.3f}")
    d["internet_yes"] = (d.internet == "yes").astype(int)
    res_int = smf.ols(FORMULA + " + internet_yes", data=d).fit()
    print(f"Internet added to the seven-predictor model (ties test 1 to the model): internet_yes = {res_int.params['internet_yes']:.3f}, {p_str(res_int.pvalues['internet_yes'])}; studytime = {res_int.params['studytime']:.3f}, {p_str(res_int.pvalues['studytime'])}")

    # --------------------------------------------- 6. diagnostics
    section("6. Diagnostics on the model in section 5")
    resid = res.resid
    print(f"residuals: mean {resid.mean():.4f}, SD {resid.std(ddof=1):.3f}, min {resid.min():.3f}, max {resid.max():.3f}, skewness {stats.skew(resid):.3f}, excess kurtosis {stats.kurtosis(resid):.3f}")
    std_resid = res.get_influence().resid_studentized_internal
    print(f"standardised residuals below -3: {int((std_resid < -3).sum())}; above +3: {int((std_resid > 3).sum())}")
    print(f"of the residuals below -3, students with G3 = 0: {int(((std_resid < -3) & (d.G3 == 0)).sum())}")
    sw = stats.shapiro(resid)
    print(f"Shapiro-Wilk on residuals: W = {sw.statistic:.3f}, {p_str(sw.pvalue)}")
    bp = het_breuschpagan(resid, res.model.exog)
    print(f"Breusch-Pagan (studentised, Koenker): LM = {bp[0]:.2f}, {p_str(bp[1])}; F = {bp[2]:.2f}, {p_str(bp[3])}")
    reset = linear_reset(res, power=3, test_type="fitted", use_f=True)
    print(f"Ramsey RESET (squared and cubed fitted values, as lmtest::resettest default): F({int(reset.df_num)}, {int(reset.df_denom)}) = {float(reset.fvalue):.2f}, {p_str(float(reset.pvalue))}")
    X = sm.add_constant(d[PREDICTORS])
    vif = pd.Series({col: variance_inflation_factor(X.values, i) for i, col in enumerate(X.columns) if col != "const"})
    print("\nVIF:")
    print(vif.round(3).to_string())
    cooks = res.get_influence().cooks_distance[0]
    print(f"\nCook's distance: max {cooks.max():.4f}; count above 4/n = {4 / len(d):.4f}: {int((cooks > 4 / len(d)).sum())}; count above 1: {int((cooks > 1).sum())}")
    top = np.argsort(cooks)[::-1][:5]
    print("five largest Cook's distances (row index, G3, studytime, failures, absences, Cook's D):")
    for i in top:
        print(f"  row {i}: G3 = {d.G3.iloc[i]}, studytime = {d.studytime.iloc[i]}, failures = {d.failures.iloc[i]}, absences = {d.absences.iloc[i]}, D = {cooks[i]:.4f}")

    z = d[d.G3 == 0]
    print(f"\nThe {len(z)} zero-grade students: {int((z.absences == 0).sum())} of {len(z)} record 0 absences; {int((z.school == 'MS').sum())} of {len(z)} are at MS; {int((z.G2 == 0).sum())} of {len(z)} already had G2 = 0; G1 runs from {int(z.G1.min())} to {int(z.G1.max())}")

    res_hc3 = smf.ols(FORMULA, data=d).fit(cov_type="HC3", use_t=True)
    coef_table(res_hc3, "Re-run 1: same model, HC3 heteroscedasticity-consistent standard errors (t-based intervals, as R's coefci):")

    d_nz = d[d.G3 > 0]
    res_nz = smf.ols(FORMULA, data=d_nz).fit()
    coef_table(res_nz, f"Re-run 2: same model without the {int((d.G3 == 0).sum())} students with G3 = 0 (n = {len(d_nz)}):")
    resid_nz = res_nz.resid
    std_nz = res_nz.get_influence().resid_studentized_internal
    sw_nz = stats.shapiro(resid_nz)
    bp_nz = het_breuschpagan(resid_nz, res_nz.model.exog)
    print(f"re-run 2 residuals: skewness {stats.skew(resid_nz):.3f}, excess kurtosis {stats.kurtosis(resid_nz):.3f}; standardised residuals below -3: {int((std_nz < -3).sum())}, above +3: {int((std_nz > 3).sum())}")
    print(f"re-run 2 Shapiro-Wilk: W = {sw_nz.statistic:.3f}, {p_str(sw_nz.pvalue)}; Breusch-Pagan (studentised): LM = {bp_nz[0]:.2f}, {p_str(bp_nz[1])}")
    reset_nz = linear_reset(res_nz, power=3, test_type="fitted", use_f=True)
    print(f"re-run 2 RESET: F({int(reset_nz.df_num)}, {int(reset_nz.df_denom)}) = {float(reset_nz.fvalue):.2f}, {p_str(float(reset_nz.pvalue))}")

    # ------------------------------ 6b. locating the curvature the RESET test found
    section("6b. Which term bends the line? Nested F-tests and RESET on re-coded models")
    print("Each candidate term is freed from its linear coding into categories; the nested F asks whether the")
    print("categories fit better than the linear coding, and RESET is repeated on the re-coded model.")

    def reset_line(r) -> str:
        t = linear_reset(r, power=3, test_type="fitted", use_f=True)
        return f"RESET F({int(t.df_num)}, {int(t.df_denom)}) = {float(t.fvalue):.2f}, {p_str(float(t.pvalue))}"

    res_st_cat = smf.ols(FORMULA.replace("studytime", "C(studytime)"), data=d).fit()
    nest_st = anova_lm(res, res_st_cat)
    print(f"\nstudytime as four categories instead of a 1 to 4 scale: nested F({int(nest_st.df_diff.iloc[1])}, {int(nest_st.df_resid.iloc[1])}) = {nest_st.F.iloc[1]:.2f}, {p_str(nest_st['Pr(>F)'].iloc[1])} (equal steps not rejected); {reset_line(res_st_cat)} (curvature unchanged)")

    res_f_cat = smf.ols(FORMULA.replace("failures", "C(failures)"), data=d).fit()
    nest_f = anova_lm(res, res_f_cat)
    print(f"failures as categories instead of a count: nested F({int(nest_f.df_diff.iloc[1])}, {int(nest_f.df_resid.iloc[1])}) = {nest_f.F.iloc[1]:.2f}, {p_str(nest_f['Pr(>F)'].iloc[1])} (the count coding is rejected); {reset_line(res_f_cat)} (curvature gone)")

    res_f_sq = smf.ols(FORMULA + " + I(failures ** 2)", data=d).fit()
    print(f"failures squared added to the original model: coefficient {res_f_sq.params['I(failures ** 2)']:.3f}, {p_str(res_f_sq.pvalues['I(failures ** 2)'])}; {reset_line(res_f_sq)}")

    print("\nRaw mean G3 by number of past failures (the whole drop comes with the first failure):")
    print(d.groupby("failures").G3.agg(n="size", mean="mean").to_string())

    coef_table(res_f_cat, "Re-run 3: same model with failures as categories (0 failures is the baseline):")
    print(reset_line(res_f_cat))

    # --------------------------------------------- 7. leakage
    section("7. Leakage: adding the period grades G1 and G2 as predictors")
    print(f"correlation G1 with G3: r = {d.G1.corr(d.G3):.3f}; G2 with G3: r = {d.G2.corr(d.G3):.3f}; G1 with G2: r = {d.G1.corr(d.G2):.3f}")
    res_leak = smf.ols(LEAK_FORMULA, data=d).fit()
    coef_table(res_leak, "OLS with G1 and G2 added:")
    print("\nSide by side (coefficient, p) for the seven original predictors:")
    side = pd.DataFrame(
        {
            "coef without G1, G2": res.params[PREDICTORS],
            "p without": res.pvalues[PREDICTORS],
            "coef with G1, G2": res_leak.params[PREDICTORS],
            "p with": res_leak.pvalues[PREDICTORS],
        }
    )
    print(side.to_string())
    print(f"adjusted R2: {res.rsquared_adj:.3f} without the period grades, {res_leak.rsquared_adj:.3f} with them")


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else "student-por.csv")