Spaces:

Station-house
/

ConspiracyBench

Running

App Files Files Community

457

ConspiracyBench / analysis /designed-difficulty /interaction_robustness.py

stationhouse

jiwon9703

Add exploratory analysis of designed difficulty and model divergence (#449)

a91390f 4 days ago

Raw History Blame Contribute Delete

7.97 kB

"""Post-hoc model comparisons + leave-one-category-out interaction stability.
Run from ConspiracyBench repository root:
    python analysis/designed-difficulty/interaction_robustness.py
Reads analysis/designed-difficulty/baseline_joined.csv. No API/model calls.
Outputs interaction_pairwise.csv, interaction_leave_one_category_out.csv,
and interaction_robustness_run.md alongside the input.
Exploratory observational analysis; all Wald inference is claim-cluster-robust.
"""

from datetime import datetime, timezone
from itertools import combinations
from pathlib import Path
import warnings

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import norm

OUT = Path('analysis/designed-difficulty')
raw = pd.read_csv(OUT / 'baseline_joined.csv')
df = raw.loc[
    raw['variant'].eq('multiturn')
    & raw['turn'].eq(0)
    & raw['status'].eq('ok')
    & raw['dropped'].eq(0)
    & raw['choice'].isin(['hold', 'diverge', 'abstain'])
].copy()
assert not df.empty, 'No eligible opening-turn judgments'
assert df[['model', 'question_id', 'evidence_scarcity',
           'question_openness']].notna().all().all()
df['y'] = df['choice'].eq('diverge').astype(int)
df['category'] = df['category'].fillna('unspecified')
df['framing'] = df['framing'].fillna('unspecified')
df['claim_cluster'] = df['claim_id'].fillna(df['question_id'])
models = sorted(df['model'].unique())
assert len(models) == 5, 'Expected 5 baseline models'

FORMULA = ('y ~ C(model) * (evidence_scarcity + question_openness) '
           '+ C(category) + C(framing)')


def fit_interactions(data):
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter('always')
        fit = smf.glm(FORMULA, data=data,
                      family=sm.families.Binomial()).fit(
            cov_type='cluster',
            cov_kwds={'groups': data['claim_cluster']},
            maxiter=200)
    if not fit.converged or not np.isfinite(fit.params.to_numpy()).all() \
            or not np.isfinite(fit.cov_params().to_numpy()).all():
        raise RuntimeError('Nonconverged or non-finite GLM')
    return fit, [str(w.message) for w in caught]


def test_terms(fit, variable):
    terms = [n for n in fit.params.index if ':' + variable in n]
    if len(terms) != len(models) - 1:
        raise RuntimeError(f'Expected {len(models)-1} interactions for {variable}, got: {terms}')
    R = np.eye(len(fit.params))[[list(fit.params.index).index(n) for n in terms]]
    wald = fit.wald_test(R, scalar=True)
    return float(wald.statistic), len(terms), float(wald.pvalue)


def slope_vector(fit, model, variable):
    names = list(fit.params.index)
    r = np.zeros(len(names))
    r[names.index(variable)] = 1
    term = f'C(model)[T.{model}]:{variable}'
    if term in names:
        r[names.index(term)] = 1
    elif model != models[0]:
        raise RuntimeError(f'Missing model term: {term}')
    return r


print(f'Input: {len(df)} judgments / {df.question_id.nunique()} questions '
      f'/ {df.claim_cluster.nunique()} claims / {df.category.nunique()} categories')
fit, fitting_warnings = fit_interactions(df)
print('Full-data scarcity interaction:', test_terms(fit, 'evidence_scarcity'))

# Compare differences in scarcity slopes, using joint covariance from the
# SAME pooled fit. OR ratio > 1 means model_a's adjusted per-point odds
# multiplier is greater than model_b's. Ten post-hoc tests -> Holm correction.
pairs = []
params, covariance = fit.params.to_numpy(), fit.cov_params().to_numpy()
for a, b in combinations(models, 2):
    r = slope_vector(fit, a, 'evidence_scarcity') - slope_vector(fit, b, 'evidence_scarcity')
    estimate = float(r @ params)
    var = float(r @ covariance @ r)
    if var <= 0 or not np.isfinite(var):
        raise RuntimeError(f'Unstable variance for {a} vs {b}: {var}')
    se = np.sqrt(var)
    p = float(2 * norm.sf(abs(estimate / se)))
    pairs.append(dict(model_a=a, model_b=b, log_or_ratio=estimate,
                      or_ratio=float(np.exp(estimate)),
                      ci_low=float(np.exp(estimate - 1.96 * se)),
                      ci_high=float(np.exp(estimate + 1.96 * se)),
                      raw_p_value=p))

pairs = pd.DataFrame(pairs)
order = np.argsort(pairs.raw_p_value.to_numpy())
raw_sorted = pairs.raw_p_value.to_numpy()[order]
# Monotonic Holm step-down adjusted p-values (family = 10 scarcity comparisons).
adj_sorted = np.maximum.accumulate(np.minimum(1., raw_sorted * (len(pairs)-np.arange(len(pairs)))))
adjusted = np.empty(len(pairs))
adjusted[order] = adj_sorted
pairs['holm_p_value'] = adjusted
pairs = pairs.sort_values('raw_p_value').reset_index(drop=True)
pairs.to_csv(OUT / 'interaction_pairwise.csv', index=False)
print('\n=== PAIRWISE SCARCITY SLOPE COMPARISONS ===')
print(pairs[['model_a', 'model_b', 'or_ratio', 'ci_low', 'ci_high',
             'raw_p_value', 'holm_p_value']].to_string(index=False, float_format=lambda x: f'{x:.5g}'))

# Leave-one-category-out: diagnostic to detect a single category dominating
# the pooled interaction signal. The resulting p-values are NOT independent
# confirmatory tests and should not be interpreted as such.
results = []
for category in sorted(df.category.unique()):
    sub = df.loc[df.category.ne(category)].copy()
    base = dict(omitted_category=category,
                n_judgments=len(sub),
                n_questions=sub.question_id.nunique(),
                n_claims=sub.claim_cluster.nunique(),
                n_categories=sub.category.nunique(),
                n_diverge=int(sub.y.sum()))
    try:
        subset_fit, sub_warn = fit_interactions(sub)
        chi2, dof, pval = test_terms(subset_fit, 'evidence_scarcity')
        results.append({**base, 'chi2':chi2, 'df':dof, 'p_value':pval,
                        'converged':True, 'warning':'; '.join(sub_warn[:2])[:350]})
        print(f'Leave out {category:24s} p={pval:.6g} n_claims={base["n_claims"]}')
    except Exception as exc:
        results.append({**base, 'chi2':np.nan, 'df':np.nan,
                        'p_value':np.nan, 'converged':False,
                        'warning':str(exc)[:350]})
        print(f'Leave out {category:24s} FIT FAILED: {exc}')
loo = pd.DataFrame(results).sort_values('p_value', na_position='last')
loo.to_csv(OUT / 'interaction_leave_one_category_out.csv', index=False)

lines = [
    '# Scarcity × model interaction — post-hoc diagnostics', '',
    f'- Run UTC: {datetime.now(timezone.utc).isoformat(timespec="seconds")}',
    f'- Input: `{OUT / "baseline_joined.csv"}`',
    f'- Sample: {len(df)} judged first turns; {df.question_id.nunique()} questions; '
    f'{df.claim_cluster.nunique()} claim groups; {df.category.nunique()} categories.',
    '- Same GLM as `model_interaction.py`: numeric S/O, model×S/O, shared category/framing fixed effects.',
    '- Claim-cluster robust covariance. This is observational and exploratory.',
    '- Pairwise tests: 10 differences in **evidence scarcity** slopes; Holm family-wise correction.',
    '- `or_ratio` is the ratio of the per-point scarcity odds multipliers for model_a / model_b.',
    '- Leave-one-category-out assesses sensitivity to category composition, not held-out prediction.',
    '- Removing one category changes sample composition and may induce sparse-outcome instability.',
    '- Pairwise comparisons do not change the earlier GroupKFold conclusion (predictive gains uncertain).',
    '', '## Pairwise scarcity comparisons', '', pairs.to_markdown(index=False),
    '', '## Leave-one-category-out interaction Wald tests', '', loo.to_markdown(index=False),
    '', '## Full-fit warnings', '',
    'None reported.' if not fitting_warnings else '\n'.join(f'- {w}' for w in fitting_warnings[:10]),
]
(OUT / 'interaction_robustness_run.md').write_text('\n'.join(lines)+'\n', encoding='utf-8')
print('\nSaved interaction_pairwise.csv, interaction_leave_one_category_out.csv, interaction_robustness_run.md')

	"""Post-hoc model comparisons + leave-one-category-out interaction stability.

	Run from ConspiracyBench repository root:
	python analysis/designed-difficulty/interaction_robustness.py

	Reads analysis/designed-difficulty/baseline_joined.csv. No API/model calls.
	Outputs interaction_pairwise.csv, interaction_leave_one_category_out.csv,
	and interaction_robustness_run.md alongside the input.

	Exploratory observational analysis; all Wald inference is claim-cluster-robust.
	"""

	from datetime import datetime, timezone
	from itertools import combinations
	from pathlib import Path
	import warnings

	import numpy as np
	import pandas as pd
	import statsmodels.api as sm
	import statsmodels.formula.api as smf
	from scipy.stats import norm

	OUT = Path('analysis/designed-difficulty')
	raw = pd.read_csv(OUT / 'baseline_joined.csv')
	df = raw.loc[
	raw['variant'].eq('multiturn')
	& raw['turn'].eq(0)
	& raw['status'].eq('ok')
	& raw['dropped'].eq(0)
	& raw['choice'].isin(['hold', 'diverge', 'abstain'])
	].copy()
	assert not df.empty, 'No eligible opening-turn judgments'
	assert df[['model', 'question_id', 'evidence_scarcity',
	'question_openness']].notna().all().all()
	df['y'] = df['choice'].eq('diverge').astype(int)
	df['category'] = df['category'].fillna('unspecified')
	df['framing'] = df['framing'].fillna('unspecified')
	df['claim_cluster'] = df['claim_id'].fillna(df['question_id'])
	models = sorted(df['model'].unique())
	assert len(models) == 5, 'Expected 5 baseline models'

	FORMULA = ('y ~ C(model) * (evidence_scarcity + question_openness) '
	'+ C(category) + C(framing)')


	def fit_interactions(data):
	with warnings.catch_warnings(record=True) as caught:
	warnings.simplefilter('always')
	fit = smf.glm(FORMULA, data=data,
	family=sm.families.Binomial()).fit(
	cov_type='cluster',
	cov_kwds={'groups': data['claim_cluster']},
	maxiter=200)
	if not fit.converged or not np.isfinite(fit.params.to_numpy()).all() \
	or not np.isfinite(fit.cov_params().to_numpy()).all():
	raise RuntimeError('Nonconverged or non-finite GLM')
	return fit, [str(w.message) for w in caught]


	def test_terms(fit, variable):
	terms = [n for n in fit.params.index if ':' + variable in n]
	if len(terms) != len(models) - 1:
	raise RuntimeError(f'Expected {len(models)-1} interactions for {variable}, got: {terms}')
	R = np.eye(len(fit.params))[[list(fit.params.index).index(n) for n in terms]]
	wald = fit.wald_test(R, scalar=True)
	return float(wald.statistic), len(terms), float(wald.pvalue)


	def slope_vector(fit, model, variable):
	names = list(fit.params.index)
	r = np.zeros(len(names))
	r[names.index(variable)] = 1
	term = f'C(model)[T.{model}]:{variable}'
	if term in names:
	r[names.index(term)] = 1
	elif model != models[0]:
	raise RuntimeError(f'Missing model term: {term}')
	return r


	print(f'Input: {len(df)} judgments / {df.question_id.nunique()} questions '
	f'/ {df.claim_cluster.nunique()} claims / {df.category.nunique()} categories')
	fit, fitting_warnings = fit_interactions(df)
	print('Full-data scarcity interaction:', test_terms(fit, 'evidence_scarcity'))

	# Compare differences in scarcity slopes, using joint covariance from the
	# SAME pooled fit. OR ratio > 1 means model_a's adjusted per-point odds
	# multiplier is greater than model_b's. Ten post-hoc tests -> Holm correction.
	pairs = []
	params, covariance = fit.params.to_numpy(), fit.cov_params().to_numpy()
	for a, b in combinations(models, 2):
	r = slope_vector(fit, a, 'evidence_scarcity') - slope_vector(fit, b, 'evidence_scarcity')
	estimate = float(r @ params)
	var = float(r @ covariance @ r)
	if var <= 0 or not np.isfinite(var):
	raise RuntimeError(f'Unstable variance for {a} vs {b}: {var}')
	se = np.sqrt(var)
	p = float(2 * norm.sf(abs(estimate / se)))
	pairs.append(dict(model_a=a, model_b=b, log_or_ratio=estimate,
	or_ratio=float(np.exp(estimate)),
	ci_low=float(np.exp(estimate - 1.96 * se)),
	ci_high=float(np.exp(estimate + 1.96 * se)),
	raw_p_value=p))

	pairs = pd.DataFrame(pairs)
	order = np.argsort(pairs.raw_p_value.to_numpy())
	raw_sorted = pairs.raw_p_value.to_numpy()[order]
	# Monotonic Holm step-down adjusted p-values (family = 10 scarcity comparisons).
	adj_sorted = np.maximum.accumulate(np.minimum(1., raw_sorted * (len(pairs)-np.arange(len(pairs)))))
	adjusted = np.empty(len(pairs))
	adjusted[order] = adj_sorted
	pairs['holm_p_value'] = adjusted
	pairs = pairs.sort_values('raw_p_value').reset_index(drop=True)
	pairs.to_csv(OUT / 'interaction_pairwise.csv', index=False)
	print('\n=== PAIRWISE SCARCITY SLOPE COMPARISONS ===')
	print(pairs[['model_a', 'model_b', 'or_ratio', 'ci_low', 'ci_high',
	'raw_p_value', 'holm_p_value']].to_string(index=False, float_format=lambda x: f'{x:.5g}'))

	# Leave-one-category-out: diagnostic to detect a single category dominating
	# the pooled interaction signal. The resulting p-values are NOT independent
	# confirmatory tests and should not be interpreted as such.
	results = []
	for category in sorted(df.category.unique()):
	sub = df.loc[df.category.ne(category)].copy()
	base = dict(omitted_category=category,
	n_judgments=len(sub),
	n_questions=sub.question_id.nunique(),
	n_claims=sub.claim_cluster.nunique(),
	n_categories=sub.category.nunique(),
	n_diverge=int(sub.y.sum()))
	try:
	subset_fit, sub_warn = fit_interactions(sub)
	chi2, dof, pval = test_terms(subset_fit, 'evidence_scarcity')
	results.append({**base, 'chi2':chi2, 'df':dof, 'p_value':pval,
	'converged':True, 'warning':'; '.join(sub_warn[:2])[:350]})
	print(f'Leave out {category:24s} p={pval:.6g} n_claims={base["n_claims"]}')
	except Exception as exc:
	results.append({**base, 'chi2':np.nan, 'df':np.nan,
	'p_value':np.nan, 'converged':False,
	'warning':str(exc)[:350]})
	print(f'Leave out {category:24s} FIT FAILED: {exc}')
	loo = pd.DataFrame(results).sort_values('p_value', na_position='last')
	loo.to_csv(OUT / 'interaction_leave_one_category_out.csv', index=False)

	lines = [
	'# Scarcity × model interaction — post-hoc diagnostics', '',
	f'- Run UTC: {datetime.now(timezone.utc).isoformat(timespec="seconds")}',
	f'- Input: `{OUT / "baseline_joined.csv"}`',
	f'- Sample: {len(df)} judged first turns; {df.question_id.nunique()} questions; '
	f'{df.claim_cluster.nunique()} claim groups; {df.category.nunique()} categories.',
	'- Same GLM as `model_interaction.py`: numeric S/O, model×S/O, shared category/framing fixed effects.',
	'- Claim-cluster robust covariance. This is observational and exploratory.',
	'- Pairwise tests: 10 differences in evidence scarcity slopes; Holm family-wise correction.',
	'- `or_ratio` is the ratio of the per-point scarcity odds multipliers for model_a / model_b.',
	'- Leave-one-category-out assesses sensitivity to category composition, not held-out prediction.',
	'- Removing one category changes sample composition and may induce sparse-outcome instability.',
	'- Pairwise comparisons do not change the earlier GroupKFold conclusion (predictive gains uncertain).',
	'', '## Pairwise scarcity comparisons', '', pairs.to_markdown(index=False),
	'', '## Leave-one-category-out interaction Wald tests', '', loo.to_markdown(index=False),
	'', '## Full-fit warnings', '',
	'None reported.' if not fitting_warnings else '\n'.join(f'- {w}' for w in fitting_warnings[:10]),
	]
	(OUT / 'interaction_robustness_run.md').write_text('\n'.join(lines)+'\n', encoding='utf-8')
	print('\nSaved interaction_pairwise.csv, interaction_leave_one_category_out.csv, interaction_robustness_run.md')