1"""Replication entry point for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads data/firms.csv only; writes out/results.json, out/table1.csv, out/figure1.png,4out/exclusions.csv, out/sensitivity.csv. Deterministic (no randomness used)."""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)15SEED = None # no stochastic step16COV = "HC1" # heteroskedasticity-robust standard errors1718raw = pd.read_csv("data/firms.csv")19n_raw = len(raw)20req = ["firm_id", "skill_index", "training_hours", "productivity_growth"]21assert all(c in raw.columns for c in req)2223excl = []24data_issues = []25d = raw.copy()2627def step(name, rule, mask_keep, rationale):28 global d29 nb = len(d); d = d[mask_keep(d)].copy(); na = len(d)30 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=rationale))3132step("1", "drop rows with missing values in required columns",33 lambda x: x[req].notna().all(axis=1), "complete cases required for OLS")34step("2", "drop duplicate firm_id",35 lambda x: ~x.firm_id.duplicated(), "one row per firm expected")36n_neg = int((d.training_hours < 0).sum())37step("3", "drop rows with training_hours < 0",38 lambda x: x.training_hours >= 0,39 "negative hours of training per employee are physically impossible; treated as recording errors")40if n_neg:41 data_issues.append({"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}); impossible values",42 "handling": "excluded from main analysis (exclusion step 3); sensitivity analyses keep them as-is and censor them to 0 (out/sensitivity.csv)"})43pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)44n_an = len(d)4546def fit(formula, df):47 return smf.ols(formula, data=df).fit(cov_type=COV)4849def rmse(m):50 return float(np.sqrt(np.mean(m.resid ** 2)))5152naive = fit("productivity_growth ~ training_hours", d)53adj = fit("productivity_growth ~ training_hours + skill_index", d)5455def row(name, m, n):56 ci = m.conf_int().loc["training_hours"]57 return dict(model=name, estimate=float(m.params["training_hours"]), se=float(m.bse["training_hours"]),58 ci_low=float(ci[0]), ci_high=float(ci[1]), n=int(n), rmse=rmse(m))5960table = pd.DataFrame([row("unadjusted: growth ~ hours", naive, n_an),61 row("adjusted: growth ~ hours + skill_index (target)", adj, n_an)])62table.to_csv(f"{OUT}/table1.csv", index=False)6364# Sensitivity analyses (not the main estimate)65sens = []66full = raw.dropna(subset=req)67sens.append(row("adjusted, all rows incl. negative hours", fit("productivity_growth ~ training_hours + skill_index", full), len(full)))68sens.append(row("unadjusted, all rows incl. negative hours", fit("productivity_growth ~ training_hours", full), len(full)))69cens = full.assign(training_hours=full.training_hours.clip(lower=0))70sens.append(row("adjusted, negative hours censored to 0", fit("productivity_growth ~ training_hours + skill_index", cens), len(cens)))71sens.append(row("adjusted + skill_index^2", fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d), n_an))72sens.append(row("adjusted + skill^2 + skill^3", fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2) + I(skill_index**3)", d), n_an))73sens.append(row("adjusted + hours^2 (linear term only reported)", fit("productivity_growth ~ training_hours + I(training_hours**2) + skill_index", d), n_an))74adj_conv = smf.ols("productivity_growth ~ training_hours + skill_index", data=d).fit()75sens.append(row("adjusted, conventional (non-robust) SE", adj_conv, n_an))76pd.DataFrame(sens).to_csv(f"{OUT}/sensitivity.csv", index=False)7778eff = adj.params["training_hours"]; lo, hi = adj.conf_int().loc["training_hours"]79conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")80nlo, nhi = naive.conf_int().loc["training_hours"]8182results = {83 "task_id": "F6-conf-04",84 "feasibility": "complete",85 "feasibility_reason": None,86 "data_issues": data_issues,87 "n_input_rows": {"firms.csv": int(n_raw)},88 "n_analysed": int(n_an),89 "random_seed": SEED,90 "naive_estimate": float(naive.params["training_hours"]),91 "naive_ci_low": float(nlo), "naive_ci_high": float(nhi),92 "naive_se": float(naive.bse["training_hours"]),93 "baseline_metric": rmse(naive),94 "model_metric": rmse(adj),95 "effect_estimate": float(eff),96 "effect_se": float(adj.bse["training_hours"]),97 "effect_ci_low": float(lo), "effect_ci_high": float(hi),98 "skill_index_coefficient": float(adj.params["skill_index"]),99 "skill_index_se": float(adj.bse["skill_index"]),100 "se_type": "HC1 heteroskedasticity-robust; 95% CI = estimate +/- 1.96*SE (normal)",101 "corr_hours_skill": float(d.training_hours.corr(d.skill_index)),102 "sens_effect_all_rows": sens[0]["estimate"],103 "sens_effect_censored0": sens[2]["estimate"],104 "sens_effect_skill_sq": sens[3]["estimate"],105 "sens_effect_skill_cubic": sens[4]["estimate"],106 "sens_effect_conventional_se": sens[6]["se"],107 "conclusion": conclusion,108}109with open(f"{OUT}/results.json", "w") as f:110 json.dump(results, f, indent=2)111112# Figure113fig, ax = plt.subplots(figsize=(7, 5))114ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.4, color="grey", label=f"firms (n={n_an})")115xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)116ax.plot(xs, naive.params["Intercept"] + naive.params["training_hours"] * xs, "b-",117 label=f"unadjusted: slope {naive.params['training_hours']:.3f}")118sbar = d.skill_index.mean()119ax.plot(xs, adj.params["Intercept"] + adj.params["skill_index"] * sbar + eff * xs, "r-",120 label=f"adjusted for skill_index (at mean): slope {eff:.3f}")121ax.set_xlabel("training hours per employee (programme year)")122ax.set_ylabel("productivity growth, following year (%)")123ax.set_title("Simulated firms: productivity growth vs training hours")124ax.legend(fontsize=8)125fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)126127print(table.to_string(index=False))128print(pd.DataFrame(sens).to_string(index=False))129print("conclusion:", conclusion)130