1"""Replication entry point for F6-conf-04. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads data/firms.csv only; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No 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 plt1213OUT = "out"14os.makedirs(OUT, exist_ok=True)15df = pd.read_csv("data/firms.csv")16n_raw = len(df)1718# ---- audit ----19issues = []20n_neg = int((df.training_hours < 0).sum())21if n_neg:22 issues.append({"issue": f"{n_neg} rows have negative training_hours (min {df.training_hours.min():.3f}), physically impossible",23 "handling": "kept in main analysis (smooth continuous tail of simulated hours-skill relation; dropping would truncate on the treatment); sensitivity analysis excluding them reported in table1.csv"})24assert df.isna().sum().sum() == 0 and not df.firm_id.duplicated().any()2526# exclusion ledger: no rows excluded in main analysis27pd.DataFrame([{"step": 1, "rule": "drop rows with missing values or duplicate firm_id", "n_before": n_raw,28 "n_after": n_raw, "n_removed": 0, "rationale": "none found; all rows retained"}]29 ).to_csv(f"{OUT}/exclusions.csv", index=False)3031def fit(formula, data):32 return smf.ols(formula, data=data).fit(cov_type="HC1")3334def rmse(m):35 return float(np.sqrt(np.mean(m.resid ** 2)))3637naive = fit("productivity_growth ~ training_hours", df)38main = fit("productivity_growth ~ training_hours + skill_index", df)39sens_pos = fit("productivity_growth ~ training_hours + skill_index", df[df.training_hours >= 0])40sens_quad = fit("productivity_growth ~ training_hours + skill_index + I(skill_index**2)", df)41# conventional-SE version of main model for comparison42main_conv = smf.ols("productivity_growth ~ training_hours + skill_index", data=df).fit()4344def row(name, m):45 ci = m.conf_int().loc["training_hours"]46 return {"model": name, "estimate": float(m.params["training_hours"]), "se": float(m.bse["training_hours"]),47 "ci_low": float(ci[0]), "ci_high": float(ci[1]), "n": int(m.nobs), "rmse": rmse(m)}4849rows = [row("unadjusted: y ~ hours (HC1)", naive),50 row("adjusted (main): y ~ hours + skill (HC1)", main),51 row("adjusted, conventional SE", main_conv),52 row("sensitivity: main excl. negative hours (HC1)", sens_pos),53 row("sensitivity: + skill^2 (HC1)", sens_quad)]54pd.DataFrame(rows).to_csv(f"{OUT}/table1.csv", index=False)5556r_n, r_m = rows[0], rows[1]57lo, hi = r_m["ci_low"], r_m["ci_high"]58conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")5960results = {61 "task_id": "F6-conf-04",62 "feasibility": "complete",63 "feasibility_reason": None,64 "data_issues": issues,65 "n_input_rows": {"firms.csv": n_raw},66 "n_analysed": int(main.nobs),67 "random_seed": None,68 "naive_estimate": r_n["estimate"],69 "naive_ci_low": r_n["ci_low"],70 "naive_ci_high": r_n["ci_high"],71 "baseline_metric": r_n["rmse"],72 "model_metric": r_m["rmse"],73 "effect_estimate": r_m["estimate"],74 "effect_se": r_m["se"],75 "effect_se_type": "HC1 heteroskedasticity-robust",76 "effect_ci_low": lo,77 "effect_ci_high": hi,78 "skill_index_coefficient": float(main.params["skill_index"]),79 "skill_index_se": float(main.bse["skill_index"]),80 "corr_skill_hours": float(df.skill_index.corr(df.training_hours)),81 "sens_excl_negative_hours_estimate": rows[3]["estimate"],82 "sens_excl_negative_hours_n": rows[3]["n"],83 "sens_skill_quadratic_estimate": rows[4]["estimate"],84 "sens_skill_quadratic_sq_coef_pvalue": float(sens_quad.pvalues["I(skill_index ** 2)"]),85 "conclusion": conclusion,86}87with open(f"{OUT}/results.json", "w") as f:88 json.dump(results, f, indent=2)8990# ---- figure ----91fig, ax = plt.subplots(figsize=(7, 5))92ax.scatter(df.training_hours, df.productivity_growth, s=6, alpha=0.4, color="grey", label="firms")93x = np.linspace(df.training_hours.min(), df.training_hours.max(), 100)94ax.plot(x, naive.params["Intercept"] + naive.params["training_hours"] * x, color="C3",95 label=f"unadjusted (slope {r_n['estimate']:.3f})")96ax.plot(x, main.params["Intercept"] + main.params["training_hours"] * x + main.params["skill_index"] * df.skill_index.mean(),97 color="C0", label=f"adjusted for skill, at mean skill (slope {r_m['estimate']:.3f})")98ax.set_xlabel("training hours per employee")99ax.set_ylabel("productivity growth next year (%)")100ax.legend()101ax.set_title("Productivity growth vs training hours (simulated)")102fig.tight_layout()103fig.savefig(f"{OUT}/figure1.png", dpi=120)104print(json.dumps(results, indent=2))105print(pd.DataFrame(rows).to_string())106