feat: statistical analysis on the whole data

1 files changed, 776 insertions, 0 deletions

diff --git a/src/statistical_analysis.py b/src/statistical_analysis.py
new file mode 100644
index 0000000..9c6d080
--- /dev/null
+++ b/src/statistical_analysis.py
@@ -0,0 +1,776 @@
+# /// script
+# requires-python = ">=3.10"
+# dependencies = [
+#     "pandas",
+#     "numpy",
+#     "scipy",
+#     "statsmodels",
+# ]
+# ///
+"""
+Statistical Analysis of LLM Deception Experiment Results
+=========================================================
+Produces bulletproof statistics for the dissertation including:
+- Chi-squared tests of independence and homogeneity
+- Fisher's exact tests (where cell counts are small)
+- Cochran-Armitage trend tests (oversight as ordinal)
+- Cramér's V effect sizes
+- Odds ratios with 95% confidence intervals
+- Bonferroni-corrected pairwise comparisons
+- Breslow-Day test for homogeneity of odds ratios
+- Logistic regression for interaction effects
+
+Usage:
+    python src/statistical_analysis.py logs/full_experiment/results.csv
+"""
+
+import sys
+import csv
+import warnings
+from collections import Counter, defaultdict
+from itertools import combinations
+
+import numpy as np
+import pandas as pd
+from scipy import stats
+from scipy.stats import chi2_contingency, fisher_exact, chi2
+from statsmodels.stats.contingency_tables import Table2x2
+from statsmodels.stats.proportion import proportion_confint
+import statsmodels.api as sm
+
+warnings.filterwarnings("ignore", category=FutureWarning)
+
+# ---------------------------------------------------------------------------
+# Constants
+# ---------------------------------------------------------------------------
+
+DECEPTIVE_CATEGORIES = {"MISLEADING", "FABRICATED", "OMITTED"}
+HONEST_CATEGORIES = {"ACCURATE"}
+OVERSIGHT_ORDER = {"low": 0, "mid": 1, "high": 2}
+SIGNIFICANCE_LEVEL = 0.05
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def load_data(csv_path: str) -> pd.DataFrame:
+    """Load results CSV and add derived columns."""
+    df = pd.read_csv(csv_path)
+    df["is_deceptive_output"] = df["blackbox_category"].isin(DECEPTIVE_CATEGORIES).astype(int)
+    df["is_compliant_reasoning"] = (df["glassbox_category"] == "COMPLIANT_DECEPTION").astype(int)
+    df["is_strategic"] = (df["glassbox_sophistication"] == "STRATEGIC").astype(int)
+    df["oversight_ordinal"] = df["oversight"].map(OVERSIGHT_ORDER)
+    return df
+
+
+def cramers_v(contingency_table: np.ndarray) -> float:
+    """Compute Cramér's V from a contingency table."""
+    chi2_val = chi2_contingency(contingency_table)[0]
+    n = contingency_table.sum()
+    min_dim = min(contingency_table.shape) - 1
+    if min_dim == 0 or n == 0:
+        return 0.0
+    return np.sqrt(chi2_val / (n * min_dim))
+
+
+def cramers_v_bias_corrected(contingency_table: np.ndarray) -> float:
+    """Compute bias-corrected Cramér's V (Bergsma 2013)."""
+    chi2_val = chi2_contingency(contingency_table)[0]
+    n = contingency_table.sum()
+    r, k = contingency_table.shape
+    phi2 = chi2_val / n
+    phi2_corrected = max(0, phi2 - ((k - 1) * (r - 1)) / (n - 1))
+    k_corrected = k - ((k - 1) ** 2) / (n - 1)
+    r_corrected = r - ((r - 1) ** 2) / (n - 1)
+    min_corrected = min(k_corrected - 1, r_corrected - 1)
+    if min_corrected <= 0:
+        return 0.0
+    return np.sqrt(phi2_corrected / min_corrected)
+
+
+def odds_ratio_ci(table_2x2: np.ndarray, alpha: float = 0.05):
+    """
+    Compute odds ratio with confidence interval from a 2x2 table.
+    Table format: [[a, b], [c, d]]
+    OR = (a*d) / (b*c)
+    """
+    a, b, c, d = table_2x2[0, 0], table_2x2[0, 1], table_2x2[1, 0], table_2x2[1, 1]
+    # Add 0.5 continuity correction if any cell is 0
+    if 0 in (a, b, c, d):
+        a, b, c, d = a + 0.5, b + 0.5, c + 0.5, d + 0.5
+    odds_ratio = (a * d) / (b * c)
+    log_or = np.log(odds_ratio)
+    se_log_or = np.sqrt(1/a + 1/b + 1/c + 1/d)
+    z = stats.norm.ppf(1 - alpha / 2)
+    ci_lower = np.exp(log_or - z * se_log_or)
+    ci_upper = np.exp(log_or + z * se_log_or)
+    return odds_ratio, ci_lower, ci_upper
+
+
+def cochran_armitage_trend_test(counts: list, nobs: list) -> tuple:
+    """
+    Cochran-Armitage test for trend in proportions.
+    counts: list of successes at each ordinal level
+    nobs: list of total observations at each ordinal level
+    Uses equally-spaced scores [0, 1, 2, ...].
+    Returns (z_statistic, p_value).
+    """
+    k = len(counts)
+    scores = np.arange(k, dtype=float)
+    n = np.array(nobs, dtype=float)
+    x = np.array(counts, dtype=float)
+    N = n.sum()
+    p_hat = x.sum() / N
+    
+    t_bar = np.sum(scores * n) / N
+    numerator = np.sum(scores * x) - x.sum() * t_bar
+    denominator_sq = p_hat * (1 - p_hat) * (np.sum(scores**2 * n) - N * t_bar**2)
+    
+    if denominator_sq <= 0:
+        return 0.0, 1.0
+    
+    z = numerator / np.sqrt(denominator_sq)
+    p_value = 2 * stats.norm.sf(abs(z))  # two-sided
+    return z, p_value
+
+
+def print_section(title: str):
+    """Print a formatted section header."""
+    print(f"\n{'='*80}")
+    print(f"  {title}")
+    print(f"{'='*80}\n")
+
+
+def print_subsection(title: str):
+    """Print a formatted subsection header."""
+    print(f"\n  --- {title} ---\n")
+
+
+def interpret_effect_size(v: float) -> str:
+    """Interpret Cramér's V effect size (Cohen's conventions)."""
+    if v < 0.1:
+        return "negligible"
+    elif v < 0.3:
+        return "small"
+    elif v < 0.5:
+        return "medium"
+    else:
+        return "large"
+
+
+def sig_stars(p: float) -> str:
+    """Return significance stars."""
+    if p < 0.001:
+        return "***"
+    elif p < 0.01:
+        return "**"
+    elif p < 0.05:
+        return "*"
+    else:
+        return "n.s."
+
+
+def format_p(p: float) -> str:
+    """Format p-value for display."""
+    if p < 0.001:
+        return f"p < .001"
+    else:
+        return f"p = {p:.4f}"
+
+
+# ---------------------------------------------------------------------------
+# Analysis Functions
+# ---------------------------------------------------------------------------
+
+def analysis_1_scenario_effect(df: pd.DataFrame):
+    """Test H1: Deception rates differ between scenarios."""
+    print_section("1. SCENARIO EFFECT ON DECEPTION")
+    
+    # --- 1a. Output-level (blackbox) deception ---
+    print_subsection("1a. Output-level deception rate by scenario (blackbox)")
+    
+    ct = pd.crosstab(df["scenario"], df["is_deceptive_output"])
+    print("  Contingency table:")
+    for scenario in ct.index:
+        n_honest = ct.loc[scenario, 0] if 0 in ct.columns else 0
+        n_deceptive = ct.loc[scenario, 1] if 1 in ct.columns else 0
+        total = n_honest + n_deceptive
+        rate = n_deceptive / total * 100
+        ci_low, ci_high = proportion_confint(n_deceptive, total, alpha=0.05, method="wilson")
+        print(f"    {scenario}: {n_deceptive}/{total} ({rate:.1f}%) "
+              f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+    
+    chi2_val, p_val, dof, expected = chi2_contingency(ct.values)
+    v = cramers_v(ct.values)
+    v_corrected = cramers_v_bias_corrected(ct.values)
+    
+    print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+    print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+    print(f"  Cramér's V (bias-corrected) = {v_corrected:.3f}")
+    
+    # Odds ratio
+    table_2x2 = ct.values
+    or_val, or_ci_low, or_ci_high = odds_ratio_ci(table_2x2)
+    print(f"  Odds ratio = {or_val:.3f}  95% CI: [{or_ci_low:.3f}, {or_ci_high:.3f}]")
+    
+    # Fisher's exact for robustness
+    if ct.shape == (2, 2):
+        _, fisher_p = fisher_exact(ct.values)
+        print(f"  Fisher's exact test: {format_p(fisher_p)} {sig_stars(fisher_p)}")
+    
+    # --- 1b. Reasoning-level (glassbox) compliance ---
+    print_subsection("1b. Reasoning-level compliance by scenario (glassbox: COMPLIANT_DECEPTION)")
+    
+    ct_gb = pd.crosstab(df["scenario"], df["is_compliant_reasoning"])
+    print("  Contingency table:")
+    for scenario in ct_gb.index:
+        n_not = ct_gb.loc[scenario, 0] if 0 in ct_gb.columns else 0
+        n_compliant = ct_gb.loc[scenario, 1] if 1 in ct_gb.columns else 0
+        total = n_not + n_compliant
+        rate = n_compliant / total * 100
+        ci_low, ci_high = proportion_confint(n_compliant, total, alpha=0.05, method="wilson")
+        print(f"    {scenario}: {n_compliant}/{total} ({rate:.1f}%) "
+              f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+    
+    chi2_val, p_val, dof, expected = chi2_contingency(ct_gb.values)
+    v = cramers_v(ct_gb.values)
+    print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+    print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+    
+    or_val, or_ci_low, or_ci_high = odds_ratio_ci(ct_gb.values)
+    print(f"  Odds ratio = {or_val:.3f}  95% CI: [{or_ci_low:.3f}, {or_ci_high:.3f}]")
+
+
+def analysis_2_oversight_effect(df: pd.DataFrame):
+    """Test H2: Oversight level affects deception rates."""
+    print_section("2. OVERSIGHT EFFECT ON DECEPTION")
+    
+    for scenario in sorted(df["scenario"].unique()):
+        print_subsection(f"2a. Output-level deception by oversight — {scenario}")
+        
+        sdf = df[df["scenario"] == scenario]
+        ct = pd.crosstab(sdf["oversight"], sdf["is_deceptive_output"])
+        # Reorder rows
+        oversight_labels = [o for o in ["low", "mid", "high"] if o in ct.index]
+        ct = ct.loc[oversight_labels]
+        
+        print("  Contingency table:")
+        counts_for_trend = []
+        nobs_for_trend = []
+        for oversight in oversight_labels:
+            n_honest = ct.loc[oversight, 0] if 0 in ct.columns else 0
+            n_deceptive = ct.loc[oversight, 1] if 1 in ct.columns else 0
+            total = n_honest + n_deceptive
+            rate = n_deceptive / total * 100
+            ci_low, ci_high = proportion_confint(n_deceptive, total, alpha=0.05, method="wilson")
+            print(f"    {oversight}: {n_deceptive}/{total} ({rate:.1f}%) "
+                  f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+            counts_for_trend.append(n_deceptive)
+            nobs_for_trend.append(total)
+        
+        # Chi-squared test (3x2)
+        chi2_val, p_val, dof, expected = chi2_contingency(ct.values)
+        v = cramers_v(ct.values)
+        print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+        print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+        
+        # Check expected cell counts
+        min_expected = expected.min()
+        print(f"  Minimum expected cell count = {min_expected:.1f}", end="")
+        if min_expected < 5:
+            print("  ⚠ WARNING: Expected count < 5, chi-squared may be unreliable")
+        else:
+            print("  ✓ All expected counts ≥ 5")
+        
+        # Cochran-Armitage trend test (treats oversight as ordinal)
+        z_trend, p_trend = cochran_armitage_trend_test(counts_for_trend, nobs_for_trend)
+        print(f"\n  Cochran-Armitage trend test: z = {z_trend:.3f}, {format_p(p_trend)} {sig_stars(p_trend)}")
+        if p_trend < 0.05:
+            direction = "increasing" if z_trend > 0 else "decreasing"
+            print(f"    → Significant linear trend: deception {direction} with oversight level")
+        else:
+            print(f"    → No significant linear trend detected")
+        
+        # Pairwise comparisons (Bonferroni-corrected)
+        print_subsection(f"2b. Pairwise oversight comparisons — {scenario} (Bonferroni-corrected, α/3 = {SIGNIFICANCE_LEVEL/3:.4f})")
+        
+        pairs = list(combinations(oversight_labels, 2))
+        n_comparisons = len(pairs)
+        
+        for o1, o2 in pairs:
+            subset = sdf[sdf["oversight"].isin([o1, o2])]
+            ct_pair = pd.crosstab(subset["oversight"], subset["is_deceptive_output"])
+            ct_pair = ct_pair.loc[[o1, o2]]
+            
+            if ct_pair.shape == (2, 2):
+                chi2_p, p_p, dof_p, _ = chi2_contingency(ct_pair.values)
+                _, fisher_p = fisher_exact(ct_pair.values)
+                or_val, or_ci_low, or_ci_high = odds_ratio_ci(ct_pair.values)
+                
+                adjusted_alpha = SIGNIFICANCE_LEVEL / n_comparisons
+                sig = "SIG" if p_p < adjusted_alpha else "n.s."
+                
+                print(f"    {o1} vs {o2}:")
+                print(f"      χ²(1) = {chi2_p:.3f}, raw {format_p(p_p)}, "
+                      f"Fisher {format_p(fisher_p)}")
+                print(f"      OR = {or_val:.3f}  95% CI: [{or_ci_low:.3f}, {or_ci_high:.3f}]  "
+                      f"[Bonferroni: {sig}]")
+
+
+def analysis_3_model_effects(df: pd.DataFrame):
+    """Test H3: Models differ in deception rates."""
+    print_section("3. MODEL-LEVEL DIFFERENCES")
+    
+    for scenario in sorted(df["scenario"].unique()):
+        print_subsection(f"3a. Output deception by model — {scenario}")
+        
+        sdf = df[df["scenario"] == scenario]
+        ct = pd.crosstab(sdf["model"], sdf["is_deceptive_output"])
+        
+        print("  Model deception rates:")
+        model_rates = []
+        for model in sorted(ct.index):
+            n_honest = ct.loc[model, 0] if 0 in ct.columns else 0
+            n_deceptive = ct.loc[model, 1] if 1 in ct.columns else 0
+            total = n_honest + n_deceptive
+            rate = n_deceptive / total * 100
+            ci_low, ci_high = proportion_confint(n_deceptive, total, alpha=0.05, method="wilson")
+            print(f"    {model}: {n_deceptive}/{total} ({rate:.1f}%) "
+                  f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+            model_rates.append((model, rate))
+        
+        chi2_val, p_val, dof, expected = chi2_contingency(ct.values)
+        v = cramers_v(ct.values)
+        print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+        print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+        
+        min_expected = expected.min()
+        if min_expected < 5:
+            print(f"  ⚠ Minimum expected count = {min_expected:.1f} < 5")
+
+
+def analysis_4_oversight_x_scenario_interaction(df: pd.DataFrame):
+    """Test H4: Oversight effect differs between scenarios (interaction)."""
+    print_section("4. OVERSIGHT × SCENARIO INTERACTION")
+    
+    print_subsection("4a. Logistic regression: deception ~ scenario + oversight + scenario×oversight")
+    
+    # Encode variables
+    model_df = df.copy()
+    model_df["scenario_code"] = (model_df["scenario"] == "child_protection").astype(int)
+    model_df["oversight_num"] = model_df["oversight_ordinal"]
+    model_df["interaction"] = model_df["scenario_code"] * model_df["oversight_num"]
+    
+    X = model_df[["scenario_code", "oversight_num", "interaction"]]
+    X = sm.add_constant(X)
+    y = model_df["is_deceptive_output"]
+    
+    try:
+        logit_model = sm.Logit(y, X).fit(disp=0)
+        print(logit_model.summary2().tables[1].to_string())
+        
+        print(f"\n  Model fit: Pseudo R² = {logit_model.prsquared:.4f}, "
+              f"AIC = {logit_model.aic:.1f}, BIC = {logit_model.bic:.1f}")
+        
+        # Likelihood ratio test for interaction term
+        X_no_interaction = model_df[["scenario_code", "oversight_num"]]
+        X_no_interaction = sm.add_constant(X_no_interaction)
+        logit_reduced = sm.Logit(y, X_no_interaction).fit(disp=0)
+        
+        lr_stat = -2 * (logit_reduced.llf - logit_model.llf)
+        lr_p = chi2.sf(lr_stat, df=1)
+        print(f"\n  Likelihood ratio test for interaction: "
+              f"LR = {lr_stat:.3f}, {format_p(lr_p)} {sig_stars(lr_p)}")
+        if lr_p < 0.05:
+            print("    → Significant interaction: oversight affects scenarios differently")
+        else:
+            print("    → No significant interaction")
+        
+    except Exception as e:
+        print(f"  Logistic regression failed: {e}")
+    
+    # Breslow-Day test for homogeneity of odds ratios across oversight strata
+    print_subsection("4b. Breslow-Day test for homogeneity of odds ratios")
+    print("  (Tests whether the scenario effect is consistent across oversight levels)")
+    
+    ors = []
+    for oversight in ["low", "mid", "high"]:
+        odf = df[df["oversight"] == oversight]
+        ct = pd.crosstab(odf["scenario"], odf["is_deceptive_output"])
+        if ct.shape == (2, 2):
+            or_val, ci_low, ci_high = odds_ratio_ci(ct.values)
+            ors.append((oversight, or_val, ci_low, ci_high))
+            print(f"    {oversight}: OR = {or_val:.3f}  95% CI: [{ci_low:.3f}, {ci_high:.3f}]")
+    
+    # Manual Breslow-Day: compare ORs across strata
+    # Use Cochran-Mantel-Haenszel common OR for comparison
+    tables = []
+    for oversight in ["low", "mid", "high"]:
+        odf = df[df["oversight"] == oversight]
+        ct = pd.crosstab(odf["scenario"], odf["is_deceptive_output"])
+        if ct.shape == (2, 2):
+            tables.append(ct.values)
+    
+    if len(tables) == 3:
+        # CMH test
+        from statsmodels.stats.contingency_tables import StratifiedTable
+        strat = StratifiedTable(tables)
+        cmh_result = strat.test_null_odds()
+        print(f"\n  Cochran-Mantel-Haenszel test:")
+        print(f"    Common OR = {strat.summary().tables[0].data[1][1] if hasattr(strat.summary(), 'tables') else 'N/A'}")
+        print(f"    CMH statistic = {cmh_result.statistic:.3f}, {format_p(cmh_result.pvalue)} {sig_stars(cmh_result.pvalue)}")
+        
+        bd_result = strat.test_equal_odds()
+        print(f"\n  Breslow-Day test for homogeneity of ORs:")
+        print(f"    Statistic = {bd_result.statistic:.3f}, {format_p(bd_result.pvalue)} {sig_stars(bd_result.pvalue)}")
+        if bd_result.pvalue < 0.05:
+            print("    → ORs differ significantly across oversight strata")
+            print("      (The scenario effect is NOT consistent across oversight levels)")
+        else:
+            print("    → ORs are consistent across oversight strata")
+
+
+def analysis_5_model_x_oversight_within_scenario(df: pd.DataFrame):
+    """Test H5: Oversight effect varies by model within each scenario."""
+    print_section("5. MODEL × OVERSIGHT INTERACTION (per scenario)")
+    
+    for scenario in sorted(df["scenario"].unique()):
+        print_subsection(f"5. Oversight effect by model — {scenario}")
+        
+        sdf = df[df["scenario"] == scenario]
+        
+        for model in sorted(sdf["model"].unique()):
+            mdf = sdf[sdf["model"] == model]
+            
+            counts = []
+            nobs = []
+            rates = []
+            for oversight in ["low", "mid", "high"]:
+                odf = mdf[mdf["oversight"] == oversight]
+                dec = odf["is_deceptive_output"].sum()
+                tot = len(odf)
+                counts.append(dec)
+                nobs.append(tot)
+                rates.append(dec / tot * 100 if tot > 0 else 0)
+            
+            z_trend, p_trend = cochran_armitage_trend_test(counts, nobs)
+            trend_sig = sig_stars(p_trend)
+            
+            model_short = model.split("/")[-1] if "/" in model else model
+            print(f"    {model_short:30s}  low={rates[0]:5.1f}%  mid={rates[1]:5.1f}%  "
+                  f"high={rates[2]:5.1f}%  "
+                  f"trend z={z_trend:+.2f} {format_p(p_trend)} {trend_sig}")
+
+
+def analysis_6_deception_type(df: pd.DataFrame):
+    """Test H6: Deception type distribution differs between scenarios."""
+    print_section("6. DECEPTION TYPE ANALYSIS")
+    
+    deceptive_df = df[df["blackbox_category"].isin(DECEPTIVE_CATEGORIES)]
+    
+    print_subsection("6a. Deception type distribution by scenario")
+    ct = pd.crosstab(deceptive_df["scenario"], deceptive_df["blackbox_category"])
+    
+    for scenario in ct.index:
+        total = ct.loc[scenario].sum()
+        print(f"  {scenario} (n={total}):")
+        for cat in sorted(ct.columns):
+            count = ct.loc[scenario, cat]
+            rate = count / total * 100
+            print(f"    {cat}: {count} ({rate:.1f}%)")
+    
+    chi2_val, p_val, dof, expected = chi2_contingency(ct.values)
+    v = cramers_v(ct.values)
+    print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+    print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+    
+    min_expected = expected.min()
+    if min_expected < 5:
+        print(f"  ⚠ Minimum expected count = {min_expected:.1f} < 5, interpreting with caution")
+
+
+def analysis_7_glassbox_categories(df: pd.DataFrame):
+    """Test H7: Glassbox intent categories differ between scenarios."""
+    print_section("7. REASONING INTENT ANALYSIS (Glassbox)")
+    
+    print_subsection("7a. Glassbox category distribution by scenario")
+    ct = pd.crosstab(df["scenario"], df["glassbox_category"])
+    
+    for scenario in ct.index:
+        total = ct.loc[scenario].sum()
+        print(f"  {scenario} (n={total}):")
+        for cat in sorted(ct.columns):
+            count = ct.loc[scenario, cat]
+            rate = count / total * 100
+            ci_low, ci_high = proportion_confint(count, total, alpha=0.05, method="wilson")
+            print(f"    {cat}: {count} ({rate:.1f}%) "
+                  f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+    
+    chi2_val, p_val, dof, expected = chi2_contingency(ct.values)
+    v = cramers_v(ct.values)
+    print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+    print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+    
+    min_expected = expected.min()
+    if min_expected < 5:
+        print(f"  ⚠ Minimum expected count = {min_expected:.1f} < 5")
+    
+    # Sophistication distribution
+    print_subsection("7b. Sophistication distribution by scenario")
+    ct_soph = pd.crosstab(df["scenario"], df["glassbox_sophistication"])
+    
+    for scenario in ct_soph.index:
+        total = ct_soph.loc[scenario].sum()
+        print(f"  {scenario} (n={total}):")
+        for tier in ["NAIVE", "REACTIVE", "TACTICAL", "STRATEGIC"]:
+            if tier in ct_soph.columns:
+                count = ct_soph.loc[scenario, tier]
+                rate = count / total * 100
+                print(f"    {tier}: {count} ({rate:.1f}%)")
+    
+    # Filter to valid tiers only for test
+    valid_tiers = [t for t in ["NAIVE", "REACTIVE", "TACTICAL", "STRATEGIC"] if t in ct_soph.columns]
+    ct_soph_valid = ct_soph[valid_tiers]
+    
+    if ct_soph_valid.shape[1] > 1:
+        chi2_val, p_val, dof, expected = chi2_contingency(ct_soph_valid.values)
+        v = cramers_v(ct_soph_valid.values)
+        print(f"\n  Chi-squared test: χ²({dof}) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+        print(f"  Cramér's V = {v:.3f} ({interpret_effect_size(v)})")
+
+
+def analysis_8_reasoning_output_alignment(df: pd.DataFrame):
+    """Test H8: Alignment between reasoning intent and output deception."""
+    print_section("8. REASONING–OUTPUT ALIGNMENT")
+    
+    for scenario in sorted(df["scenario"].unique()):
+        print_subsection(f"8. Alignment matrix — {scenario}")
+        
+        sdf = df[df["scenario"] == scenario]
+        ct = pd.crosstab(sdf["glassbox_category"], sdf["blackbox_category"])
+        
+        # Print nicely
+        all_bb = ["ACCURATE", "MISLEADING", "OMITTED", "FABRICATED"]
+        all_gb = ["ROBUST_HONEST", "CONSIDERED_HONEST", "COMPLIANT_DECEPTION", "FAILED_DECEPTION", "UNINTENTIONAL_FAILURE"]
+        
+        present_bb = [c for c in all_bb if c in ct.columns]
+        present_gb = [g for g in all_gb if g in ct.index]
+        
+        # Header
+        header = f"    {'':30s}" + "".join(f"  {c:>12s}" for c in present_bb) + f"  {'TOTAL':>8s}"
+        print(header)
+        print("    " + "-" * (len(header) - 4))
+        
+        for gb in present_gb:
+            row_vals = []
+            for bb in present_bb:
+                val = ct.loc[gb, bb] if bb in ct.columns and gb in ct.index else 0
+                row_vals.append(val)
+            total = sum(row_vals)
+            row_str = f"    {gb:30s}" + "".join(f"  {v:>12d}" for v in row_vals) + f"  {total:>8d}"
+            print(row_str)
+        
+        # Key misalignment metrics
+        compliant = sdf[sdf["glassbox_category"] == "COMPLIANT_DECEPTION"]
+        if len(compliant) > 0:
+            accurate_despite_compliant = (compliant["blackbox_category"] == "ACCURATE").sum()
+            rate = accurate_despite_compliant / len(compliant) * 100
+            ci_low, ci_high = proportion_confint(accurate_despite_compliant, len(compliant), alpha=0.05, method="wilson")
+            print(f"\n  Compliant reasoning → accurate output: {accurate_despite_compliant}/{len(compliant)} ({rate:.1f}%)"
+                  f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+        
+        honest_cats = sdf[sdf["glassbox_category"].isin(["ROBUST_HONEST", "CONSIDERED_HONEST"])]
+        if len(honest_cats) > 0:
+            deceptive_despite_honest = honest_cats["is_deceptive_output"].sum()
+            rate = deceptive_despite_honest / len(honest_cats) * 100
+            ci_low, ci_high = proportion_confint(deceptive_despite_honest, len(honest_cats), alpha=0.05, method="wilson")
+            print(f"  Honest reasoning → deceptive output: {deceptive_despite_honest}/{len(honest_cats)} ({rate:.1f}%)"
+                  f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+
+
+def analysis_9_safety_training(df: pd.DataFrame):
+    """Test H9: Safety-trained model vs base model comparison."""
+    print_section("9. SAFETY TRAINING EFFECT (gpt-oss-20b vs gpt-oss-safeguard-20b)")
+    
+    base_name = "openai/gpt-oss-20b"
+    safe_name = "openai/gpt-oss-safeguard-20b"
+    
+    pair_df = df[df["model"].isin([base_name, safe_name])]
+    
+    if len(pair_df) == 0:
+        print("  No data found for safety training comparison.")
+        return
+    
+    for scenario in sorted(pair_df["scenario"].unique()):
+        print_subsection(f"9. Safety training effect — {scenario}")
+        
+        sdf = pair_df[pair_df["scenario"] == scenario]
+        ct = pd.crosstab(sdf["model"], sdf["is_deceptive_output"])
+        
+        for model in [base_name, safe_name]:
+            if model in ct.index:
+                n_honest = ct.loc[model, 0] if 0 in ct.columns else 0
+                n_deceptive = ct.loc[model, 1] if 1 in ct.columns else 0
+                total = n_honest + n_deceptive
+                rate = n_deceptive / total * 100
+                ci_low, ci_high = proportion_confint(n_deceptive, total, alpha=0.05, method="wilson")
+                label = "base" if "safeguard" not in model else "safeguard"
+                print(f"    {label}: {n_deceptive}/{total} ({rate:.1f}%) "
+                      f"  95% CI: [{ci_low*100:.1f}%, {ci_high*100:.1f}%]")
+        
+        if ct.shape == (2, 2):
+            chi2_val, p_val, dof, _ = chi2_contingency(ct.values)
+            _, fisher_p = fisher_exact(ct.values)
+            or_val, or_ci_low, or_ci_high = odds_ratio_ci(ct.values)
+            
+            print(f"\n  Chi-squared: χ²(1) = {chi2_val:.3f}, {format_p(p_val)} {sig_stars(p_val)}")
+            print(f"  Fisher's exact: {format_p(fisher_p)} {sig_stars(fisher_p)}")
+            print(f"  Odds ratio = {or_val:.3f}  95% CI: [{or_ci_low:.3f}, {or_ci_high:.3f}]")
+        else:
+            print(f"\n  ⚠ Cannot compute 2x2 test — table shape is {ct.shape}")
+
+
+def analysis_10_full_logistic_model(df: pd.DataFrame):
+    """Run a full logistic regression with all factors."""
+    print_section("10. FULL LOGISTIC REGRESSION MODEL")
+    
+    print_subsection("10a. Deceptive output ~ scenario + oversight + model")
+    
+    model_df = df.copy()
+    model_df["scenario_code"] = (model_df["scenario"] == "child_protection").astype(int)
+    model_df["oversight_num"] = model_df["oversight_ordinal"]
+    
+    # Create model dummies
+    model_dummies = pd.get_dummies(model_df["model"], prefix="model", drop_first=True, dtype=float)
+    
+    X = pd.concat([
+        model_df[["scenario_code", "oversight_num"]],
+        model_dummies
+    ], axis=1)
+    X = sm.add_constant(X)
+    y = model_df["is_deceptive_output"]
+    
+    try:
+        logit_model = sm.Logit(y, X).fit(disp=0, maxiter=100)
+        print(logit_model.summary2().tables[1].to_string())
+        print(f"\n  Pseudo R² = {logit_model.prsquared:.4f}")
+        print(f"  AIC = {logit_model.aic:.1f}")
+        print(f"  n = {logit_model.nobs:.0f}")
+        
+        # Print odds ratios for interpretability
+        print("\n  Odds ratios (exponentiated coefficients):")
+        for name, coef, pval in zip(X.columns, logit_model.params, logit_model.pvalues):
+            if name == "const":
+                continue
+            or_val = np.exp(coef)
+            ci = logit_model.conf_int().loc[name]
+            or_ci = (np.exp(ci[0]), np.exp(ci[1]))
+            print(f"    {name:45s}  OR = {or_val:.3f}  "
+                  f"95% CI: [{or_ci[0]:.3f}, {or_ci[1]:.3f}]  {sig_stars(pval)}")
+            
+    except Exception as e:
+        print(f"  Logistic regression failed: {e}")
+    
+    # Also the interaction model
+    print_subsection("10b. With scenario × oversight interaction")
+    
+    model_df["interaction"] = model_df["scenario_code"] * model_df["oversight_num"]
+    X2 = pd.concat([
+        model_df[["scenario_code", "oversight_num", "interaction"]],
+        model_dummies
+    ], axis=1)
+    X2 = sm.add_constant(X2)
+    
+    try:
+        logit_interaction = sm.Logit(y, X2).fit(disp=0, maxiter=100)
+        print(logit_interaction.summary2().tables[1].to_string())
+        
+        # LR test for interaction
+        lr_stat = -2 * (logit_model.llf - logit_interaction.llf)
+        lr_p = chi2.sf(lr_stat, df=1)
+        print(f"\n  LR test for interaction: LR = {lr_stat:.3f}, {format_p(lr_p)} {sig_stars(lr_p)}")
+    except Exception as e:
+        print(f"  Interaction model failed: {e}")
+
+
+def analysis_11_summary_table(df: pd.DataFrame):
+    """Print a publication-ready summary table."""
+    print_section("11. SUMMARY TABLE (Publication-Ready)")
+    
+    print(f"  {'Model':35s} | {'Scenario':25s} | {'Low':>10s} | {'Mid':>10s} | {'High':>10s} | {'Overall':>10s}")
+    print("  " + "-" * 105)
+    
+    for model in sorted(df["model"].unique()):
+        for scenario in sorted(df["scenario"].unique()):
+            rates = []
+            for oversight in ["low", "mid", "high"]:
+                sub = df[(df["model"] == model) & (df["scenario"] == scenario) & (df["oversight"] == oversight)]
+                if len(sub) > 0:
+                    rate = sub["is_deceptive_output"].mean() * 100
+                    rates.append(f"{rate:.0f}%")
+                else:
+                    rates.append("—")
+            
+            overall_sub = df[(df["model"] == model) & (df["scenario"] == scenario)]
+            overall_rate = overall_sub["is_deceptive_output"].mean() * 100
+            
+            model_display = model.split("/")[-1] if "/" in model else model
+            print(f"  {model_display:35s} | {scenario:25s} | {rates[0]:>10s} | {rates[1]:>10s} | {rates[2]:>10s} | {overall_rate:>9.0f}%")
+        print("  " + "-" * 105)
+    
+    # Marginals
+    print(f"\n  {'MARGINAL':35s} | {'':25s} | {'Low':>10s} | {'Mid':>10s} | {'High':>10s} | {'Overall':>10s}")
+    print("  " + "-" * 105)
+    for scenario in sorted(df["scenario"].unique()):
+        rates = []
+        for oversight in ["low", "mid", "high"]:
+            sub = df[(df["scenario"] == scenario) & (df["oversight"] == oversight)]
+            rate = sub["is_deceptive_output"].mean() * 100
+            rates.append(f"{rate:.0f}%")
+        overall = df[df["scenario"] == scenario]["is_deceptive_output"].mean() * 100
+        print(f"  {'ALL MODELS':35s} | {scenario:25s} | {rates[0]:>10s} | {rates[1]:>10s} | {rates[2]:>10s} | {overall:>9.0f}%")
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+
+def main():
+    if len(sys.argv) < 2:
+        csv_path = "logs/full_experiment/results.csv"
+    else:
+        csv_path = sys.argv[1]
+    
+    print(f"Loading data from: {csv_path}")
+    df = load_data(csv_path)
+    print(f"Loaded {len(df)} rows: "
+          f"{df['model'].nunique()} models, "
+          f"{df['scenario'].nunique()} scenarios, "
+          f"{df['oversight'].nunique()} oversight levels")
+    print(f"Significance level: α = {SIGNIFICANCE_LEVEL}")
+    
+    analysis_1_scenario_effect(df)
+    analysis_2_oversight_effect(df)
+    analysis_3_model_effects(df)
+    analysis_4_oversight_x_scenario_interaction(df)
+    analysis_5_model_x_oversight_within_scenario(df)
+    analysis_6_deception_type(df)
+    analysis_7_glassbox_categories(df)
+    analysis_8_reasoning_output_alignment(df)
+    analysis_9_safety_training(df)
+    analysis_10_full_logistic_model(df)
+    analysis_11_summary_table(df)
+    
+    print_section("ANALYSIS COMPLETE")
+    print(f"  Total runs analysed: {len(df)}")
+    print(f"  Models: {sorted(df['model'].unique())}")
+    print(f"  Scenarios: {sorted(df['scenario'].unique())}")
+    print(f"  Oversight levels: {sorted(df['oversight'].unique())}")
+    print(f"  Judge categories used:")
+    print(f"    Blackbox: {DECEPTIVE_CATEGORIES} vs {HONEST_CATEGORIES}")
+    print(f"    Glassbox: {sorted(df['glassbox_category'].unique())}")
+    print(f"    Sophistication: {sorted(df['glassbox_sophistication'].unique())}")
+
+
+if __name__ == "__main__":
+    main()