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.4Deterministic: no 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 plt1213DATA = "data/firms.csv"14OUT = "out"15os.makedirs(OUT, exist_ok=True)1617raw = pd.read_csv(DATA)18n_raw = len(raw)19cols = ["firm_id", "skill_index", "training_hours", "productivity_growth"]20assert list(raw.columns) == cols2122excl = []23data_issues = []24d = raw.copy()2526# Step 1: missing values (none expected)27nb = len(d); d = d.dropna(subset=cols); excl.append(("1", "drop rows with missing values", nb, len(d), "required variables must be observed"))28# Step 2: duplicate firm_id (none expected)29nb = len(d); d = d.drop_duplicates(subset="firm_id"); excl.append(("2", "drop duplicate firm_id", nb, len(d), "one row per firm"))30# Step 3: negative training hours (physically impossible)31n_neg = int((d.training_hours < 0).sum())32nb = len(d); d = d[d.training_hours >= 0]; excl.append(("3", "drop training_hours < 0", nb, len(d),33 "hours of training per employee cannot be negative; values are invalid (likely data error). Sensitivity with them retained is reported."))34if n_neg:35 data_issues.append({"issue": f"{n_neg} rows have negative training_hours (min {raw.training_hours.min():.3f}), which is physically impossible",36 "handling": "excluded from the main analysis; sensitivity analysis re-estimates both models on all 1200 rows (keys sens_all_rows_*)"})37data_issues.append({"issue": "training_hours is strongly correlated with pre-programme skill_index (r = %.3f in all 1200 rows), and skill_index is strongly negatively related to productivity_growth: the unadjusted slope is confounded" % raw[["skill_index","training_hours"]].corr().iloc[0,1],38 "handling": "skill_index (measured before the programme) is included as a control in the model that targets the estimand"})3940pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(41 n_removed=lambda x: x.n_before - x.n_after)[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]].to_csv(42 os.path.join(OUT, "exclusions.csv"), index=False)4344COV = "HC1" # heteroskedasticity-robust standard errors45def fit(formula, df):46 return smf.ols(formula, df).fit(cov_type=COV)4748def rmse(m):49 return float(np.sqrt(np.mean(m.resid ** 2)))5051naive = fit("productivity_growth ~ training_hours", d)52adj = fit("productivity_growth ~ training_hours + skill_index", d)53# robustness checks (not the main estimate)54quad = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d)55inter = fit("productivity_growth ~ training_hours * skill_index", d)56naive_all = fit("productivity_growth ~ training_hours", raw)57adj_all = fit("productivity_growth ~ training_hours + skill_index", raw)5859def row(name, m, var="training_hours"):60 ci = m.conf_int().loc[var]61 return {"model": name, "estimate": float(m.params[var]), "se": float(m.bse[var]),62 "ci_low": float(ci[0]), "ci_high": float(ci[1])}6364table = pd.DataFrame([65 row("M1 unadjusted: growth ~ hours", naive),66 row("M2 adjusted (main): growth ~ hours + skill_index", adj),67 row("R1 adds skill_index^2", quad),68 row("R2 adds hours x skill_index (effect at skill_index=0)", inter),69 row("S1 unadjusted, all 1200 rows incl. negative hours", naive_all),70 row("S2 adjusted, all 1200 rows incl. negative hours", adj_all),71])72table.to_csv(os.path.join(OUT, "table1.csv"), index=False)7374nci = naive.conf_int().loc["training_hours"]; aci = adj.conf_int().loc["training_hours"]75lo, hi = float(aci[0]), float(aci[1])76conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")7778res = {79 "task_id": "F6-conf-04",80 "feasibility": "complete",81 "feasibility_reason": None,82 "data_issues": data_issues,83 "n_input_rows": {"firms.csv": int(n_raw)},84 "n_analysed": int(adj.nobs),85 "random_seed": None,86 "naive_estimate": float(naive.params["training_hours"]),87 "naive_se": float(naive.bse["training_hours"]),88 "naive_ci_low": float(nci[0]), "naive_ci_high": float(nci[1]),89 "baseline_metric": rmse(naive),90 "model_metric": rmse(adj),91 "effect_estimate": float(adj.params["training_hours"]),92 "effect_se": float(adj.bse["training_hours"]),93 "effect_ci_low": lo, "effect_ci_high": hi,94 "skill_index_coefficient": float(adj.params["skill_index"]),95 "skill_index_se": float(adj.bse["skill_index"]),96 "se_type": "HC1 heteroskedasticity-robust; 95% CI = estimate +/- 1.96*SE (normal)",97 "conclusion": conclusion,98 "corr_skill_training": float(d[["skill_index", "training_hours"]].corr().iloc[0, 1]),99 "robust_quadratic_skill_effect": float(quad.params["training_hours"]),100 "robust_quadratic_skill_sq_pvalue": float(quad.pvalues["I(skill_index ** 2)"]),101 "robust_interaction_coef": float(inter.params["training_hours:skill_index"]),102 "robust_interaction_pvalue": float(inter.pvalues["training_hours:skill_index"]),103 "sens_all_rows_naive_estimate": float(naive_all.params["training_hours"]),104 "sens_all_rows_effect_estimate": float(adj_all.params["training_hours"]),105 "sens_all_rows_effect_se": float(adj_all.bse["training_hours"]),106 "n_excluded_negative_hours": n_neg,107}108with open(os.path.join(OUT, "results.json"), "w") as f:109 json.dump(res, f, indent=2)110111# Figure112fig, ax = plt.subplots(figsize=(7, 5))113ax.scatter(d.training_hours, d.productivity_growth, s=8, alpha=0.35, color="grey", label="firms (n=%d)" % len(d))114xs = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)115ax.plot(xs, naive.params["Intercept"] + naive.params["training_hours"] * xs, color="tab:red",116 label="unadjusted OLS (slope %.3f)" % naive.params["training_hours"])117sbar = d.skill_index.mean()118ax.plot(xs, adj.params["Intercept"] + adj.params["skill_index"] * sbar + adj.params["training_hours"] * xs,119 color="tab:blue", label="adjusted for skill_index, at mean skill (slope %.3f)" % adj.params["training_hours"])120ax.set_xlabel("training hours per employee (programme year)")121ax.set_ylabel("productivity growth, following year (%)")122ax.set_title("Training hours vs productivity growth (simulated data)")123ax.legend(fontsize=8)124fig.tight_layout()125fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)126127print(f"n_raw={n_raw} n_analysed={res['n_analysed']}")128print(f"naive slope={res['naive_estimate']:.4f} [{res['naive_ci_low']:.4f}, {res['naive_ci_high']:.4f}]")129print(f"adjusted slope={res['effect_estimate']:.4f} (HC1 se {res['effect_se']:.4f}) [{lo:.4f}, {hi:.4f}] -> {conclusion}")130