1"""Replication script for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np, pandas as pd6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt9import statsmodels.formula.api as smf1011os.makedirs("out", exist_ok=True)12raw = pd.read_csv("data/firms.csv")13n_raw = len(raw)14issues = []1516# ---- audit ----17assert raw.isna().sum().sum() == 018assert not raw.firm_id.duplicated().any()19neg = raw.training_hours < 020issues.append({21 "issue": f"{int(neg.sum())} rows have negative training_hours (min {raw.training_hours.min():.3f}), which is physically impossible (hours per employee cannot be < 0); likely an untruncated simulation draw or recording error.",22 "handling": "Excluded from the main analysis; sensitivity analysis re-estimates with them included (results reported in table1.csv and report)."})23d = raw[~neg].copy()2425# ---- models (HC1 robust SEs) ----26naive = smf.ols("productivity_growth ~ training_hours", d).fit(cov_type="HC1")27adj = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit(cov_type="HC1")28adj_all = smf.ols("productivity_growth ~ training_hours + skill_index", raw).fit(cov_type="HC1")29naive_all = smf.ols("productivity_growth ~ training_hours", raw).fit(cov_type="HC1")30# diagnostic: nonlinearity / heterogeneity checks31quad = smf.ols("productivity_growth ~ training_hours + I(training_hours**2) + skill_index + I(skill_index**2)", d).fit(cov_type="HC1")32inter = smf.ols("productivity_growth ~ training_hours * skill_index", d).fit(cov_type="HC1")3334rmse = lambda m: float(np.sqrt(np.mean(m.resid ** 2)))35def row(name, m):36 ci = m.conf_int().loc["training_hours"]37 return {"model": name, "estimate": m.params["training_hours"], "se": m.bse["training_hours"],38 "ci_low": ci[0], "ci_high": ci[1]}39table = pd.DataFrame([40 row("unadjusted (main sample)", naive),41 row("adjusted for skill_index (main sample) [targets estimand]", adj),42 row("unadjusted (all rows incl. negative hours)", naive_all),43 row("adjusted for skill_index (all rows incl. negative hours)", adj_all),44 row("adjusted + quadratics (diagnostic)", quad),45])46table.to_csv("out/table1.csv", index=False)4748eci = adj.conf_int().loc["training_hours"]; nci = naive.conf_int().loc["training_hours"]49if eci[0] > 0: concl = "positive_effect"50elif eci[1] < 0: concl = "negative_effect"51else: concl = "no_detectable_effect"5253res = {54 "task_id": "F6-conf-04",55 "naive_estimate": float(naive.params["training_hours"]),56 "naive_ci_low": float(nci[0]), "naive_ci_high": float(nci[1]),57 "baseline_metric": rmse(naive),58 "model_metric": rmse(adj),59 "effect_estimate": float(adj.params["training_hours"]),60 "effect_se": float(adj.bse["training_hours"]),61 "effect_ci_low": float(eci[0]), "effect_ci_high": float(eci[1]),62 "skill_index_coefficient": float(adj.params["skill_index"]),63 "conclusion": concl,64 "feasibility": "complete", "feasibility_reason": None,65 "data_issues": issues,66 "n_input_rows": {"firms.csv": int(n_raw)},67 "n_analysed": int(adj.nobs),68 "random_seed": None,69 "diagnostics": {70 "se_type": "HC1 heteroskedasticity-robust",71 "corr_training_skill": float(d[["training_hours", "skill_index"]].corr().iloc[0, 1]),72 "sensitivity_all_rows_effect": float(adj_all.params["training_hours"]),73 "sensitivity_all_rows_se": float(adj_all.bse["training_hours"]),74 "quad_training_sq_coef": float(quad.params["I(training_hours ** 2)"]),75 "quad_training_sq_p": float(quad.pvalues["I(training_hours ** 2)"]),76 "quad_skill_sq_p": float(quad.pvalues["I(skill_index ** 2)"]),77 "interaction_coef": float(inter.params["training_hours:skill_index"]),78 "interaction_p": float(inter.pvalues["training_hours:skill_index"]),79 "conventional_se_effect": float(smf.ols("productivity_growth ~ training_hours + skill_index", d).fit().bse["training_hours"]),80 "r2_naive": float(naive.rsquared), "r2_adjusted": float(adj.rsquared),81 },82}83with open("out/results.json", "w") as f:84 json.dump(res, f, indent=2)8586# ---- figure ----87fig, ax = plt.subplots(figsize=(7, 5))88ax.scatter(raw.training_hours[~neg], raw.productivity_growth[~neg], s=8, alpha=0.4, label="firms (analysed)")89ax.scatter(raw.training_hours[neg], raw.productivity_growth[neg], s=14, marker="x", c="k", label="excluded (negative hours)")90x = np.linspace(d.training_hours.min(), d.training_hours.max(), 100)91ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, "r-",92 label=f"unadjusted: slope {naive.params['training_hours']:.3f}")93sm = d.skill_index.mean()94ax.plot(x, adj.params["Intercept"] + adj.params["training_hours"] * x + adj.params["skill_index"] * sm, "g-",95 label=f"adjusted for skill (at mean skill): slope {adj.params['training_hours']:.3f}")96ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth (%)")97ax.set_title("Productivity growth vs training hours (simulated data)")98ax.legend(fontsize=8)99fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)100print(json.dumps(res, indent=2))101