1"""Replication entry point for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5import statsmodels.formula.api as smf6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = None # no randomness used11os.makedirs("out", exist_ok=True)12raw = pd.read_csv("data/firms.csv")13n_raw = len(raw)1415# ---- audit ----16issues = []17assert raw.isna().sum().sum() == 018assert raw.firm_id.duplicated().sum() == 0 and raw.duplicated().sum() == 019neg = raw.training_hours < 020issues.append({21 "issue": f"{int(neg.sum())} rows have negative training_hours (min {raw.training_hours.min():.3f}), which is physically impossible (hours per employee >= 0).",22 "handling": "Excluded from the main estimates; sensitivity analysis re-estimates both models on all rows (results essentially unchanged)."})23d = raw.loc[~neg].copy()2425def fit(formula, data):26 return smf.ols(formula, data).fit(cov_type="HC1")2728naive = fit("productivity_growth ~ training_hours", d)29adj = fit("productivity_growth ~ training_hours + skill_index", d)30naive_all = fit("productivity_growth ~ training_hours", raw)31adj_all = fit("productivity_growth ~ training_hours + skill_index", raw)32flex = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2) + I(training_hours**2) + training_hours:skill_index", d)33rmse = lambda m: float(np.sqrt(np.mean(m.resid ** 2)))3435def row(name, m):36 ci = m.conf_int().loc["training_hours"]37 return dict(model=name, estimate=m.params["training_hours"], se=m.bse["training_hours"],38 ci_low=ci[0], ci_high=ci[1])3940table = pd.DataFrame([41 row("unadjusted (main sample)", naive),42 row("adjusted for skill_index (main sample)", adj),43 row("unadjusted (all rows, sensitivity)", naive_all),44 row("adjusted for skill_index (all rows, sensitivity)", adj_all)])45table.to_csv("out/table1.csv", index=False)4647nci = naive.conf_int().loc["training_hours"]; aci = adj.conf_int().loc["training_hours"]48lo, hi = float(aci[0]), float(aci[1])49concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")50res = {51 "task_id": "F6-conf-04",52 "naive_estimate": float(naive.params["training_hours"]),53 "naive_ci_low": float(nci[0]), "naive_ci_high": float(nci[1]),54 "baseline_metric": rmse(naive), "model_metric": rmse(adj),55 "effect_estimate": float(adj.params["training_hours"]),56 "effect_se": float(adj.bse["training_hours"]),57 "effect_ci_low": lo, "effect_ci_high": hi,58 "skill_index_coefficient": float(adj.params["skill_index"]),59 "conclusion": concl,60 "feasibility": "complete", "feasibility_reason": None,61 "data_issues": issues,62 "n_input_rows": {"firms.csv": int(n_raw)},63 "n_analysed": int(adj.nobs),64 "random_seed": SEED,65 "se_type": "HC1 heteroskedasticity-robust",66 "corr_skill_training": float(d[["skill_index", "training_hours"]].corr().iloc[0, 1]),67 "sensitivity_all_rows_effect_estimate": float(adj_all.params["training_hours"]),68 "sensitivity_all_rows_naive_estimate": float(naive_all.params["training_hours"]),69 "diag_flexible_training_coef": float(flex.params["training_hours"]),70 "diag_flexible_min_pvalue_nonlinear_terms": float(min(flex.pvalues[["I(skill_index ** 2)", "I(training_hours ** 2)", "training_hours:skill_index"]])),71}72with open("out/results.json", "w") as f:73 json.dump(res, f, indent=2)7475# ---- figure ----76fig, ax = plt.subplots(figsize=(7, 5))77ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.4, color="grey", label="firms")78x = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)79ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, "b-",80 label=f"unadjusted OLS (slope {res['naive_estimate']:.3f})")81ax.plot(x, adj.params["Intercept"] + adj.params["training_hours"] * x + adj.params["skill_index"] * d.skill_index.mean(),82 "r-", label=f"adjusted for skill_index at its mean (slope {res['effect_estimate']:.3f})")83ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth next year (%)")84ax.set_title("Productivity growth vs training hours (simulated data)"); ax.legend()85fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)86print(json.dumps(res, indent=2)); print(table)87