1"""Replication script for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd7import statsmodels.formula.api as smf8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112SEED = None # no randomness used13os.makedirs("out", exist_ok=True)14raw = pd.read_csv("data/firms.csv")15n_raw = len(raw)1617issues = []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({"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}); "24 "hours per employee cannot be negative, so these values are invalid",25 "handling": "excluded from all main estimates; sensitivity analysis re-estimates with them included "26 "and with them set to 0 (reported in report/table1)"})27d = raw[~neg].copy()2829def fit(formula, data, cov="HC1"):30 return smf.ols(formula, data).fit(cov_type=cov)3132def rmse(m):33 return float(np.sqrt(np.mean(m.resid ** 2)))3435naive = fit("productivity_growth ~ training_hours", d)36adj = fit("productivity_growth ~ training_hours + skill_index", d)37# sensitivity38adj_all = fit("productivity_growth ~ training_hours + skill_index", raw)39rc = raw.copy(); rc["training_hours"] = rc.training_hours.clip(lower=0)40adj_clip = fit("productivity_growth ~ training_hours + skill_index", rc)41adj_conv = d.pipe(lambda x: smf.ols("productivity_growth ~ training_hours + skill_index", x).fit())42adj_quad = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d)4344def row(name, m):45 ci = m.conf_int().loc["training_hours"]46 return {"model": name, "estimate": m.params["training_hours"], "se": m.bse["training_hours"],47 "ci_low": ci[0], "ci_high": ci[1], "n": int(m.nobs), "rmse": rmse(m)}4849rows = [row("unadjusted (HC1)", naive), row("adjusted for skill_index (HC1) [main]", adj),50 row("adjusted, conventional SE", adj_conv),51 row("adjusted, incl. negative-hours rows (HC1)", adj_all),52 row("adjusted, negative hours set to 0 (HC1)", adj_clip),53 row("adjusted + skill_index^2 (HC1)", adj_quad)]54pd.DataFrame(rows).to_csv("out/table1.csv", index=False)5556nci = naive.conf_int().loc["training_hours"]; aci = adj.conf_int().loc["training_hours"]57lo, hi = float(aci[0]), float(aci[1])58concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")59res = {60 "task_id": "F6-conf-04",61 "naive_estimate": float(naive.params["training_hours"]),62 "naive_ci_low": float(nci[0]), "naive_ci_high": float(nci[1]),63 "baseline_metric": rmse(naive), "model_metric": rmse(adj),64 "effect_estimate": float(adj.params["training_hours"]),65 "effect_se": float(adj.bse["training_hours"]),66 "effect_se_type": "HC1 heteroskedasticity-robust",67 "effect_ci_low": lo, "effect_ci_high": hi,68 "skill_index_coefficient": float(adj.params["skill_index"]),69 "conclusion": concl,70 "feasibility": "complete", "feasibility_reason": None,71 "data_issues": issues,72 "n_input_rows": {"firms.csv": n_raw},73 "n_analysed": int(adj.nobs),74 "random_seed": SEED,75}76with open("out/results.json", "w") as f:77 json.dump(res, f, indent=2)7879fig, ax = plt.subplots(figsize=(7, 5))80ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.4, color="grey", label="firms (analysed)")81x = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)82ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, "r-",83 label=f"unadjusted OLS (slope {naive.params['training_hours']:.3f})")84sm = d.skill_index.mean()85ax.plot(x, adj.params["Intercept"] + adj.params["training_hours"] * x + adj.params["skill_index"] * sm, "b-",86 label=f"adjusted for skill_index at mean (slope {adj.params['training_hours']:.3f})")87ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth next year (%)")88ax.set_title("Productivity growth vs training hours (simulated data)")89ax.legend(fontsize=8); fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)90print(json.dumps(res, indent=2)); print(pd.DataFrame(rows).to_string())91