1"""Replication entry point: run from workspace root with ~/wp29/venv/bin/python out/analysis.py"""2import json, os3import numpy as np, pandas as pd4import matplotlib5matplotlib.use("Agg")6import matplotlib.pyplot as plt7import statsmodels.formula.api as smf89SEED = None # no randomness used10OUT = "out"11d_raw = pd.read_csv("data/firms.csv")12n_raw = len(d_raw)1314# Exclusions15excl = []16d = d_raw.dropna()17excl.append(("1", "drop rows with missing values", n_raw, len(d), n_raw - len(d), "required fields"))18n0 = len(d); d = d[d.training_hours >= 0].copy()19excl.append(("2", "drop training_hours < 0", n0, len(d), n0 - len(d), "negative hours per employee are physically impossible (recording error)"))20pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{OUT}/exclusions.csv", index=False)2122COV = "HC1" # heteroskedasticity-robust SEs23def fit(formula, data):24 return smf.ols(formula, data).fit(cov_type=COV)2526def rmse(m):27 return float(np.sqrt(np.mean(m.resid ** 2)))2829naive = fit("productivity_growth ~ training_hours", d)30adj = fit("productivity_growth ~ training_hours + skill_index", d)31# sensitivity models32adj_all = fit("productivity_growth ~ training_hours + skill_index", d_raw)33naive_all = fit("productivity_growth ~ training_hours", d_raw)34adj_q = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d)3536def row(name, m):37 ci = m.conf_int().loc["training_hours"]38 return dict(model=name, estimate=float(m.params["training_hours"]), se=float(m.bse["training_hours"]),39 ci_low=float(ci[0]), ci_high=float(ci[1]), n=int(m.nobs), rmse=rmse(m))40rows = [row("unadjusted (main sample)", naive), row("adjusted for skill_index (main sample) [estimand]", adj),41 row("adjusted + skill_index^2 (sensitivity)", adj_q), row("unadjusted (all 1200 rows, sensitivity)", naive_all),42 row("adjusted for skill_index (all 1200 rows, sensitivity)", adj_all)]43pd.DataFrame(rows).to_csv(f"{OUT}/table1.csv", index=False)4445r_n, r_a = rows[0], rows[1]46lo, hi = r_a["ci_low"], r_a["ci_high"]47conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")48n_neg = int((d_raw.training_hours < 0).sum())49results = {50 "task_id": "F6-conf-04",51 "naive_estimate": r_n["estimate"], "naive_ci_low": r_n["ci_low"], "naive_ci_high": r_n["ci_high"],52 "naive_se": r_n["se"],53 "baseline_metric": r_n["rmse"], "model_metric": r_a["rmse"],54 "effect_estimate": r_a["estimate"], "effect_se": r_a["se"], "effect_se_type": "HC1 heteroskedasticity-robust",55 "effect_ci_low": lo, "effect_ci_high": hi,56 "skill_index_coefficient": float(adj.params["skill_index"]),57 "skill_index_se": float(adj.bse["skill_index"]),58 "conclusion": conclusion,59 "feasibility": "complete", "feasibility_reason": None,60 "data_issues": [{"issue": f"{n_neg} rows have negative training_hours (impossible values)",61 "handling": "excluded from main analysis; sensitivity estimate on all 1200 rows reported (sensitivity_all_rows_effect)"},62 {"issue": "strong confounding: skill_index correlates 0.69 with training_hours and -0.79 with productivity_growth",63 "handling": "adjusted for pre-programme skill_index by OLS"}],64 "n_input_rows": {"firms.csv": n_raw},65 "n_analysed": int(adj.nobs),66 "random_seed": SEED,67 "sensitivity_all_rows_effect": rows[4]["estimate"], "sensitivity_all_rows_se": rows[4]["se"],68 "sensitivity_quadratic_skill_effect": rows[2]["estimate"], "sensitivity_quadratic_skill_se": rows[2]["se"],69 "corr_skill_training": float(d.skill_index.corr(d.training_hours)),70}71with open(f"{OUT}/results.json", "w") as f:72 json.dump(results, f, indent=2)7374# Figure75fig, ax = plt.subplots(figsize=(7, 5))76ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.4, color="grey", label="firms (n=%d)" % len(d))77x = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)78ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, "b-", label=f"unadjusted OLS (slope {r_n['estimate']:.3f})")79sm = d.skill_index.mean()80ax.plot(x, adj.params["Intercept"] + adj.params["training_hours"] * x + adj.params["skill_index"] * sm, "r-",81 label=f"adjusted for skill_index at mean (slope {r_a['estimate']:.3f})")82ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth next year (%)")83ax.legend(); ax.set_title("Productivity growth vs training hours (simulated data)")84fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)85print(json.dumps({k: results[k] for k in ["naive_estimate", "effect_estimate", "effect_se", "effect_ci_low", "effect_ci_high", "conclusion", "n_analysed"]}, indent=1))86