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.copy()21def step(name, rule, mask, why):22 global d23 nb = len(d); d = d[mask(d)].copy(); na = len(d)24 excl.append(dict(step=name, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))25step("1", "drop rows with missing values in any column", lambda x: x.notna().all(axis=1),26 "complete cases required (none expected)")27step("2", "drop duplicated firm_id", lambda x: ~x.firm_id.duplicated(), "one row per firm")28step("3", "drop training_hours < 0", lambda x: x.training_hours >= 0,29 "negative hours of training are logically impossible; treated as invalid treatment records")30pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)31n_neg = int((raw.training_hours < 0).sum())3233# ---- models -----------------------------------------------------------------34COV = "HC1" # heteroskedasticity-robust SEs35def fit(formula, data):36 m = smf.ols(formula, data).fit(cov_type=COV)37 ci = m.conf_int(0.05).loc["training_hours"]38 rmse = float(np.sqrt(np.mean(m.resid ** 2)))39 return m, float(m.params["training_hours"]), float(m.bse["training_hours"]), float(ci[0]), float(ci[1]), rmse4041m0, b0, se0, l0, h0, rmse0 = fit("productivity_growth ~ training_hours", d)42m1, b1, se1, l1, h1, rmse1 = fit("productivity_growth ~ training_hours + skill_index", d)43# sensitivity: all raw rows (including negative hours)44_, bs, ses, ls, hs, _ = fit("productivity_growth ~ training_hours + skill_index", raw)45_, bn_all, sen_all, ln_all, hn_all, _ = fit("productivity_growth ~ training_hours", raw)46# diagnostics (not the estimand): quadratic in skill / hours, interaction47_, bq, seq, lq, hq, rmseq = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d)48mi = smf.ols("productivity_growth ~ training_hours * skill_index", d).fit(cov_type=COV)49mh = smf.ols("productivity_growth ~ training_hours + I(training_hours**2) + skill_index", d).fit(cov_type=COV)50mconv = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit()5152if l1 > 0: concl = "positive_effect"53elif h1 < 0: concl = "negative_effect"54else: concl = "no_detectable_effect"5556table = pd.DataFrame([57 dict(model="unadjusted OLS (growth ~ hours), HC1", estimate=b0, se=se0, ci_low=l0, ci_high=h0),58 dict(model="adjusted OLS (growth ~ hours + skill_index), HC1 [main]", estimate=b1, se=se1, ci_low=l1, ci_high=h1),59 dict(model="sensitivity: adjusted OLS incl. 12 negative-hours rows, HC1", estimate=bs, se=ses, ci_low=ls, ci_high=hs),60 dict(model="sensitivity: adjusted OLS + skill_index^2, HC1", estimate=bq, se=seq, ci_low=lq, ci_high=hq),61])62table.to_csv(f"{OUT}/table1.csv", index=False)6364# ---- figure -------------------------------------------------------------------65fig, ax = plt.subplots(figsize=(7, 5))66ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.35, color="grey", label="firms (analysed)")67xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)68ax.plot(xs, m0.params["Intercept"] + b0 * xs, color="C0", lw=2, label=f"unadjusted: slope {b0:.3f}")69sk = d.skill_index.mean()70ax.plot(xs, m1.params["Intercept"] + b1 * xs + m1.params["skill_index"] * sk, color="C3", lw=2,71 label=f"adjusted for skill_index (at mean): slope {b1:.3f}")72ax.set_xlabel("training hours per employee (programme year)")73ax.set_ylabel("productivity growth, following year (%)")74ax.set_title("Productivity growth vs training hours")75ax.legend(fontsize=8)76fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)7778res = {79 "task_id": "F6-conf-04",80 "feasibility": "complete",81 "feasibility_reason": None,82 "data_issues": [83 {"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}), which is impossible",84 "handling": "excluded from main analysis (exclusions.csv step 3); sensitivity estimate including them reported (sensitivity_all_rows_*)"},85 {"issue": "training_hours strongly correlated with pre-programme skill_index (r=%.2f), and skill_index strongly negatively related to growth: confounding" % d[["skill_index","training_hours"]].corr().iloc[0,1],86 "handling": "adjusted for skill_index (selection-on-observables model)"},87 ],88 "n_input_rows": {"firms.csv": int(n_raw)},89 "n_analysed": int(len(d)),90 "random_seed": None,91 "naive_estimate": b0, "naive_se": se0, "naive_ci_low": l0, "naive_ci_high": h0,92 "baseline_metric": rmse0,93 "model_metric": rmse1,94 "effect_estimate": b1, "effect_se": se1, "effect_se_type": "HC1 heteroskedasticity-robust",95 "effect_ci_low": l1, "effect_ci_high": h1,96 "effect_se_conventional": float(mconv.bse["training_hours"]),97 "skill_index_coefficient": float(m1.params["skill_index"]),98 "skill_index_se": float(m1.bse["skill_index"]),99 "conclusion": concl,100 "corr_skill_hours": float(d[["skill_index", "training_hours"]].corr().iloc[0, 1]),101 "sensitivity_all_rows_estimate": bs, "sensitivity_all_rows_ci_low": ls, "sensitivity_all_rows_ci_high": hs,102 "sensitivity_all_rows_naive_estimate": bn_all,103 "sensitivity_quad_skill_estimate": bq, "sensitivity_quad_skill_ci_low": lq, "sensitivity_quad_skill_ci_high": hq,104 "diag_interaction_coef": float(mi.params["training_hours:skill_index"]),105 "diag_interaction_pvalue": float(mi.pvalues["training_hours:skill_index"]),106 "diag_hours_sq_coef": float(mh.params["I(training_hours ** 2)"]),107 "diag_hours_sq_pvalue": float(mh.pvalues["I(training_hours ** 2)"]),108 "r2_naive": float(m0.rsquared), "r2_adjusted": float(m1.rsquared),109}110with open(f"{OUT}/results.json", "w") as f:111 json.dump(res, f, indent=2)112print(json.dumps({k: res[k] for k in ["naive_estimate", "effect_estimate", "effect_ci_low", "effect_ci_high", "conclusion", "n_analysed"]}, indent=1))113