1"""Replication entry point for F6-conf-04.2Run from workspace root: ~/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 recorded as null)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.formula.api as smf9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))14OUT = os.path.join(ROOT, "out")15raw = pd.read_csv(os.path.join(ROOT, "data", "firms.csv"))16n_raw = len(raw)1718# ---- exclusions ----19excl = []20d = raw.dropna(subset=["skill_index", "training_hours", "productivity_growth"])21excl.append(("1", "drop rows with missing values in analysis columns", n_raw, len(d),22 "estimation requires complete cases (none found)"))23n0 = len(d)24d = d[d.training_hours >= 0].copy()25excl.append(("2", "drop rows with training_hours < 0", n0, len(d),26 "negative training hours are physically impossible (invalid treatment value)"))27pd.DataFrame([{"step": s, "rule": r, "n_before": a, "n_after": b, "n_removed": a - b,28 "rationale": why} for s, r, a, b, why in excl]).to_csv(29 os.path.join(OUT, "exclusions.csv"), index=False)30n_neg = n0 - len(d)3132COV = "HC1" # heteroskedasticity-robust standard errors333435def 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]), rmse404142# ---- baseline (P2): unadjusted OLS ----43mn, b_n, se_n, lo_n, hi_n, rmse_n = fit("productivity_growth ~ training_hours", d)44# ---- target model (P3): adjust for pre-programme skill_index ----45ma, b_a, se_a, lo_a, hi_a, rmse_a = fit("productivity_growth ~ training_hours + skill_index", d)46# ---- sensitivity: full sample incl. negative hours; flexible skill specification ----47_, b_nf, se_nf, lo_nf, hi_nf, _ = fit("productivity_growth ~ training_hours", raw)48_, b_af, se_af, lo_af, hi_af, _ = fit("productivity_growth ~ training_hours + skill_index", raw)49_, b_q, se_q, lo_q, hi_q, _ = fit(50 "productivity_growth ~ training_hours + skill_index + I(skill_index**2) + I(skill_index**3)", d)51# conventional (non-robust) SE for the target model, for comparison52se_a_conv = float(smf.ols("productivity_growth ~ training_hours + skill_index", d).fit().bse["training_hours"])5354if lo_a > 0:55 concl = "positive_effect"56elif hi_a < 0:57 concl = "negative_effect"58else:59 concl = "no_detectable_effect"6061table = pd.DataFrame([62 ("unadjusted (main sample, n=%d)" % len(d), b_n, se_n, lo_n, hi_n),63 ("adjusted for skill_index (main sample, n=%d) [target]" % len(d), b_a, se_a, lo_a, hi_a),64 ("adjusted, cubic in skill_index (main sample, n=%d)" % len(d), b_q, se_q, lo_q, hi_q),65 ("unadjusted (all rows incl. negative hours, n=%d)" % n_raw, b_nf, se_nf, lo_nf, hi_nf),66 ("adjusted for skill_index (all rows incl. negative hours, n=%d)" % n_raw, b_af, se_af, lo_af, hi_af),67], columns=["model", "estimate", "se", "ci_low", "ci_high"])68table.to_csv(os.path.join(OUT, "table1.csv"), index=False)6970results = {71 "task_id": "F6-conf-04",72 "feasibility": "complete",73 "feasibility_reason": None,74 "data_issues": [75 {"issue": "%d rows have negative training_hours (min %.3f), which is physically impossible" % (n_neg, raw.training_hours.min()),76 "handling": "excluded from main analysis (exclusions.csv step 2); sensitivity with all %d rows reported (sensitivity_* keys, table1)" % n_raw},77 {"issue": "skill_index strongly correlated with both training_hours (r=%.2f) and productivity_growth (r=%.2f): confounding of the unadjusted slope" % (raw.skill_index.corr(raw.training_hours), raw.skill_index.corr(raw.productivity_growth)),78 "handling": "target model adjusts for skill_index (pre-programme covariate)"},79 ],80 "n_input_rows": {"firms.csv": int(n_raw)},81 "n_analysed": int(len(d)),82 "random_seed": None,83 "naive_estimate": b_n, "naive_se": se_n, "naive_ci_low": lo_n, "naive_ci_high": hi_n,84 "baseline_metric": rmse_n,85 "model_metric": rmse_a,86 "effect_estimate": b_a, "effect_se": se_a, "effect_ci_low": lo_a, "effect_ci_high": hi_a,87 "effect_se_type": "HC1 heteroskedasticity-robust",88 "effect_se_conventional": se_a_conv,89 "skill_index_coefficient": float(ma.params["skill_index"]),90 "skill_index_se": float(ma.bse["skill_index"]),91 "adjusted_r2": float(ma.rsquared), "naive_r2": float(mn.rsquared),92 "conclusion": concl,93 "sensitivity_cubic_skill_estimate": b_q, "sensitivity_cubic_skill_ci_low": lo_q, "sensitivity_cubic_skill_ci_high": hi_q,94 "sensitivity_fullsample_naive_estimate": b_nf,95 "sensitivity_fullsample_effect_estimate": b_af, "sensitivity_fullsample_effect_ci_low": lo_af, "sensitivity_fullsample_effect_ci_high": hi_af,96 "n_excluded_negative_hours": int(n_neg),97}98with open(os.path.join(OUT, "results.json"), "w") as f:99 json.dump(results, f, indent=2)100101# ---- figure ----102fig, ax = plt.subplots(figsize=(7, 5))103sc = ax.scatter(d.training_hours, d.productivity_growth, c=d.skill_index, cmap="viridis", s=8, alpha=0.6)104plt.colorbar(sc, ax=ax, label="skill_index (pre-programme)")105xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)106ax.plot(xs, mn.params["Intercept"] + b_n * xs, "r-", lw=2, label="unadjusted OLS: %.3f pp/h" % b_n)107sbar = d.skill_index.mean()108ax.plot(xs, ma.params["Intercept"] + b_a * xs + ma.params["skill_index"] * sbar, "k--", lw=2,109 label="adjusted for skill (at mean skill): %.3f pp/h" % b_a)110ax.set_xlabel("training hours per employee (programme year)")111ax.set_ylabel("productivity growth, following year (%)")112ax.set_title("Productivity growth vs training hours (n=%d)" % len(d))113ax.legend()114fig.tight_layout()115fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)116117print("n_analysed=%d naive=%.4f [%.4f, %.4f] effect=%.4f (HC1 se %.4f) [%.4f, %.4f] conclusion=%s"118 % (len(d), b_n, lo_n, hi_n, b_a, se_a, lo_a, hi_a, concl))119