diff --git a/README.md b/README.md index 1c603c8..d78295e 100644 --- a/README.md +++ b/README.md @@ -217,6 +217,9 @@ src/behavioral_decoding/ docs/ ├── literature.md the neuroforecasting canon, with verified DOIs ├── design.md each decision, its rejected alternative, and what would falsify it +├── estimators.md per-modality base learners, and why mixed-effects/Riemannian need their own path +├── deap.md DEAP loader: format traps, circularity, the market route +├── narps.md NARPS loader: format, the individual-vs-aggregate honesty point └── data_sources.md candidate datasets per modality, and the gap between them ``` diff --git a/configs/experiment_default.yaml b/configs/experiment_default.yaml index ef526c9..28f04bb 100644 --- a/configs/experiment_default.yaml +++ b/configs/experiment_default.yaml @@ -42,9 +42,12 @@ model: weight_floor: 0.0 base_learner: - fmri: logistic # ~5 ROI features; a regularised linear problem - eeg: logistic - face: logistic # high-dimensional embeddings; heavily regularised + # elasticnet (L1+L2 logistic) shares weight across correlated features and + # still drops dead ones; see docs/estimators.md. linear_svm is a supported + # face alternative (comparable, slower). + fmri: elasticnet # correlated ROI betas (NAcc_L/NAcc_R) + eeg: elasticnet # collinear band-power columns + face: elasticnet # high-dimensional embeddings; linear_svm also fine behavior: gradient_boosting # tabular, mixed monotone effects n_bags: diff --git a/configs/experiment_majority_vote.yaml b/configs/experiment_majority_vote.yaml index 1bf3d5f..1171c30 100644 --- a/configs/experiment_majority_vote.yaml +++ b/configs/experiment_majority_vote.yaml @@ -32,9 +32,9 @@ model: drop_below_chance: true weight_floor: 0.0 base_learner: - fmri: logistic - eeg: logistic - face: logistic + fmri: elasticnet + eeg: elasticnet + face: elasticnet behavior: gradient_boosting n_bags: fmri: 25 diff --git a/docs/estimators.md b/docs/estimators.md new file mode 100644 index 0000000..e4deb13 --- /dev/null +++ b/docs/estimators.md @@ -0,0 +1,125 @@ +# Per-modality estimators + +Which base learner each modality uses, why, and how to change it. Also an honest +account of two stronger models that were suggested and do **not** drop into this +framework's pipeline unchanged, and what it would take to add them properly. + +## The pipeline each estimator lives in + +Every modality model is a bag of identical pipelines: + +``` +BaggingClassifier + └── Pipeline(StandardScaler -> AdaptiveOverSampler(SMOTE) -> base estimator) +``` + +So a base estimator here always sees **standardised, class-balanced, flat +feature vectors**. Two consequences that decide what fits: + +- The estimator must expose `predict_proba`. The ensemble weights modalities on + out-of-fold probabilities (`evaluation/cv.out_of_fold_proba`), so a learner + without calibrated probabilities is unusable as-is. +- The features are a flat `(n_trials, n_features)` matrix, column-standardised. + Anything whose feature space is *not* a flat Euclidean vector (a covariance + matrix on its manifold, say) is distorted by per-column scaling and by SMOTE's + straight-line interpolation before the estimator ever sees it. + +## Defaults + +| modality | default | why | +|---|---|---| +| fMRI | `elasticnet` | ROI betas are correlated (NAcc_L/NAcc_R move together); L2 shares weight across them, L1 drops dead ROIs | +| EEG | `elasticnet` | band-power and ERP columns are collinear; regularised logistic is the right baseline | +| face | `elasticnet` | high-dimensional embeddings; `linear_svm` is an equally good, slower alternative | +| behaviour | `gradient_boosting` | low-dimensional, mixed, monotone economic/self-report features | + +Set them per experiment in the config: + +```yaml +model: + base_learner: + fmri: elasticnet # or logistic, linear_svm, svm, random_forest + face: linear_svm +``` + +## The learners + +- **`elasticnet`**: L1+L2 logistic (`saga` solver, `l1_ratio=0.5`). The default + for every dense neural/embedding block. The L2 term shares weight across + correlated features instead of arbitrarily picking one, which plain L1 (lasso) + does badly with collinear ROIs; the L1 term still zeroes out uninformative + features, which plain L2 never does. Needs standardised input, which the + pipeline provides. +- **`linear_svm`**: linear-kernel SVC with Platt-scaled probabilities. A strong + baseline for high-dimensional embeddings and often within noise of elastic-net + on the face block. Slower, because `probability=True` fits an internal CV for + calibration; that is why elastic-net is the face default and this is the + documented alternative. +- **`logistic`**: plain L2 logistic. Kept as a simple, fast reference. +- **`svm`**: RBF-kernel SVC. Non-linear; rarely the right call at these sample + sizes, but available. +- **`random_forest`**, **`gradient_boosting`**: tree ensembles for the tabular + behaviour block. + +## Two stronger models, and why they are not drop-ins + +Both were suggested and both are genuinely better *for the right setup*. Neither +fits the flat-feature, subject-grouped-CV pipeline without a dedicated path, and +bolting them in naively would produce numbers that look principled and are not. +Recorded here so the decision is deliberate, not forgotten. + +### Hierarchical / mixed-effects model (fMRI) + +A random-intercept-per-subject logistic model is the textbook way to handle +subject-to-subject differences, and within a single sample it usually beats a +pooled model. The catch here is the **cross-subject cross-validation**. This +framework groups folds by subject, so every test subject is unseen at fit time. +A per-subject random intercept has no estimate for a subject the model never +saw; the best it can do at test is fall back to the population mean, which is +what a pooled model already gives. So the partial-pooling benefit that makes +mixed models shine largely evaporates for *held-out-subject* generalisation, +which is the quantity this framework reports. + +Mixed effects would help if the evaluation were within-subject (predicting new +trials for subjects already seen). It is not, on purpose: individuating a person +from their own repeated trials is a much easier and less interesting claim than +generalising across people. Elastic-net with `class_weight="balanced"` plus +subject-grouped CV is the honest default here. + +If you still want it: fit `statsmodels` `BinomialBayesMixedGLM` (a random +intercept per subject) inside each training fold, and at test time drop the +random-effect term and predict from fixed effects only. Wrap it to expose +`predict_proba`, and report it next to elastic-net rather than replacing it. Do +not expect a cross-subject gain. + +### Riemannian methods (EEG) + +For EEG built from **channel covariance / connectivity** features, Riemannian +tangent-space classification is a real upgrade over band power: it respects the +geometry of the space of covariance matrices instead of treating their entries +as independent numbers. But two steps in this pipeline are wrong for covariance +features: + +- **StandardScaler** z-scores each covariance entry independently, which + destroys the positive-definite structure the Riemannian method depends on. +- **SMOTE** interpolates along straight lines between covariance vectors. The + geodesic between two covariance matrices is *not* a straight line in entry + space, so the synthetic minority points are off-manifold. + +A correct Riemannian EEG arm therefore needs its own pipeline: emit per-trial +covariance matrices (not band power), project to the tangent space at the +Riemannian mean of the *training* fold, and only then standardise / resample / +classify in that tangent space. That is a separate feature family and a separate +pipeline, best added as a distinct `riemann` path (with `pyriemann`, or a +scipy-`logm` tangent map for a dependency-free version) rather than a base +learner slotted into the existing flat pipeline. It is a good next step; it is +not a one-line default change, and pretending it were would corrupt the manifold +structure it exists to exploit. + +## What did change + +`elasticnet` and `linear_svm` were added to the learner factory and the +dense-block defaults moved from plain logistic to elastic-net. Everything else +about the pipeline (bagging, in-fold SMOTE, out-of-fold weighting, subject- +grouped CV) is unchanged, so these are honest swaps of the base estimator, not a +change to how anything is evaluated. diff --git a/src/behavioral_decoding/models/modality_models.py b/src/behavioral_decoding/models/modality_models.py index 6ea368d..035a158 100644 --- a/src/behavioral_decoding/models/modality_models.py +++ b/src/behavioral_decoding/models/modality_models.py @@ -37,10 +37,17 @@ logger = get_logger(__name__) +# Elastic-net logistic is the default for every dense-feature neural/embedding +# block: it shares weight across correlated features and still drops dead ones, +# which suits collinear ROI betas, collinear band-power columns, and +# high-dimensional face embeddings alike. `linear_svm` is a supported +# alternative for the face block (comparable accuracy, slower because of the +# probability calibration). Behaviour stays on gradient boosting for its mixed, +# monotone, low-dimensional features. See docs/estimators.md. DEFAULT_BASE_LEARNER: Dict[str, str] = { - FMRI: "logistic", - EEG: "logistic", - FACE: "logistic", + FMRI: "elasticnet", + EEG: "elasticnet", + FACE: "elasticnet", BEHAVIOR: "gradient_boosting", } @@ -82,6 +89,39 @@ def make_base_learner( params = {"n_estimators": 200, "max_depth": 3, "random_state": seed} params.update(kwargs) return GradientBoostingClassifier(**params) + if kind == "elasticnet": + # L1 + L2 logistic. The L2 part shares weight across correlated features + # (NAcc_L / NAcc_R move together; band-power columns are collinear), + # while the L1 part still drops dead ones. This is usually a better fit + # than plain L2 for the fMRI and EEG blocks, and than plain L1 for the + # high-dimensional face block. Needs the saga solver, and needs scaled + # input, which the pipeline's StandardScaler provides. + params = { + "penalty": "elasticnet", + "solver": "saga", + "l1_ratio": 0.5, + "C": 1.0, + "max_iter": 5000, + "tol": 1e-3, + "class_weight": class_weight, + "random_state": seed, + } + params.update(kwargs) + return LogisticRegression(**params) + if kind == "linear_svm": + # Linear-kernel SVM, a strong baseline for high-dimensional embeddings + # (the face block). probability=True adds Platt scaling via an internal + # CV so the ensemble can read predict_proba; that calibration costs time + # but the SMOTE step upstream keeps the classes balanced enough for it. + params = { + "kernel": "linear", + "C": 1.0, + "probability": True, + "class_weight": class_weight, + "random_state": seed, + } + params.update(kwargs) + return SVC(**params) if kind == "svm": params = { "kernel": "rbf", @@ -93,8 +133,8 @@ def make_base_learner( params.update(kwargs) return SVC(**params) raise ValueError( - f"unknown base learner {kind!r}; choose from logistic, random_forest, " - "gradient_boosting, svm" + f"unknown base learner {kind!r}; choose from logistic, elasticnet, " + "linear_svm, random_forest, gradient_boosting, svm" ) diff --git a/tests/test_estimators.py b/tests/test_estimators.py new file mode 100644 index 0000000..6f65830 --- /dev/null +++ b/tests/test_estimators.py @@ -0,0 +1,81 @@ +"""Base-learner factory: the estimators each modality can use. + +These check that every named learner instantiates with the right regularisation, +exposes ``predict_proba`` (the ensemble weights on out-of-fold probabilities, so +a learner without it is unusable), and survives a fit inside the leakage-safe +pipeline on imbalanced data. +""" + +from __future__ import annotations + +import numpy as np +import pytest +from sklearn.linear_model import LogisticRegression +from sklearn.svm import SVC + +from behavioral_decoding.balance.smote import make_balanced_pipeline +from behavioral_decoding.models.modality_models import build_modality_model, make_base_learner + +pytest.importorskip("imblearn") + + +def _imbalanced(seed: int = 0, n: int = 160, d: int = 8): + rng = np.random.default_rng(seed) + X = rng.normal(size=(n, d)) + y = np.array([0] * (n - n // 5) + [1] * (n // 5)) + # Plant a little separable signal so predict_proba is not degenerate. + X[y == 1] += 0.6 + return X, y + + +def test_elasticnet_is_l1_l2_logistic_with_saga(): + est = make_base_learner("elasticnet") + assert isinstance(est, LogisticRegression) + assert est.penalty == "elasticnet" + assert est.solver == "saga" + assert 0.0 < est.l1_ratio < 1.0 + + +def test_linear_svm_is_linear_kernel_with_probability(): + est = make_base_learner("linear_svm") + assert isinstance(est, SVC) + assert est.kernel == "linear" + assert est.probability is True + + +def test_unknown_learner_lists_the_valid_names(): + with pytest.raises(ValueError, match="elasticnet.*linear_svm"): + make_base_learner("transformer") + + +@pytest.mark.parametrize("kind", ["logistic", "elasticnet", "linear_svm", "random_forest"]) +def test_learner_fits_and_predicts_proba_in_pipeline(kind): + X, y = _imbalanced() + pipe = make_balanced_pipeline( + make_base_learner(kind, class_weight="balanced"), + sampler="smote", + k_neighbors=5, + random_state=0, + ) + pipe.fit(X, y) + proba = pipe.predict_proba(X) + assert proba.shape == (len(y), 2) + assert np.all((proba >= 0) & (proba <= 1)) + # Above chance on planted signal, so the column is informative. + assert proba[y == 1, 1].mean() > proba[y == 0, 1].mean() + + +def test_kwargs_override_defaults(): + est = make_base_learner("elasticnet", l1_ratio=0.2, C=0.5) + assert est.l1_ratio == 0.2 + assert est.C == 0.5 + + +@pytest.mark.parametrize("kind", ["elasticnet", "linear_svm"]) +def test_build_modality_model_accepts_new_learners(kind): + y = np.array([0] * 80 + [1] * 20) + model = build_modality_model("fmri", y=y, base_learner=kind, n_bags=5, seed=0) + assert model.n_estimators == 5 + assert model.bd_spec_["base_learner"] == kind + # class_weight balanced is passed through for both (neither is gradient boosting). + assert model.bd_spec_["class_weight"] == "balanced"