Metadata-Version: 2.4
Name: confound-controls
Version: 0.4.0
Summary: Confound-matched controls for classifier claims: 1:1 matching, bootstrap AUROC, anchored verdicts, and a refusal when the control is vacuous.
Author: Jaret Arnold
License: GPL-3.0-or-later
Project-URL: Homepage, https://github.com/musharna/confound-controls
Project-URL: Repository, https://github.com/musharna/confound-controls
Project-URL: Changelog, https://github.com/musharna/confound-controls/blob/master/CHANGELOG.md
Project-URL: Issues, https://github.com/musharna/confound-controls/issues
Keywords: confounding,matching,auroc,bootstrap,classifier-evaluation,research-integrity
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: scipy>=1.10
Requires-Dist: scikit-learn>=1.3
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Requires-Dist: ruff>=0.6; extra == "dev"
Dynamic: license-file

# confound-controls

Did your classifier learn the signal, or a confound?

For each confound, restrict the negatives to a subset matched 1:1 to the
positives on that confound, then recompute AUROC. If the signal survives
matching, it is not that confound. If it collapses to chance, it was.

```python
from confound_controls import run_battery, format_battery, bootstrap_auroc

anchor = bootstrap_auroc(df["label"], probs).point  # uncontrolled baseline

results = run_battery(
    df,
    prob_map,
    spec={"gc": ["gc"], "expression": ["log_expr"], "joint": ["gc", "log_expr"]},
    anchor=anchor,
)
print(format_battery(results, anchor))
```

Nothing is read from disk and nothing is written to it. The frame, the
probabilities, the confound spec and the anchor are all arguments.

## The failure this exists to prevent

**1:1 matching without replacement from a negative pool the same size as the
positive set consumes every negative, whatever the confound values are.** The
matched comparison is then identical to the unmatched one — and a vacuous
control does not announce itself as inconclusive. It returns the original
result, which reads as the confound having been ruled out.

This is not hypothetical. It is how the first version of this library's own
test suite behaved: 150 positives against 150 negatives, a score computed from
nothing but GC, matched on GC, verdict `confound-robust`, AUROC 0.9312.

So `evaluate_confound` refuses by default:

```
ValueError: matching consumed the entire negative pool (150 matched from 150
available), so it selected nothing and the control is vacuous -- the matched
comparison equals the unmatched one and will report the confound as ruled out.
```

Pass `require_selective_match=False` if you mean it; the result then carries
`match_selective=False` and `format_battery` prints `[VACUOUS CONTROL]`.

## The same failure, one size up

A pool bigger than the positive set passes that check and can fail the same
way. Greedy matching spends the negatives that sit near the positives first;
the rest of the positives are paired with whatever is left. The match is
complete and selective, and still confounded. Through 0.3.2 this library
scored it anyway: 300 positives, 900 negatives, a score computed from nothing
but GC, matched on GC — verdict `confound-robust`, AUROC 0.807, battery PASS.

`evaluate_confound` now measures what the match achieved and refuses when it
achieved too little:

```
ValueError: gc: matching did not balance 'gc': |SMD| = 0.987 after matching
(was 2.010) exceeds 0.1. The matched set is still confounded, so a classifier
that separates it has not been shown to be robust to this confound.
```

The standardized mean difference is the difference in group means in units of
a standard deviation, measured against the unmatched groups' spread before and
after. `max_smd=0.1` is the conventional line; `max_smd=None` reports instead of
raising (`result.balance`, and `[IMBALANCED]` in `format_battery`).

Balance is bought with a **caliper**: a positive with no negative within
`caliper` standard deviations is left unmatched and reported, rather than paired
with something that does not resemble it. Only matched positives are scored.
The same data with `caliper=0.1`:

```
gc   AUROC=0.5015 CI[0.4350,0.5652] recovery=0.00 npos=148 -> confound-driven
     SMD before 2.010   SMD after 0.005   var ratio 1.017
```

With many confounds every sample is far from every other, and nearest-neighbour
matching stops finding anyone. `method="propensity"` matches on the logit of
P(positive | confounds) instead, hardest positives first, with the caliper in
standard deviations of that logit (0.2 is the usual width). On the `lalonde`
data it fits the same logits as R's MatchIt (to 1e-7), picks the same 184 pairs
up to exact ties, and reports the same post-match SMDs to six decimals;
`tests/data/` holds MatchIt's output and `tests/test_balance.py` compares
against it. `standardized_mean_differences`, `balance_report` and
`propensity_logit` are usable on their own.

## Verdicts

| verdict           | meaning                                                                             |
| ----------------- | ----------------------------------------------------------------------------------- |
| `confound-robust` | enough of the anchor's above-chance signal survived, and the interval clears chance |
| `confound-driven` | the interval straddles chance — matching removed the signal                         |
| `partial`         | signal is real but diminished                                                       |
| `inconclusive`    | the interval is too wide to support any of the above (opt in via `max_ci_width`)    |
| `inverted`        | the interval clears chance from BELOW — the control reversed the ranking            |

`inconclusive` has no counterpart in the source this came from, which folded
those cases into `confound-driven` — reporting an absent measurement as a
finding.

`inverted` likewise had no counterpart: an interval entirely below chance
_excludes_ chance, so it fell through to `partial` — reporting a control that
reliably reversed the ranking as a diminished-but-real signal.

## Two more controls

**Incremental validity** — does the new feature add anything _beyond_ the
confounds? Fit confound-only and confound+feature models, evaluate both on the
same held-out rows, and put a **paired** bootstrap interval on the AUROC
difference.

Paired is the point. Two independent intervals on two AUROCs overlap far more
often than the interval on their difference excludes zero, because the pair is
computed on the same resampled rows and the shared variance cancels. Comparing
two separately-reported AUROCs by eye is exactly the error this prevents.

```python
from confound_controls import incremental_validity

res = incremental_validity(conf_train, conf_test, feat_train, feat_test, y_train, y_test)
res.verdict  # adds-signal | harms | no-added-signal | underpowered
```

`feat_train` must be out-of-fold. Fitting the feature on the rows the confound
model is evaluated against leaks the label into the augmented model and
manufactures the very lift the control is testing for.

The source's verdict was two-valued (`"FM-adds-signal" if lo > 0 else
"no-added-signal"`), so an interval straddling zero, an interval too wide to
say anything, and an interval lying entirely _below_ zero all reported the
same. The last of those is an augmented model that is reliably **worse** — a
finding, not an absence of one.

**Ablation** — destroy the structure of interest, keep the nuisance properties,
score with the same model, and see what survives.

```python
from confound_controls import assert_ablation_changed_input, ablation_control

assert_ablation_changed_input(real_inputs, ablated_inputs)  # do this first
res = ablation_control(y, p_real, p_ablated)
res.verdict  # structure-dependent | structure-independent | partial | inconclusive
```

That first call is not optional politeness. **A broken ablation leaves the
input unchanged, so the scores match exactly, the delta is zero, and the
control reports the model as surviving** — the strongest possible result,
obtained by doing nothing. No score-level statistic can catch it, which is why
the check belongs on the input. The source had this instinct
(`assert not np.allclose(shuffled, real)`); here it is a first-class function
that also handles non-float inputs like sequences, where `np.allclose` raises.

This is the same shape as the vacuous matched control above: an inert control
returns the original answer, which reads as the hypothesis surviving.

## What changed from the source

Extracted from an internal matched-negatives script, present byte-identically
in a second internal project. The code had been copied between two independent
repositories, which is what marked it as worth extracting:
duplication across repos is revealed reuse demand rather than a proxy for it.

The matching and the statistics are unchanged. What changed:

- **`ANCHOR = 0.629` was a module-level constant** — one study's own baseline
  AUROC compiled into the library, against which every imported use silently
  scored. Now a required argument, because there is no defensible default for
  "what does good look like in your problem".
- **Input paths were derived from `__file__`** and the confound columns
  (`gc`, `log_basemean`, `og_size`, `n_frac`, 64 `kmer3_*`) were hardcoded, so
  the module ran against one dataset in one repo layout. Both are arguments.
- **Incomplete matches were silent.** When the pool ran out the original
  returned fewer negatives than positives and said nothing — and the drop is
  not random, it hits the positives in the densest region of confound space.
  `MatchResult.complete` now answers this, and it is enforced by default.
- **The pass flag counted a hardcoded subset** of the study's own confound
  names, so a confound added to the spec was scored, printed, and then left out
  of the decision it existed to inform. `battery_passes` counts every result.
- **`recovery()` divided by `anchor - 0.5` without checking it.** An anchor at
  or below chance makes that ratio a non-quantity rather than a large number;
  it is now refused.
- **The bootstrap silently skipped single-class resamples.** With few positives
  most resamples are skipped and the interval comes from a handful of
  replicates. It now raises rather than returning a confident-looking number.

## Sequence ablation

For sequence models, the ablation is a knockout: scramble a span, keep its
dinucleotide composition, re-score, and ask whether the drop is real and larger
than a length-matched control knockout elsewhere.

```python
from confound_controls import (
    knockout_span,
    sample_control_span,
    grouped_delta_ci,
    confirm_knockout,
)

ko = knockout_span(seq, start, end, seed=1).require_changed()
ctl = sample_control_span(len(seq), end - start, motif_spans, rng).require_disjoint()

pooled = grouped_delta_ci(real_scores, ko_scores, family_ids)
control = grouped_delta_ci(real_scores, ctl_scores, family_ids)
confirm_knockout(pooled, control)  # confirmed | not-confirmed
```

**`grouped_delta_ci` resamples whole groups, not rows.** Promoters from one
paralog family are not independent observations; a row-level bootstrap treats
96 correlated rows as 96 independent ones and reports an interval far narrower
than the data supports. Measured on the bundled fixture: grouped 0.189 wide vs
row-level 0.071.

**Three silent fallbacks, all pointing the same way.** Each returned a
plausible value on failure, and each failure makes the ablation weaker or
absent — which reads as the model _surviving_ it:

- the dinucleotide shuffle ended `return seq` on failure, so an un-shuffleable
  span came back **unchanged** — a knockout that knocked nothing out. The
  comment called this "rare"; measured on a real 29-mer, knocking out `[8:20]`
  leaves it unchanged for **2 of the first 8 seeds**.
- the control-span sampler tried 50 times to avoid the real motif spans, then
  returned an overlapping start regardless — a partial real knockout posing as
  the negative control it is compared against.
- the cluster bootstrap never checked it had enough distinct groups to resample.

All three now report. `ShuffleResult.changed`, `ControlSpan.disjoint`, and a
`min_groups` floor, each with a `require_*` that raises.

The Altschul-Erikson shuffle itself is ported unchanged — **including its
determinism fix**. Using `set(graph)` for the vertex list made `rng.choice`
consume the PRNG in a different order per process (set iteration depends on the
randomized hash of string keys), so the same `(seq, seed)` produced different
shuffles run to run. A same-process test cannot catch that; the suite runs four
subprocesses and compares.

## Not extracted

`aim3_confound_table.py` builds the confound table itself from FASTA, DE
results and orthogroup sizes. That is study data-prep rather than a control,
and its sequence-composition code would duplicate what `seqbench` already
does — so it stays where it is.

## The source repos are untouched

Both source projects are active. This is
extract-and-leave-source-alone: nothing was refactored there, and neither repo
imports this package.

## Install

```bash
pip install confound-controls
```

From a checkout, for development:

```bash
pip install -e ".[dev]"
pytest -q
```

Requires numpy, pandas, scipy, scikit-learn.
