"""Replication entry point for F6-conf-04. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py Reads only data/firms.csv; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv. No randomness is used (random_seed = null).""" import json, os import numpy as np import pandas as pd import statsmodels.formula.api as smf import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt OUT = "out" os.makedirs(OUT, exist_ok=True) raw = pd.read_csv("data/firms.csv") n_raw = len(raw) # ---- exclusions ------------------------------------------------------------- excl = [] d = raw.copy() def step(name, rule, mask, why): global d nb = len(d); d = d[mask(d)].copy(); na = len(d) excl.append(dict(step=name, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why)) step("1", "drop rows with missing values in any column", lambda x: x.notna().all(axis=1), "complete cases required (none expected)") step("2", "drop duplicated firm_id", lambda x: ~x.firm_id.duplicated(), "one row per firm") step("3", "drop training_hours < 0", lambda x: x.training_hours >= 0, "negative hours of training are logically impossible; treated as invalid treatment records") pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False) n_neg = int((raw.training_hours < 0).sum()) # ---- models ----------------------------------------------------------------- COV = "HC1" # heteroskedasticity-robust SEs def fit(formula, data): m = smf.ols(formula, data).fit(cov_type=COV) ci = m.conf_int(0.05).loc["training_hours"] rmse = float(np.sqrt(np.mean(m.resid ** 2))) return m, float(m.params["training_hours"]), float(m.bse["training_hours"]), float(ci[0]), float(ci[1]), rmse m0, b0, se0, l0, h0, rmse0 = fit("productivity_growth ~ training_hours", d) m1, b1, se1, l1, h1, rmse1 = fit("productivity_growth ~ training_hours + skill_index", d) # sensitivity: all raw rows (including negative hours) _, bs, ses, ls, hs, _ = fit("productivity_growth ~ training_hours + skill_index", raw) _, bn_all, sen_all, ln_all, hn_all, _ = fit("productivity_growth ~ training_hours", raw) # diagnostics (not the estimand): quadratic in skill / hours, interaction _, bq, seq, lq, hq, rmseq = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d) mi = smf.ols("productivity_growth ~ training_hours * skill_index", d).fit(cov_type=COV) mh = smf.ols("productivity_growth ~ training_hours + I(training_hours**2) + skill_index", d).fit(cov_type=COV) mconv = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit() if l1 > 0: concl = "positive_effect" elif h1 < 0: concl = "negative_effect" else: concl = "no_detectable_effect" table = pd.DataFrame([ dict(model="unadjusted OLS (growth ~ hours), HC1", estimate=b0, se=se0, ci_low=l0, ci_high=h0), dict(model="adjusted OLS (growth ~ hours + skill_index), HC1 [main]", estimate=b1, se=se1, ci_low=l1, ci_high=h1), dict(model="sensitivity: adjusted OLS incl. 12 negative-hours rows, HC1", estimate=bs, se=ses, ci_low=ls, ci_high=hs), dict(model="sensitivity: adjusted OLS + skill_index^2, HC1", estimate=bq, se=seq, ci_low=lq, ci_high=hq), ]) table.to_csv(f"{OUT}/table1.csv", index=False) # ---- figure ------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(7, 5)) ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.35, color="grey", label="firms (analysed)") xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100) ax.plot(xs, m0.params["Intercept"] + b0 * xs, color="C0", lw=2, label=f"unadjusted: slope {b0:.3f}") sk = d.skill_index.mean() ax.plot(xs, m1.params["Intercept"] + b1 * xs + m1.params["skill_index"] * sk, color="C3", lw=2, label=f"adjusted for skill_index (at mean): slope {b1:.3f}") ax.set_xlabel("training hours per employee (programme year)") ax.set_ylabel("productivity growth, following year (%)") ax.set_title("Productivity growth vs training hours") ax.legend(fontsize=8) fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig) res = { "task_id": "F6-conf-04", "feasibility": "complete", "feasibility_reason": None, "data_issues": [ {"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}), which is impossible", "handling": "excluded from main analysis (exclusions.csv step 3); sensitivity estimate including them reported (sensitivity_all_rows_*)"}, {"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], "handling": "adjusted for skill_index (selection-on-observables model)"}, ], "n_input_rows": {"firms.csv": int(n_raw)}, "n_analysed": int(len(d)), "random_seed": None, "naive_estimate": b0, "naive_se": se0, "naive_ci_low": l0, "naive_ci_high": h0, "baseline_metric": rmse0, "model_metric": rmse1, "effect_estimate": b1, "effect_se": se1, "effect_se_type": "HC1 heteroskedasticity-robust", "effect_ci_low": l1, "effect_ci_high": h1, "effect_se_conventional": float(mconv.bse["training_hours"]), "skill_index_coefficient": float(m1.params["skill_index"]), "skill_index_se": float(m1.bse["skill_index"]), "conclusion": concl, "corr_skill_hours": float(d[["skill_index", "training_hours"]].corr().iloc[0, 1]), "sensitivity_all_rows_estimate": bs, "sensitivity_all_rows_ci_low": ls, "sensitivity_all_rows_ci_high": hs, "sensitivity_all_rows_naive_estimate": bn_all, "sensitivity_quad_skill_estimate": bq, "sensitivity_quad_skill_ci_low": lq, "sensitivity_quad_skill_ci_high": hq, "diag_interaction_coef": float(mi.params["training_hours:skill_index"]), "diag_interaction_pvalue": float(mi.pvalues["training_hours:skill_index"]), "diag_hours_sq_coef": float(mh.params["I(training_hours ** 2)"]), "diag_hours_sq_pvalue": float(mh.pvalues["I(training_hours ** 2)"]), "r2_naive": float(m0.rsquared), "r2_adjusted": float(m1.rsquared), } with open(f"{OUT}/results.json", "w") as f: json.dump(res, f, indent=2) print(json.dumps({k: res[k] for k in ["naive_estimate", "effect_estimate", "effect_ci_low", "effect_ci_high", "conclusion", "n_analysed"]}, indent=1))