1"""Replication entry point for F6-conf-04. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import matplotlib5matplotlib.use("Agg")6import matplotlib.pyplot as plt7import statsmodels.formula.api as smf89SEED = None # no randomness used10OUT = "out"11os.makedirs(OUT, exist_ok=True)12raw = pd.read_csv("data/firms.csv")13n_raw = len(raw)1415# Exclusion ledger16rows = []17d = raw.dropna()18rows.append(("1", "drop rows with missing values", n_raw, len(d), "no imputation allowed"))19n0 = len(d)20d = d[d["training_hours"] >= 0].copy()21rows.append(("2", "drop training_hours < 0", n0, len(d), "negative hours are physically impossible; treated as data errors"))22excl = pd.DataFrame(rows, columns=["step", "rule", "n_before", "n_after", "rationale"])23excl["n_removed"] = excl["n_before"] - excl["n_after"]24excl[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]].to_csv(f"{OUT}/exclusions.csv", index=False)25n_neg = int((raw["training_hours"] < 0).sum())2627COV = "HC1"28naive = smf.ols("productivity_growth ~ training_hours", d).fit(cov_type=COV)29adj = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit(cov_type=COV)30# sensitivity analyses31adj_full = smf.ols("productivity_growth ~ training_hours + skill_index", raw).fit(cov_type=COV)32adj_quad = smf.ols("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", d).fit(cov_type=COV)33adj_conv = smf.ols("productivity_growth ~ training_hours + skill_index", d).fit()3435def row(name, m):36 ci = m.conf_int().loc["training_hours"]37 return dict(model=name, estimate=float(m.params["training_hours"]), se=float(m.bse["training_hours"]),38 ci_low=float(ci[0]), ci_high=float(ci[1]))3940table = pd.DataFrame([41 row("unadjusted (baseline): growth ~ hours", naive),42 row("adjusted (main): growth ~ hours + skill_index", adj),43 row("sensitivity: adjusted, full sample incl. negative hours", adj_full),44 row("sensitivity: adjusted + skill_index^2", adj_quad),45 row("sensitivity: adjusted, conventional SE", adj_conv),46])47table.to_csv(f"{OUT}/table1.csv", index=False)4849rmse = lambda m: float(np.sqrt(np.mean(m.resid ** 2)))50lo, hi = float(adj.conf_int().loc["training_hours", 0]), float(adj.conf_int().loc["training_hours", 1])51concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")5253res = {54 "task_id": "F6-conf-04",55 "feasibility": "complete",56 "feasibility_reason": None,57 "data_issues": [{"issue": f"{n_neg} rows have negative training_hours (physically impossible)",58 "handling": "excluded from main analysis; full-sample estimate reported as sensitivity in table1.csv"}],59 "n_input_rows": {"firms.csv": int(n_raw)},60 "n_analysed": int(adj.nobs),61 "random_seed": SEED,62 "naive_estimate": float(naive.params["training_hours"]),63 "naive_ci_low": float(naive.conf_int().loc["training_hours", 0]),64 "naive_ci_high": float(naive.conf_int().loc["training_hours", 1]),65 "naive_se": float(naive.bse["training_hours"]),66 "baseline_metric": rmse(naive),67 "model_metric": rmse(adj),68 "effect_estimate": float(adj.params["training_hours"]),69 "effect_se": float(adj.bse["training_hours"]),70 "effect_ci_low": lo,71 "effect_ci_high": hi,72 "skill_index_coefficient": float(adj.params["skill_index"]),73 "se_type": "heteroskedasticity-robust HC1",74 "conclusion": concl,75 "sensitivity_full_sample_estimate": float(adj_full.params["training_hours"]),76 "sensitivity_quadratic_skill_estimate": float(adj_quad.params["training_hours"]),77 "sensitivity_conventional_se": float(adj_conv.bse["training_hours"]),78 "corr_skill_training": float(d["skill_index"].corr(d["training_hours"])),79}80with open(f"{OUT}/results.json", "w") as f:81 json.dump(res, f, indent=2)8283fig, ax = plt.subplots(figsize=(7, 5))84ax.scatter(d["training_hours"], d["productivity_growth"], s=8, alpha=0.4, color="grey", label="firms")85xs = np.linspace(d["training_hours"].min(), d["training_hours"].max(), 100)86ax.plot(xs, naive.params["Intercept"] + naive.params["training_hours"] * xs, "r-",87 label=f"unadjusted OLS (slope {res['naive_estimate']:.3f})")88sm = d["skill_index"].mean()89ax.plot(xs, adj.params["Intercept"] + adj.params["training_hours"] * xs + adj.params["skill_index"] * sm, "b-",90 label=f"adjusted for skill_index at mean (slope {res['effect_estimate']:.3f})")91ax.set_xlabel("training hours per employee"); ax.set_ylabel("productivity growth (%)")92ax.set_title("Productivity growth vs training hours (n=%d)" % res["n_analysed"]); ax.legend()93fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)94print(json.dumps({k: res[k] for k in ["naive_estimate", "effect_estimate", "effect_se", "conclusion", "n_analysed"]}, indent=1))95