1"""Agent G — preregistered estimation from data/validated/runs_<tag>.parquet.2Outputs: analysis/tables/*.csv (+ .tex) and analysis/derived/estimates_<tag>.json. All numbers in the manuscript come from here.3Usage: .venv/bin/python analysis/econometrics/estimate.py --tag pilot [--q 70] [--B 5000] [--perm 20000]"""4from __future__ import annotations5import argparse, json, os, warnings6import numpy as np, pandas as pd7warnings.filterwarnings("ignore")8REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))9TAB = os.path.join(REPO, "analysis", "tables"); DER = os.path.join(REPO, "analysis", "derived")10ARMS = ["A0", "A1", "A2", "A3", "A4"]11CONTRASTS = [("A3", "A1", "primary"), ("A1", "A0", "secondary"), ("A4", "A3", "secondary"), ("A2", "A1", "secondary")]1213def load(tag, q):14 df = pd.read_parquet(os.path.join(REPO, "data", "validated", f"runs_{tag}.parquet"))15 df = df[df["status"].isin(["collected", "evaluated"])].copy() # metrics are computed before the run is marked evaluated16 df["S"] = ((df["det_completion"].astype(bool)) & (~df["critical_error"].astype(bool)) & (df["quality"].astype(float) >= q)).astype(float)17 df.loc[df["quality"].isna(), "S"] = np.nan18 df["D"] = df["det_completion"].astype(float); df["E"] = df["critical_error"].astype(float)19 df["C"] = df["total_cost_usd"].astype(float); df["T"] = df["wall_total_s"].astype(float)20 df["censored"] = df["timed_out"].astype(bool)21 # Killed sessions emit no `result` event, so their marginal cost is unobserved (cost_usd = 0). Two cases:22 # (a) an A4 run whose recovery session resumed the same conversation: the client reports a CUMULATIVE total_cost_usd for the23 # resumed session (WP29 gotcha), stored in recovery_cost_usd → that total already covers both sessions, no imputation;24 # (b) a killed session without resume (A3 interrupted control): impute from wall time with a per-model slope fitted by least25 # squares on the completed pilot runs (cost ≈ β·wall_s; R² = 0.93 for Opus; token counts streamed per event are not usable26 # because per-event output usage is partial — see data section). Flagged in cost_imputed_flag.27 df["cost_observed"] = df["cost_usd"].fillna(0).astype(float) + df["recovery_cost_usd"].fillna(0).astype(float); df["cost_imputed_flag"] = False; df["cost_cumulative_resume"] = False28 try:29 ref = pd.read_parquet(os.path.join(REPO, "data", "validated", "runs_pilot.parquet")); ref = ref[ref["cost_usd"].fillna(0) > 0]30 except Exception: ref = df[df["cost_usd"].fillna(0) > 0] # noqa: BLE00131 for mdl, g in df.groupby("model"):32 fit = ref[ref["model"] == mdl]33 resumed = g[(g["cost_usd"].fillna(0) == 0) & (g["recovery_cost_usd"].fillna(0) > 0)]34 df.loc[resumed.index, "C"] = resumed["recovery_cost_usd"].astype(float); df.loc[resumed.index, "cost_cumulative_resume"] = True35 miss = g[(g["cost_usd"].fillna(0) == 0) & (g["recovery_cost_usd"].fillna(0) == 0) & g["censored"]]36 if len(fit) >= 5 and len(miss):37 w = fit["wall_s"].astype(float).values; beta = float((w * fit["cost_usd"].values).sum() / (w * w).sum())38 df.loc[miss.index, "C"] = beta * miss["wall_s"].astype(float); df.loc[miss.index, "cost_imputed_flag"] = True39 return df4041def task_means(df, col, arm):42 return df[df.arm == arm].groupby("task_id")[col].mean()4344def weighted_rate(df, col, arm, weights=None):45 tm = task_means(df, col, arm).dropna()46 if weights is None: return float(tm.mean())47 w = weights.reindex(tm.index).fillna(0); return float((tm * w).sum() / w.sum())4849def cluster_bootstrap(df, col, a, b, B, rng, weights=None):50 ta = task_means(df, col, a); tb = task_means(df, col, b); tasks = sorted(set(ta.index) & set(tb.index))51 d = np.array([ta[t] - tb[t] for t in tasks]); w = np.array([weights[t] if weights is not None else 1.0 for t in tasks])52 est = float((d * w).sum() / w.sum()); boots = []53 for _ in range(B):54 idx = rng.integers(0, len(tasks), len(tasks)); boots.append((d[idx] * w[idx]).sum() / w[idx].sum())55 lo, hi = np.percentile(boots, [2.5, 97.5]); return est, float(lo), float(hi), len(tasks)5657def permutation_p(df, col, a, b, n_perm, rng):58 ta = task_means(df, col, a); tb = task_means(df, col, b); tasks = sorted(set(ta.index) & set(tb.index))59 d = np.array([ta[t] - tb[t] for t in tasks]); obs = d.mean(); cnt = 060 for _ in range(n_perm):61 s = rng.choice([-1, 1], len(d)); cnt += abs((d * s).mean()) >= abs(obs) - 1e-1262 return float(cnt / n_perm)6364def holm(pvals):65 idx = np.argsort(pvals); m = len(pvals); adj = np.empty(m); run = 066 for r, i in enumerate(idx): run = max(run, (m - r) * pvals[i]); adj[i] = min(1.0, run)67 return adj.tolist()6869def mixed_logit(df, col):70 try:71 import statsmodels.formula.api as smf72 d = df.dropna(subset=[col]).copy(); d["y"] = d[col].astype(int)73 md = smf.mixedlm if False else None74 import statsmodels.api as sm75 X = pd.get_dummies(d["arm"], drop_first=False).astype(float); X = X[[a for a in ARMS if a in X.columns]]76 m = sm.GLM(d["y"], X, family=sm.families.Binomial()).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["task_id"])[0]})77 out = {a: {"logit": float(m.params[a]), "se_cluster": float(m.bse[a]), "p": float(1 / (1 + np.exp(-m.params[a])))} for a in X.columns}78 return {"model": "cluster-robust logit (arm dummies, no intercept), clustered by scenario", "params": out}79 except Exception as e: # noqa: BLE00180 return {"error": repr(e)}8182def survival(df):83 try:84 from lifelines import KaplanMeierFitter85 out = {}86 for a in ARMS:87 d = df[df.arm == a]88 if len(d) == 0: continue89 km = KaplanMeierFitter(); km.fit(d["T"], event_observed=~d["censored"])90 out[a] = {"median_s": float(km.median_survival_time_), "n": int(len(d)), "censored": int(d["censored"].sum()),91 "p_done_600": float(1 - km.survival_function_at_times(600).iloc[0]), "p_done_1200": float(1 - km.survival_function_at_times(1200).iloc[0])}92 return out93 except Exception as e: # noqa: BLE00194 return {"error": repr(e)}9596def alf_curves(df, budgets, weights=None, h_allow=0.0, q_col="S"):97 """ALF_o(B,H,q) = Σ_j w_oj P(S=1, C≤B, H≤H). H = 5 min per unanswered question (assumption) + logged interventions."""98 rows = []99 for occ, g in df.groupby("occupation"):100 for a in ARMS:101 ga = g[g.arm == a]102 if len(ga) == 0: continue103 for B in budgets:104 ok = (ga[q_col] == 1) & (ga["C"] <= B) & (ga["H_min"] <= h_allow)105 tm = ok.groupby(ga["task_id"]).mean()106 w = (weights.reindex(tm.index).fillna(0) if weights is not None else pd.Series(1.0, index=tm.index))107 rows.append({"occupation": occ, "arm": a, "B": B, "H": h_allow, "ALF": float((tm * w).sum() / w.sum())})108 return pd.DataFrame(rows)109110def agreement(jdg):111 """Krippendorff's alpha (interval) across judges on the per-run quality scores + pairwise Pearson."""112 piv = jdg.pivot_table(index="run_id", columns="judge_idx", values="score")113 try:114 vals = piv.values.astype(float); n_units = vals.shape[0]115 obs = []; allv = vals[~np.isnan(vals)]116 De = np.var(allv) * 2 if len(allv) > 1 else np.nan117 Do = 0.0; m = 0118 for row in vals:119 r = row[~np.isnan(row)]120 if len(r) > 1:121 Do += np.sum((r[:, None] - r[None, :]) ** 2) / (len(r) - 1); m += len(r)122 Do = Do / m if m else np.nan123 alpha = 1 - Do / De if De and not np.isnan(De) else np.nan124 except Exception: alpha = np.nan # noqa: BLE001125 corr = piv.corr().round(3).to_dict() if piv.shape[1] > 1 else {}126 return {"krippendorff_alpha_interval": float(alpha) if alpha == alpha else None, "pairwise_pearson": corr, "n_runs": int(piv.shape[0])}127128def main():129 ap = argparse.ArgumentParser(); ap.add_argument("--tag", default="pilot"); ap.add_argument("--q", type=float, default=70); ap.add_argument("--B", type=int, default=5000); ap.add_argument("--perm", type=int, default=20000); ap.add_argument("--seed", type=int, default=42)130 a = ap.parse_args(); rng = np.random.default_rng(a.seed); os.makedirs(TAB, exist_ok=True); os.makedirs(DER, exist_ok=True)131 df = load(a.tag, a.q)132 # human-time proxy: 5 min per question in COMPLETION.json (assumption) + logged interventions133 def nq(x):134 try: return len(json.loads(x)) if isinstance(x, str) else 0135 except Exception: return 0 # noqa: BLE001136 df["H_min"] = df["human_intervention_min"].astype(float) + 5.0 * df.get("n_questions", pd.Series(0, index=df.index)).fillna(0).astype(float) if "n_questions" in df else df["human_intervention_min"].astype(float)137 # task weights: equal (primary) and professional minutes (sensitivity)138 tasks = {t: json.load(open(os.path.join(REPO, "benchmark", "scenarios", t, "task.json"))) for t in df.task_id.unique()}139 w_min = pd.Series({t: float(v.get("professional_minutes_estimate", 60)) for t, v in tasks.items()})140 # --- Table 1: outcomes by arm141 rows = []142 for arm in ARMS:143 d = df[df.arm == arm]144 if len(d) == 0: continue145 S = d["S"]; acc = S.sum(); rows.append({"arm": arm, "n_runs": len(d), "n_tasks": d.task_id.nunique(), "acceptance_taskw": weighted_rate(df, "S", arm), "acceptance_minw": weighted_rate(df, "S", arm, w_min),146 "det_completion": weighted_rate(df, "D", arm), "critical_error_rate": float(d["E"].mean()), "quality_mean": float(d["quality"].mean()), "quality_sd": float(d["quality"].std()),147 "timeout_rate": float(d["censored"].mean()), "cost_mean": float(d["C"].mean()), "cost_sd": float(d["C"].std()), "cost_per_accepted": float(d["C"].sum() / acc) if acc else np.nan,148 "wall_median_s": float(d["T"].median()), "tool_calls_mean": float(d["tool_calls"].mean()), "tool_errors_mean": float(d["tool_errors"].mean()), "app_transitions_mean": float(d["app_transitions"].mean()),149 "abstention_rate": float((d["abstention_status"].isin(["needs_clarification", "abstained"])).mean())})150 t1 = pd.DataFrame(rows); t1.to_csv(os.path.join(TAB, f"table1_outcomes_by_arm_{a.tag}.csv"), index=False)151 # --- Table 2: contrasts152 crows = []; pv = []153 for x, y, kind in CONTRASTS:154 for col in ("S", "D", "E", "quality"):155 if col == "quality":156 est, lo, hi, J = cluster_bootstrap(df, "quality", x, y, a.B, rng)157 p = permutation_p(df, "quality", x, y, a.perm, rng)158 else:159 est, lo, hi, J = cluster_bootstrap(df, col, x, y, a.B, rng); p = permutation_p(df, col, x, y, a.perm, rng)160 est_w = cluster_bootstrap(df, col, x, y, 200, rng, w_min)[0]161 crows.append({"contrast": f"{x}-{y}", "kind": kind, "outcome": col, "estimate": est, "ci_lo": lo, "ci_hi": hi, "perm_p": p, "J": J, "estimate_minutes_weighted": est_w})162 if col == "S" and kind == "secondary": pv.append(p)163 t2 = pd.DataFrame(crows)164 if pv:165 adj = holm(pv); k = 0166 t2["holm_p"] = np.nan167 for i, r in t2.iterrows():168 if r["outcome"] == "S" and r["kind"] == "secondary": t2.loc[i, "holm_p"] = adj[k]; k += 1169 t2.to_csv(os.path.join(TAB, f"table2_contrasts_{a.tag}.csv"), index=False)170 # --- Table 3: by occupation × arm171 t3 = df.groupby(["occupation", "arm"]).agg(n=("S", "size"), acceptance=("S", "mean"), det=("D", "mean"), quality=("quality", "mean"), critical=("E", "mean"), cost=("C", "mean"), tool_calls=("tool_calls", "mean")).reset_index()172 t3.to_csv(os.path.join(TAB, f"table3_occupation_arm_{a.tag}.csv"), index=False)173 # --- dependency: acceptance vs designed depth and vs tool calls (within arm)174 df["depth"] = df.task_id.map({t: v.get("designed_dependency_depth") for t, v in tasks.items()})175 dep = df.groupby(["arm", "depth"]).agg(n=("S", "size"), acceptance=("S", "mean"), tool_calls=("tool_calls", "mean")).reset_index(); dep.to_csv(os.path.join(TAB, f"table4_dependency_{a.tag}.csv"), index=False)176 # --- recovery index (A4)177 a4 = df[df.arm == "A4"]; fail_first = a4[a4["selfcheck_pass_first"] == 0]178 ari = {"n_A4": int(len(a4)), "n_failed_first_selfcheck": int(len(fail_first)), "n_recovery_attempted": int(a4["recovery_attempted"].astype(bool).sum()),179 "n_recovered": int((fail_first["selfcheck_pass_final"] == 1).sum()), "ARI": float((fail_first["selfcheck_pass_final"] == 1).mean()) if len(fail_first) else None,180 "acceptance_given_recovery": float(fail_first["S"].mean()) if len(fail_first) else None, "mean_recovery_cost": float(a4["recovery_cost_usd"].fillna(0).mean())}181 # --- ALF curves182 budgets = [0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4]183 alf0 = alf_curves(df, budgets, None, 0.0); alf0["weights"] = "equal"; alf1 = alf_curves(df, budgets, w_min, 0.0); alf1["weights"] = "minutes"184 alfH = alf_curves(df, budgets, None, 15.0); alfH["weights"] = "equal_H15"185 alf = pd.concat([alf0, alf1, alfH]); alf.to_csv(os.path.join(TAB, f"table5_alf_{a.tag}.csv"), index=False)186 # --- sensitivity to q187 qs = []188 for q in (50, 60, 70, 80, 90):189 d2 = df.copy(); d2["S"] = ((d2["D"] == 1) & (d2["E"] == 0) & (d2["quality"] >= q)).astype(float)190 for arm in ARMS:191 if (d2.arm == arm).any(): qs.append({"q": q, "arm": arm, "acceptance": weighted_rate(d2, "S", arm)})192 pd.DataFrame(qs).to_csv(os.path.join(TAB, f"table6_q_sensitivity_{a.tag}.csv"), index=False)193 # --- judges agreement194 jdg = pd.read_parquet(os.path.join(REPO, "data", "validated", f"judgments_{a.tag}.parquet"))195 agr = agreement(jdg)196 jm = jdg.groupby("judge_model")["score"].agg(["mean", "std", "count"]).reset_index(); jm.to_csv(os.path.join(TAB, f"table7_judges_{a.tag}.csv"), index=False)197 # det vs judge-acceptable agreement198 df["judge_majority_acceptable"] = df["judge_acceptable_votes"].apply(lambda s: (sum(1 for v in json.loads(s) if v) >= 2) if isinstance(s, str) else None)199 cross = pd.crosstab(df["D"], df["judge_majority_acceptable"]); cross.to_csv(os.path.join(TAB, f"table8_det_vs_judges_{a.tag}.csv"))200 out = {"tag": a.tag, "q": a.q, "n_runs": int(len(df)), "n_tasks": int(df.task_id.nunique()), "total_cost_usd": float(df["C"].sum()), "by_arm": t1.to_dict("records"), "contrasts": t2.to_dict("records"),201 "mixed_logit_S": mixed_logit(df, "S"), "survival": survival(df), "recovery": ari, "agreement": agr, "models": sorted(df["model"].dropna().unique().tolist()),202 "claude_code_versions": sorted(df["claude_code_version"].dropna().unique().tolist())}203 json.dump(out, open(os.path.join(DER, f"estimates_{a.tag}.json"), "w"), indent=1, default=str)204 print(t1.round(3).to_string()); print(t2.round(3).to_string()); print(json.dumps(ari, indent=1)); print(json.dumps(agr, indent=1, default=str))205206if __name__ == "__main__":207 main()208