RFC: stop cleaning the design matrix inside fit(); warn on rank deficiency

[ljchang]

Opening this as an RFC — the diff is one concrete proposal, but the design question is worth discussing before merging. Alternatives and open questions are at the bottom; I'm happy to dial this back to a narrower change.

The problem

fit(model='glm') runs DesignMatrix.clean() on X by default, dropping any column that correlates >= 0.95 with an earlier one. I hit this porting the dartbrains course to 0.6, where it surfaced as a baffling error:

ValueError: Contrast vector length (48) must match number of regressors (46)

Digging in turned up three separable issues.

1. It is a correlation heuristic, not a rank test. On a real first-level design from the course:

cosine_1 and poly_1 correlated at 0.99 >= threshold 0.95. Dropping poly_1
cosine_2 and poly_2 correlated at 0.96 >= threshold 0.95. Dropping poly_2
full: 48 -> cleaned: 46
actual matrix rank: 48 of 48 columns

The design was full rank and perfectly estimable. The docstring justifies cleaning in terms of rank deficiency, but that is not the condition being tested.

2. It is order-dependent. clean() keeps the first column of a correlated pair and drops the second, so which regressor survives depends on the order the design happened to be built in:

base.add_dct_basis(duration=128).add_poly(order=2)   # -> drops poly_1, poly_2
base.add_poly(order=2).add_dct_basis(duration=128)   # -> drops cosine_1, cosine_2

Same regressors, same data, two different fitted models. With design_clean_exclude_confounds=False as the default, regressors of interest are eligible too — not just nuisance.

3. It is silent. The only reporting is behind progress_bar, and fit() passes progress_bar=None. So the estimated model differs from the specified one with no indication.

The other half: the genuinely dangerous case was silent too

While testing design_clean=False I found the reverse gap. A truly singular design fits without complaint:

dm = DesignMatrix({'a': a, 'a_dup': a.copy(), 'c': np.ones(n)})   # rank 2 of 3
bd.fit(model='glm', X=dm, design_clean=False)
#   fit succeeded -- no error
#   betas finite: True
#   beta[a] == beta[a_dup]: True      <- effect split evenly, not unique
#   rank/singular warnings emitted: NONE

nilearn falls back to a pseudo-inverse and splits the effect across the dependent columns. The betas come back finite and plausible-looking, and any contrast touching that subspace is uninterpretable. So the status quo is unsafe in both directions: it intervenes where it shouldn't, and stays quiet where it should speak up.

What this PR does

• Removes design_clean, design_clean_thresh, design_clean_exclude_confounds, design_clean_fill_na from fit(). It estimates exactly the design it is given. DesignMatrix.clean() already exists as the explicit, discoverable way to drop columns — having an explicit method and a silent implicit call with different defaults is the confusing part.
• Adds a rank check that warns with the rank, column count, column names, and the next step:

UserWarning: Design matrix is rank deficient: rank 2 of 3 columns
(Intercept, condA, condA_dup). At least 1 column(s) are linear combinations of
the others, so the betas are not uniquely determined and contrasts involving
them are not interpretable. Inspect collinearity with `DesignMatrix.vif()` and
drop redundant regressors explicitly with `DesignMatrix.clean()` before fitting.

Net effect: fit() stops making modeling decisions for you, and starts telling you when your model is actually broken.

Open questions — input wanted

1. Warn or raise on rank deficiency? I chose warn, because over-parameterized designs can still have estimable contrasts, and raising would break pipelines currently relying (silently) on the pseudo-inverse. Raising is a one-line change if you'd rather fail hard. There's also a middle option: raise only when the deficiency is exact duplication, warn on near-singularity.

2. Full removal, or keep the flag with default=False? This PR removes the kwargs entirely. Keeping design_clean=False as an opt-in is less disruptive and preserves a convenience path — but it also keeps two ways to do the same thing.

3. Rank tolerance. I use np.linalg.matrix_rank defaults. fMRI designs with DCT bases can be numerically borderline; a condition-number check might be more informative than a hard rank threshold. I did not want to invent a threshold without input.

4. Should BrainCollection.fit() warn too? It doesn't currently take these kwargs and routes through a different path — worth confirming the per-subject designs get the same diagnostic, since that's where an unnoticed rank problem would be most costly.

5. Is anything relying on the implicit cleaning? Only one test did (test_datasets.py passed design_clean=False to defeat it); it now just calls fit(). If real pipelines depend on it, the deprecation path matters more than I've assumed.

Testing

Red-first. The four new tests fail on master:

• test_fit_estimates_every_column_given — a full-rank design with two regressors at r=0.99 keeps all three columns
• test_design_clean_kwargs_are_rejected — all four removed kwargs raise TypeError
• test_rank_deficient_design_warns / ..._names_the_columns — singular design warns, actionably
• test_full_rank_design_does_not_warn — no spurious warnings

Suites: nltools/tests/data/braindata, nltools/tests/models, nltools/tests/data/collection — 742 passed. nltools/tests/support — 88 passed. uv run poe lint clean.

test_ridge.py::TestRidgeCore::test_vs_sklearn fails, but it also fails on unmodified origin/master — pre-existing and unrelated.

Migration guide updated. API docs not regenerated: docs/api has pre-existing drift from #465 and poe docs-generate is a separate step, so including it would have mixed ~390 unrelated lines into this diff.

🤖 Generated with Claude Code

https://claude.ai/code/session_01YVDnsKutg6Yqv99dYmUeuP

[ljchang]

Update: the warning now recommends regularization, not deletion

Review feedback (@ljchang) — regularization should be encouraged over clean() when a design is rank deficient. That's right, and the first version of this PR pointed users at the weaker tool. Fixed in 10da6d2.

The argument in one line: ridge has a unique solution even when X'X is singular, because (X'X + αI) is always invertible — and that solution doesn't depend on column order. Deletion resolves the deficiency by discarding information and assigning the shared variance to whichever column sorted first.

Measured on a design with two regressors at r=0.99:

	behavior
fit(model='ridge')	swap the collinear columns → weights swap exactly (equal to float32 solver precision, ~1e-6)
clean()	order [a, b, c] keeps ['a', 'c']; order [b, a, c] keeps ['b', 'c'] — a different model

That property is now pinned by test_ridge_is_order_invariant_where_clean_is_not, and the warning-content test asserts the guidance names ridge and vif so it can't silently regress.

New warning text:

UserWarning: Design matrix is rank deficient: rank 2 of 3 columns
(Intercept, condA, condA_dup). At least 1 column(s) are linear combinations of
the others, so the OLS betas are not uniquely determined and contrasts involving
them are not interpretable. Prefer regularization: `fit(model='ridge')` keeps
every regressor and shrinks them, giving a unique solution that does not depend
on column order. Inspect the collinearity first with `DesignMatrix.vif()`.
Dropping columns with `DesignMatrix.clean()` also removes the deficiency, but it
discards information and which column survives depends on the order the design
was built in.

One caveat I made explicit in the migration guide rather than overselling ridge: regularization fixes the estimation problem, not the identifiability one. For exactly collinear regressors, no method can separate their individual contributions — that information isn't in the data. Ridge buys a stable, reproducible answer in place of an arbitrary one; it doesn't recover something that was never measured. Worth saying out loud so nobody reads the warning as "add ridge and the collinearity is handled."

This also sharpens open question 1 above. If regularization is the recommended path out of rank deficiency, warn-don't-raise looks more clearly correct: the warning now names a concrete fix that keeps the full model, rather than telling people to delete regressors.

Suites re-run: nltools/tests/data/braindata + nltools/tests/support — 605 passed. poe lint clean.

[ljchang]

Added: find_spikes(clean=True) — stop generating duplicate spike regressors

Follow-up from review (@ljchang). Pushed in 2b1dbf4.

find_spikes() runs two independent detectors (per-TR global signal, and mean absolute frame-to-frame difference) and turns every detection into its own one-hot column. A single bad volume is routinely caught by both, so the same TR gets flagged twice — and two one-hot columns marking the same TR are bitwise identical:

find_spikes(img, global_spike_cutoff=0.8, diff_spike_cutoff=0.8)
  clean=False   cols= 23  unique= 16  rank= 16   full_rank=False
  clean=True    cols= 16  unique= 16  rank= 16   full_rank=True

On real localizer data at looser cutoffs, 22 of 23 spike columns collided with another. So nltools was manufacturing exactly the degeneracy this PR teaches fit() to warn about — from a function whose entire job is handing you nuisance regressors.

clean: bool = True is added to nltools.stats.find_spikes and plumbed through BrainData.find_spikes / find_spikes_data. When a TR is flagged more than once the global_spike column wins, so the tie-break is deterministic rather than insertion-ordered. clean=False restores one column per detection.

Why this defaults on when the design cleaning removed above did not

Worth stating, since the two look superficially similar:

	fit(design_clean=True) (removed)	find_spikes(clean=True) (added)
what it drops	regressors the user chose	columns the function itself just created
are they identical?	no — merely correlated at ≥ 0.95	yes — bitwise identical
is the choice determined?	no — depends on column order	yes — nothing distinguishes the duplicates
information lost	real, and reassigned to a survivor	none

Dropping a regressor that correlates 0.95 with another is a modeling decision and belongs to the caller. Refusing to emit the same regressor twice is just not producing garbage.

Tests: no duplicates by default; full rank by default; clean=False preserves every detection; global_spike wins the tie-break; distinct TRs never merged; confound marking survives.

Suites: nltools/tests/stats + nltools/tests/data/braindata — 505 passed, poe lint clean.

One unrelated pre-existing failure noticed

TestFindSpikes::test_find_spikes_brain_data fails on unmodified origin/master (verified by checking out master's copy of both files). When no spikes are found, the empty-frame branch does pl.DataFrame({"_no_spikes": [0]*n}).drop("_no_spikes"), and dropping the only column loses the height — so the returned DM is (0, 0) rather than (n_tr, 0) and dm.shape[0] == len(data) fails.

Polars can't express "n rows, 0 columns", so a real fix needs DesignMatrix to carry an explicit height. I left it alone rather than widen this PR — happy to file it separately.

[ljchang]

Also fixed: no-spike subjects no longer break the design build

Picked up the pre-existing failure I flagged rather than filing it separately (@ljchang). Pushed in 5f35067.

It turned out not to be cosmetic. Finding no spikes is a perfectly normal outcome, but it took down the entire first-level design for that subject:

spikes = bold.find_spikes(...)      # clean subject, nothing detected
task.append(spikes, axis=1)
# ValueError: All Design Matrices must have the same number of rows!

In a group loop, one well-behaved participant crashes the pipeline — and the error says nothing about spikes.

Cause

Polars derives a frame's height from its columns, so a frame with no columns always reports 0 rows. find_spikes built its empty result as:

pl.DataFrame({"_no_spikes": [0] * len(global_mn)}).drop("_no_spikes")

Dropping the only column discards the height, so the result was (0, 0) and then failed append()'s row check against the rest of the design.

Fix

A design matrix with no regressors still describes a specific number of timepoints, so that length has to be carried explicitly:

• DesignMatrix takes an optional n_rows, used only when the frame has no columns; shape and __len__ consult it
• find_spikes passes the TR count through, so the empty result reports (n_tr, 0)
• append() additionally drops regressor-less matrices before the row check, so this composes even for a column-less matrix built without an explicit height, and appending only empty matrices returns a copy instead of raising

Belt and braces on purpose: the n_rows path makes the object self-describing, the append guard makes the operation robust regardless.

Result

TestFindSpikes::test_find_spikes_brain_data — which fails on unmodified origin/master — now passes, along with new tests covering: the empty result reporting input length, appending as a no-op, an end-to-end clean-subject design build, and .is_empty still being True.

Suites: nltools/tests/stats + nltools/tests/data — 1047 passed. poe lint clean.

This PR now carries four related changes, all circling the same theme of designs being silently malformed:

1. fit() stops silently cleaning the design
2. fit() warns on rank deficiency, recommending regularization
3. find_spikes(clean=True) stops emitting duplicate regressors
4. a regressor-less design matrix keeps its row count

Happy to split (3) and (4) into their own PR if you'd rather review the fit() API change on its own.

[ljchang]

Split: the find_spikes fixes moved to #472

Per review (@ljchang), the two find_spikes commits are now their own PR — #472 — so the fit() API change can be reviewed on its own.

This PR (#470) is now just:

1. fit() stops implicitly cleaning the design (removes the design_clean* kwargs)
2. fit() warns on rank deficiency, recommending regularization

#472 carries:

3. find_spikes(clean=True) — stop emitting duplicate spike regressors
4. a regressor-less design matrix keeps its row count

The branch was force-pushed back to 10da6d2; the two commits were cherry-picked onto master for #472 and verified to stand alone there (1046 tests, lint clean, no design_clean in the diff).

The migration-guide entry moved with them, and its cross-reference back to this PR's fit() section was dropped so #472's docs make sense independently. The two are otherwise unrelated and can merge in either order.

Re-verified here after the trim: nltools/tests/data/braindata + nltools/tests/support — 605 passed, poe lint clean.

One consequence worth noting: TestFindSpikes::test_find_spikes_brain_data fails on master and is fixed in #472, not here — so it will still fail against this branch alone.

[ljchang]

@ejolly. I feel pretty strongly that we should never have a default enabled that automatically drops regressors in our .fit(). I do not think the clean kwargs belong in the fit. it can be a useful function for a design matrix, but not for model fitting. I think this will create way more unintentional problems than protecting people from rare issues. The ordering issue is a problem where different regressors will be dropped based on the order in which they are run. not true for regularization.

Here is what i propose and have started the PRs for:

1. remove clean from fit kwargs entirely.
2. we should instead provide warnings if there is a degenerate design matrix so less technically minded users are aware (like most other stats software).
3. we should try to remove the issues in other ways. For example, I think a common one is duplicate spikes being detected in .findspikes(). I tried to add a deduplicator by default there.
4. I think documentation should encourage regularization rather than using .clean()

[ejolly]

Lets make a loud warning and provide possible solutions and suggestions:

• "Try using .clean() because your design matrix is ...."
• "Try using regularization..."

[ejolly]

Superseded by #483, which merges this branch as-is (all design_clean* kwargs removed; fit() estimates exactly the design it's given) and reworks the warning per the discussion here: diagnosis first with likely-involved columns named via pivoted QR (truncated), the p > n case now warns instead of being skipped, all three next steps offered (.vif(), .clean() with the order-dependence caveat, fit(model='ridge')), and a dedicated RankDeficientDesignWarning category. This branch's history is preserved in #483.