1"""Replication entry point for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Deterministic (no randomness used)."""4import json, os5import numpy as np, pandas as pd6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt9import statsmodels.formula.api as smf1011OUT = "out"12os.makedirs(OUT, exist_ok=True)13raw = pd.read_csv("data/firms.csv")14n_raw = len(raw)1516# ---- data audit ----17issues = []18assert raw.isna().sum().sum() == 019assert not raw.firm_id.duplicated().any()20neg = raw.training_hours < 021n_neg = int(neg.sum())22if n_neg:23 issues.append({24 "issue": f"{n_neg} firms have negative training_hours (min {raw.training_hours.min():.3f}); "25 "hours of training per employee cannot be negative, so these values are invalid.",26 "handling": "Excluded from the main analysis (no imputation). Sensitivity: re-estimated with them "27 "included; the training coefficient changes by <0.002 pp/hour."})28d = raw.loc[~neg].copy()2930def rmse(res):31 return float(np.sqrt(np.mean(res.resid ** 2)))3233COV = "HC1" # heteroskedasticity-robust SEs throughout34naive = smf.ols("productivity_growth ~ training_hours", d).fit(cov_type=COV)35adj = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit(cov_type=COV)36# sensitivity checks37adj_all = smf.ols("productivity_growth ~ training_hours + skill_index", raw).fit(cov_type=COV)38adj_q = smf.ols("productivity_growth ~ training_hours + skill_index + I(skill_index**2) + I(training_hours**2)",39 d).fit(cov_type=COV)40adj_conv = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit()4142def row(name, m, var="training_hours"):43 ci = m.conf_int().loc[var]44 return {"model": name, "estimate": float(m.params[var]), "se": float(m.bse[var]),45 "ci_low": float(ci[0]), "ci_high": float(ci[1])}4647rows = [row("unadjusted OLS (training_hours only)", naive),48 row("adjusted OLS (training_hours + skill_index) [main]", adj),49 row("sensitivity: adjusted, incl. negative-hours rows", adj_all),50 row("sensitivity: adjusted + quadratic terms (linear coef)", adj_q),51 row("sensitivity: adjusted, conventional SE", adj_conv)]52pd.DataFrame(rows).to_csv(f"{OUT}/table1.csv", index=False)53print(pd.DataFrame(rows).to_string())54print("quadratic-term p-values:", adj_q.pvalues.to_dict())5556# ---- figure ----57fig, ax = plt.subplots(figsize=(7, 5))58ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.4, color="grey", label="firms")59x = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)60ax.plot(x, naive.params.Intercept + naive.params.training_hours * x, color="C0",61 label=f"unadjusted OLS (slope {naive.params.training_hours:.3f})")62sm = d.skill_index.mean()63ax.plot(x, adj.params.Intercept + adj.params.training_hours * x + adj.params.skill_index * sm, color="C3",64 label=f"adjusted for skill_index at mean (slope {adj.params.training_hours:.3f})")65ax.set_xlabel("training hours per employee (programme year)")66ax.set_ylabel("productivity growth, following year (%)")67ax.set_title("Productivity growth vs training hours (simulated data, n=%d)" % len(d))68ax.legend(fontsize=8)69fig.tight_layout()70fig.savefig(f"{OUT}/figure1.png", dpi=150)7172lo, hi = rows[1]["ci_low"], rows[1]["ci_high"]73concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")74res = {75 "task_id": "F6-conf-04",76 "naive_estimate": rows[0]["estimate"], "naive_ci_low": rows[0]["ci_low"], "naive_ci_high": rows[0]["ci_high"],77 "baseline_metric": rmse(naive), "model_metric": rmse(adj),78 "effect_estimate": rows[1]["estimate"], "effect_se": rows[1]["se"],79 "effect_ci_low": lo, "effect_ci_high": hi,80 "skill_index_coefficient": float(adj.params.skill_index),81 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,82 "data_issues": issues, "n_input_rows": {"firms.csv": n_raw},83 "n_analysed": int(adj.nobs), "random_seed": None,84}85with open(f"{OUT}/results.json", "w") as f:86 json.dump(res, f, indent=2)87print(json.dumps(res, indent=2))88