1"""Replication entry point for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/firms.csv; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.formula.api as smf9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"14os.makedirs(OUT, exist_ok=True)15raw = pd.read_csv("data/firms.csv")16n_raw = len(raw)1718# --- exclusions ---19excl = []20d = raw.dropna()21excl.append(("1", "drop rows with missing values", n_raw, len(d), "none expected; complete-case"))22n0 = len(d)23d = d[d["training_hours"] >= 0].copy()24excl.append(("2", "drop training_hours < 0", n0, len(d), "negative hours of training are physically impossible"))25pd.DataFrame([dict(step=s, rule=r, n_before=a, n_after=b, n_removed=a - b, rationale=why)26 for s, r, a, b, why in excl]).to_csv(f"{OUT}/exclusions.csv", index=False)27n_neg = int((raw["training_hours"] < 0).sum())2829COV = "HC1" # heteroskedasticity-robust SEs3031def fit(formula, data):32 m = smf.ols(formula, data).fit(cov_type=COV)33 ci = m.conf_int(0.05)34 rmse = float(np.sqrt(np.mean(m.resid ** 2)))35 return m, ci, rmse3637naive, naive_ci, naive_rmse = fit("productivity_growth ~ training_hours", d)38adj, adj_ci, adj_rmse = fit("productivity_growth ~ training_hours + skill_index", d)39# sensitivity analyses (not the main estimate)40naive_all, _, _ = fit("productivity_growth ~ training_hours", raw)41sens_all, sens_all_ci, _ = fit("productivity_growth ~ training_hours + skill_index", raw)42sens_int, sens_int_ci, _ = fit("productivity_growth ~ training_hours * skill_index", d)43sens_q, sens_q_ci, _ = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d)44adj_conv = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit()4546b = float(adj.params["training_hours"]); lo, hi = map(float, adj_ci.loc["training_hours"])47conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")4849rows = []50def addrow(name, m, ci):51 rows.append(dict(model=name, estimate=float(m.params["training_hours"]), se=float(m.bse["training_hours"]),52 ci_low=float(ci.loc["training_hours", 0]), ci_high=float(ci.loc["training_hours", 1])))53addrow("unadjusted (growth ~ hours)", naive, naive_ci)54addrow("adjusted (growth ~ hours + skill) [main]", adj, adj_ci)55addrow("sensitivity: adjusted, incl. negative-hours rows", sens_all, sens_all_ci)56addrow("sensitivity: adjusted + hours x skill (coef at skill=0)", sens_int, sens_int_ci)57addrow("sensitivity: adjusted + skill^2", sens_q, sens_q_ci)58pd.DataFrame(rows).to_csv(f"{OUT}/table1.csv", index=False)5960results = {61 "task_id": "F6-conf-04",62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": [65 {"issue": f"{n_neg} rows have negative training_hours (physically impossible; all low skill_index)",66 "handling": "excluded from main estimate; sensitivity estimate including them reported in table1.csv"},67 ],68 "n_input_rows": {"firms.csv": int(n_raw)},69 "n_analysed": int(len(d)),70 "random_seed": None,71 "naive_estimate": float(naive.params["training_hours"]),72 "naive_ci_low": float(naive_ci.loc["training_hours", 0]),73 "naive_ci_high": float(naive_ci.loc["training_hours", 1]),74 "baseline_metric": naive_rmse,75 "model_metric": adj_rmse,76 "effect_estimate": b,77 "effect_se": float(adj.bse["training_hours"]),78 "effect_ci_low": lo,79 "effect_ci_high": hi,80 "skill_index_coefficient": float(adj.params["skill_index"]),81 "conclusion": conclusion,82 "se_type": "HC1 heteroskedasticity-robust",83 "effect_se_conventional": float(adj_conv.bse["training_hours"]),84 "naive_incl_negative_hours_estimate": float(naive_all.params["training_hours"]),85 "sensitivity_incl_negative_hours_estimate": float(sens_all.params["training_hours"]),86 "sensitivity_interaction_term": float(sens_int.params["training_hours:skill_index"]),87 "sensitivity_interaction_term_p": float(sens_int.pvalues["training_hours:skill_index"]),88 "sensitivity_skill_sq_estimate": float(sens_q.params["training_hours"]),89 "corr_skill_training": float(d["skill_index"].corr(d["training_hours"])),90 "adjusted_r2": float(adj.rsquared),91 "naive_r2": float(naive.rsquared),92}93with open(f"{OUT}/results.json", "w") as f:94 json.dump(results, f, indent=2)9596# --- figure ---97fig, ax = plt.subplots(figsize=(7, 5))98ax.scatter(d["training_hours"], d["productivity_growth"], s=8, alpha=0.35, color="grey", label="firms")99x = np.linspace(d["training_hours"].min(), d["training_hours"].max(), 100)100ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, color="tab:red",101 label=f"unadjusted: slope {results['naive_estimate']:.3f}")102sm = d["skill_index"].mean()103ax.plot(x, adj.params["Intercept"] + b * x + adj.params["skill_index"] * sm, color="tab:blue",104 label=f"adjusted for skill (at mean skill): slope {b:.3f}")105ax.set_xlabel("training hours per employee (programme year)")106ax.set_ylabel("productivity growth, following year (%)")107ax.set_title("Productivity growth vs training hours (simulated data)")108ax.legend()109fig.tight_layout()110fig.savefig(f"{OUT}/figure1.png", dpi=120)111112print(json.dumps({k: results[k] for k in ["naive_estimate", "effect_estimate", "effect_se",113 "effect_ci_low", "effect_ci_high", "conclusion", "n_analysed"]}, indent=1))114