A Python pipeline for activation and connectome predictive modeling using elastic net regularization. Designed for SLURM high-performance computing environments, with a YAML configuration interface and a 3-stage job orchestration pattern (main analysis → permutation workers → aggregation).
The pipeline implements nested cross validation, elastic net regression (single- or multi-task) or classification (binary or multi-class) and optional feature dimensionality reduction. Two analysis modes are available to suit different sample sizes and inferential goals. Coefficient significance is based on bootstrap confidence intervals, probability of direction (pd), and Benjamini-Hochberg FDR correction. Block permutation testing quantifies the unique contribution of user-defined feature subsets.
conda env create -f environment.yaml
conda activate fmri-elastic-netRequirements (pinned to tested versions):
- Python 3.10
- numpy 2.2.x, pandas 2.3.x, scipy 1.15.x
- scikit-learn 1.7.x, hdbscan 0.8.x, joblib 1.5.x, pyyaml 6.0.x
cp config_template.yaml my_project/config.yaml
cp run_template.sh my_project/run.shAt minimum, set:
analysis_type:"regression"or"classification"analysis_mode:"predict"or"correlate"(see Analysis Modes below)covariate_method:"none","incorporate", or"pre_regress"feature_reduction_method:"none","cluster_pca","apriori", or"ica"paths.data_file: absolute path to input CSVpaths.output_dir: absolute path to output directorycv_params.n_outer_folds,cv_params.n_inner_folds: fold counts (or"loo")cv_params.n_random_search_iter: hyperparameter search iterations (no default; required)
python fmri-elastic-net.py --config /path/to/config.yamlEdit run.sh to set paths and resource parameters, then submit:
sh run.shThe orchestrator (run_fmri-elastic-net.sh) submits three dependent SLURM jobs:
| Stage | Job name | Description |
|---|---|---|
| 1 | EN_Main |
Nested CV, selection frequency, bootstrap, block permutation |
| 2 | EN_Worker (array) |
Permutation null distribution (parallelized across array jobs) |
| 3 | EN_Agg |
Aggregates permutation chunks, computes final p-value |
The input data file must be a CSV with one row per subject. Required columns:
- Subject ID: any string or integer identifier (specified by
data_cols.subject_id_col) - Outcome: numeric continuous (regression) or integer class labels (classification); specified by
data_cols.post_score_col - Brain features: numeric columns identified by a substring match (e.g., all columns named
brain_*); specified bydata_cols.brain_feature_substr - Covariates (optional): numeric columns; specified by
data_cols.covariate_cols
Rows with any missing values are removed by listwise deletion before analysis.
The pipeline supports four distinct outcome types, determined by analysis_type and the
shape of the outcome column(s):
| Outcome type | analysis_type |
Outcome format | sklearn estimator |
|---|---|---|---|
| Single-task regression | "regression" |
One continuous numeric column | ElasticNet |
| Multi-task regression | "regression" |
Multiple continuous numeric columns | MultiTaskElasticNet |
| Binary classification | "classification" |
One integer column, 2 unique labels | LogisticRegression(penalty='elasticnet', solver='saga') |
| Multi-class classification | "classification" |
One integer column, 3+ unique labels | LogisticRegression(penalty='elasticnet', solver='saga') with OVR decomposition |
For multi-class classification, data_cols.reference_class must be specified. This parameter designates the reference class label; all coefficients and confidence intervals are reported as class-vs-reference contrasts (K-1 contrasts for K classes). The pipeline halts with an informative error if reference_class is absent for a multi-class outcome. For binary classification and regression, reference_class is ignored.
Multi-task regression (MultiTaskElasticNet) has two important constraints:
- Shared sparsity: all tasks share the same feature sparsity. Features are jointly selected (non-zero) or jointly excluded (zero) for all tasks simultaneously (assumes a shared 'system' for predictive features). If different tasks are driven by different feature subsets, consider running separate single-task analyses instead.
- Sample weights via sqrt(w) transformation:
sample_weight_colis supported forMultiTaskElasticNetvia aWeightTransformerpipeline step that scales both X and Y by sqrt(w_i) before fitting. This is algebraically equivalent to weighted least squares in the loss term. The L1+L2 penalty is invariant to this transformation; alpha is re-tuned on the transformed data.
Multi-task regression is triggered automatically when data_cols.post_score_col refers
to multiple columns in the data file.
Choosing the right analysis mode is the most consequential configuration decision. It controls both the L1 ratio search space and the implicit regularization.
- L1 ratio search space: 0.01–0.99 (full elastic net; data-driven regularization balance)
- Model performance evaluated via external validity (nested CV R² or AUC)
- Emphasis on generalization
- Best use case: large samples (N ~ 1,000s) where generalization to unseen data is the primary scientific question; nested CV selects the optimal L1/L2 balance for the data, including Ridge-dominant solutions when features are multicollinear or signal is diffuse
- L1 ratio search space: 0.001–0.2 (dense, shrinkage-focused solutions)
- Emphasis on stable coefficient estimation; Ridge-dominant regularization reduces variance at the cost of sparsity, suitable for multicollinear feature sets
- Model performance is still evaluated via nested CV on held-out folds, but the primary goal is reliable, non-zero feature attribution (internal validity)
- Best use case: small-to-medium samples (N ~ 100s) where overfitting to noise or P>>N is a concern; scenarios where interpretability of all features is desired rather than sparse selection
cv_params.cv_strategy controls how outer and inner CV folds are constructed:
| Strategy | Splitter (classification) | Splitter (regression) | Best use case |
|---|---|---|---|
stratified (default) |
StratifiedKFold |
KFold |
Independent observations; class balance preserved across folds |
uniform |
KFold |
KFold |
Independent observations where stratification is unnecessary |
group |
StratifiedGroupKFold |
GroupKFold |
Clustered or longitudinal data (e.g., repeated measures per subject) where observations sharing a group label must stay in the same fold |
group requires data_cols.group_column (e.g., a subject ID for repeated-measures
data). All observations sharing a group label are assigned to the same outer and
inner fold, preventing leakage across repeated measures. StratifiedGroupKFold
balances class proportions on a best-effort basis and may not achieve perfect
stratification with small or imbalanced groups; the pipeline logs a warning when a
fold's test set is missing a class or has fewer than 2 instances of a class.
Group assignment extends to every stage: permutation testing uses size-stratified
whole-group block permutation (Winkler et al., 2015) instead of global label
shuffling, selection frequency subsamples roughly half the groups (not individual
observations), and bootstrap importance resamples groups with replacement (cluster
bootstrap). fold_group_mapping.csv records each group's fold assignment.
The effective sample size relative to group count (ESS / N_groups) is logged at load time; a warning is triggered below 0.5, indicating that sample weight heterogeneity may be inflating apparent group-level power.
cv_params.use_calibration: true enables fold-local post-prediction calibration.
Parameters are fit on out-of-fold (OOF) predictions collected from an inner CV loop
within each outer training fold — never on the outer test fold — then applied to
that fold's held-out predictions. The method is selected automatically from
analysis_type:
analysis_type |
Outcome | Method | Transformation |
|---|---|---|---|
regression |
any | Linear | y_cal = a * y_pred + b (least-squares fit) |
classification |
binary | Platt scaling | p_cal = sigmoid(a * p_raw + b) (NLL fit) |
classification |
multi-class (3+) | Temperature scaling | p_cal = softmax(logits / T) (NLL fit) |
Calibration parameters are fit per outer fold and stored in
calibration_parameters.csv. Calibration affects only reported predictions and
probabilities (predict_ensemble, model_performance.csv); it has no effect on
model coefficients, so Tier 1/Tier 2 inference and selection frequency are
unaffected by use_calibration.
| Method | Description | Best use case |
|---|---|---|
none |
No covariates included | No nuisance variables to control |
incorporate |
Covariates entered as features with tunable penalty weight | Covariates are substantively interesting predictors alongside brain features |
pre_regress |
Outcome residualized on covariates fold-locally before prediction | Covariates are pure nuisance variables; their unique contribution should be removed |
When incorporate is used, a covariate_penalty_weight hyperparameter is searched
via a loguniform distribution over [model_params.covariate_penalty_weight_min, model_params.covariate_penalty_weight_max] to control regularization applied to
only covariate columns. Because this adds a third continuous hyperparameter to the
search space, cv_params.n_random_search_iter should be increased.
An optional data_cols.sample_weight_col column of non-negative observation weights
is supported for all estimators:
- ElasticNet / LogisticRegression: weights passed natively to sklearn's
sample_weightargument; sklearn normalizes internally. - MultiTaskElasticNet: weights applied via sqrt(w_i) pre-transformation of X and Y
(algebraically equivalent to weighted least squares in the loss term; alpha is
re-tuned on the transformed data via
WeightTransformer).
The pipeline normalizes weights to sum to N at load time (Hajek estimator convention; Lumley, 2010). Only relative weights matter; the absolute scale of the raw weight column has no effect on model behavior. The effective sample size ESS = (Σwᵢ)² / Σwᵢ² is logged at runtime. ESS/N < 0.5 triggers a warning: extreme weight heterogeneity reduces effective sample size below half the nominal N, which may compromise regularization path stability and bootstrap coverage. No built-in weight trimming is applied — users should examine their weight distribution and consider trimming extreme weights as a sensitivity analysis before running the full pipeline.
An optional moderating variable can be specified via data_cols.moderator_col and
data_cols.moderator_type to test whether the association between brain features and
the outcome varies as a function of a subject-level characteristic (e.g., treatment
group, age, symptom severity).
| Type | Coding | Columns produced |
|---|---|---|
continuous |
Mean-centered on training-split mean (fold-local) | 1 column |
nominal |
Deviation (effect) coding relative to the last sorted level (fold-local frequencies) | K-1 columns (K = number of levels) |
Interactions are constructed post-dimensionality-reduction (after the reducer
transforms brain features into components or passes them through for none). For each
fold or resampled iteration, the pipeline:
- Codes the moderator fold-locally via
_code_moderator. - Constructs element-wise products of each brain feature column with each moderator
column via
_construct_interactions. - Prepends the moderator main effect columns to the feature matrix.
- Appends the interaction columns after the brain features.
The resulting feature matrix layout is:
[covariates | moderator main effect(s) | brain features | brain x moderator interactions]
The moderator main effect columns are protected from regularization via
ModeratorScaler, which applies a fixed 1000x amplification (equivalent to
penalty_weight = 0.001). This ensures the elastic net's penalty effectively ignores
the moderator main effect, analogous to CovariateScaler for covariates but with a
fixed (non-tuned) scale factor. The main effect is retained to satisfy the heredity
principle (Bien, Taylor, and Tibshirani, 2013): interaction terms should only enter the
model when the corresponding main effects are present.
Interaction coefficients are included in all downstream inference stages:
- Tier 1 (fold-level t-test): per-feature t-test applied to both main and interaction coefficients. For K>2 nominal moderators, Hotelling's T-squared provides an omnibus test across the K-1 contrast coefficients per brain feature.
- Tier 2 (bootstrap CIs): Partial Ridge refit preserves both main and interaction coefficients. For classification, the logistic adaptation uses sqrt(n) column scaling (see Known Limitations).
- Selection frequency and bootstrap importance: both stages construct interactions
per-iteration using full-sample levels (
levels_override) for nominal moderators, ensuring consistent coding dimensions across resampled iterations.
Main-effect and interaction-effect visualization data are written to separate CSVs
(report_{level}_plotting.csv and report_{level}_interaction_plotting.csv,
respectively). Interaction visualization includes a moderator_value column recording
each subject's moderator value for downstream plotting. Interaction visualization is
skipped for K>2 nominal moderators (partial-dependence decomposition is not well-defined
for multi-contrast interactions) and for apriori reduction (cluster-level interaction
visualization is not produced). Interaction-effect partial associations reflect the full conditional relationship (brain main effect + moderator main effect + interaction effect), computed via out-of-fold (OOF) ensemble linear predictions: each subject's contribution uses the pipeline (reducer, scaler, model coefficients) from the fold in which that subject was held out, rather than a single representative fold's model.
Multi-task/multi-class interaction Tier 1, Tier 2, and selection frequency reporting is not yet implemented. The pipeline currently produces correct coefficients for these configurations but does not produce interaction-specific output files.
All reduction is applied fold-locally inside the CV loop to prevent data leakage.
| Method | Description | Best use case |
|---|---|---|
none |
Raw features passed directly to model | Small–medium feature sets; no assumed structure |
cluster_pca |
HDBSCAN clustering + 1-component PCA per cluster (fit inside CV) | Questionnaire items, genomic data, or any features expected to form discrete non-overlapping groups |
apriori |
Externally-defined cluster map + 1-component PCA per cluster (fit inside CV) | Pre-defined brain networks (e.g., atlas-based parcellation); network structure is theoretically motivated |
ica |
FastICA decomposition (fit inside CV); back-projection via activation patterns (Haufe et al., 2014) | Brain activation or connectivity data where regions participate in multiple overlapping networks |
Note on leakage prevention: Reduction is strictly fold-local throughout: in the nested CV stage, reduction is fit on training data only. In the descriptive reporting stages (selection frequency, bootstrap importance), each iteration fits a fresh reducer clone on the resampled/subsampled data, preserving the conditional bootstrap framework.
Label permutation test: the full nested CV is repeated n_permutations times with
shuffled outcome labels. P-value uses Laplace correction:
(count(null ≥ observed) + 1) / (n_permutations + 1).
Bootstrap confidence intervals (conditional bootstrap: Efron & Tibshirani, 1993): each iteration fits a fresh reducer clone on resampled brain features, fits the model using fold-specific hyperparameters (fixed from each fold's inner-CV tuning), and back-projects coefficients to the original feature space for aggregation. Tuning variance from the nested CV is therefore propagated into the bootstrap CIs. This ensures meaningful CI and pd computation across iterations with different reduced spaces.
For classification, the Partial Ridge method adapts Liu et al. (2020) from linear to logistic regression via differential L2 penalization (selected features scaled by sqrt(n), C=1). This is a project-specific extension; see Known Limitations.
is_significant: primary criterion — CI does not cross zerois_significant_fdr: survives Benjamini-Hochberg FDR correction at q = 0.05pd: probability of direction; p-value approximationp = 2*(1 - pd)assumes a continuous coefficient distribution. For sparse features with high L1 regularization, zero-inflated bootstrap distributions causepdto be near 0.5; useis_significantas the primary criterion in that case.
Repeated 50% subsampling (n_fold_bootstraps per fold). Per feature: proportion of
iterations with a non-zero coefficient. Descriptive only — no significance threshold.
For each user-defined feature block: only that block's columns are row-permuted and the full nested CV is rerun. Quantifies the unique predictive contribution of a feature subset (e.g., activation from a specific brain network) beyond the remaining features. Requires the user to label feature columns with a "block" identifier (e.g. if we want to know the unique predictive contribution of fronto-limbic "FL" circuits relative to the rest of the brain, block permutation is applied to features that contain the string "FL" in its column header).
All output files are written to paths.output_dir:
| File | Description |
|---|---|
nested_cv_scores.csv |
Observed model performance (R² or AUC) |
model_performance.csv |
Comprehensive evaluation metrics (regression: RMSE, MAE, R², Pearson r; classification: AUC-ROC, Log-Loss, Sensitivity, Specificity, Balanced Accuracy) |
model_performance_per_fold.csv |
Per-fold performance metric (R2 for regression, AUC_ROC for classification) with fold index and held-out sample size |
model_performance_fold_summary.csv |
Summary statistics (mean, SD, min, max) of per-fold performance across all K folds |
confusion_matrix.csv |
Confusion matrix (multi-class classification only) |
permutation_null_distribution_{metric}.csv |
Null distribution from label permutation |
permutation_result.csv |
Observed score, p-value, n_permutations (aggregate mode) |
report_selection_frequency.csv |
Subsampling-based selection frequency per feature/component |
report_feature_importance.csv |
Bootstrap CIs, pd, FDR flags per feature/component (all feature_reduction_method values) |
report_cluster_importance.csv |
Bootstrap CIs, pd, FDR flags at cluster level (apriori) |
report_block_permutation.csv |
Block-specific observed score and p-value |
report_fold_ensemble_importance.csv |
Tier 1 inference: fold-wise t-test mean/SD/CV, t-statistic, p-value, CI, and significance flags per feature |
report_fold_diagnostics.csv |
Per-fold hyperparameter records (alpha/C, l1_ratio, penalty_weight) from nested CV |
report_fold_params_summary.csv |
Mean ± SD, min, max across K folds for each hyperparameter |
report_fold_bootstrap_ci.csv |
Tier 2 inference: pooled fold-wise bootstrap percentile CIs and significance flags per feature |
fold_group_mapping.csv |
Group-to-fold assignment (when cv_params.cv_strategy: "group") |
calibration_parameters.csv |
Per-fold calibration parameters (when cv_params.use_calibration: true) |
cv_predictions.csv |
Per-observation out-of-fold predictions (subject ID, true value, predicted value, fold, calibration status) |
cluster_loadings.csv |
PCA loadings per cluster (cluster_pca or apriori) |
cluster_loadings_fold_{n}.csv |
Per-fold PCA loadings for transparency |
ica_mixing_matrix.csv |
ICA mixing matrix A (P × K), activation pattern basis |
ica_mixing_matrix_fold_{n}.csv |
Per-fold ICA mixing matrix for transparency |
report_interaction_importance.csv |
Bootstrap CIs, pd, FDR flags per interaction term (when moderator is configured) |
report_{level}_plotting.csv |
Subject-level feature vs. outcome data for visualization (main effects) |
report_{level}_interaction_plotting.csv |
Subject-level interaction partial associations with moderator_value column (when moderator is configured; K>2 nominal excluded) |
bootstrap_coef_distribution.npz |
Full bootstrap coefficient array (when save_distributions: true) |
block_perm_null_{label}.csv |
Block-specific permutation null scores (when save_distributions: true) |
pipeline.log |
Full logging output with timing and diagnostics |
For multi-task regression or multi-class classification, per-task/per-class output
files are written to output_dir/task_{label}/ subdirectories, with aggregate
summaries written to the top-level output_dir.
See config_template.yaml for all parameters with inline comments and defaults.
See INPUT_SPECIFICATION.md for the complete exhaustive specification including
parameter types, ranges, constraints, output schemas, and known edge cases.
- Bootstrap importance uses fold-specific hyperparameters from each fold's inner-CV tuning (not full-dataset hyperparameters). Tuning variance is propagated into Tier 2 CIs. Each bootstrap iteration re-fits a fresh reducer clone (per-iteration re-reduction), preserving the conditional bootstrap framework across iterations with different reduced spaces.
- Tier 2 bootstrap CIs are computed from a pooled mixture of iterations with different fold-specific hyperparameter configurations. Percentile CIs from this mixture may have sub-nominal coverage for threshold-adjacent features; they are best interpreted as sensitivity diagnostics (Efron & Tibshirani, 1993, Ch. 13).
std_coef_meanreports the fully standardized coefficient (dimensionless: SDs of Y per 1 SD of X). For regression, the divisor is per-fold SD(Y). For classification, the divisor is per-fold SD(Y*) using the latent variable approach (Long, 1997; Menard, 2004, 2011): SD(Y*) = sqrt(Var(cross-validated logits) + pi^2/3). Bothraw_coef_mean(change in Y per unit change in X) andstd_coef_mean(fully standardized) are reported in all output files.raw_coef_meanfor reduction methods (cluster_pca, apriori, ica) is approximate: the back-projected coefficient divided by original-feature SD is not equivalent to a standardized beta from direct regression on original features.std_coef_meanandpdare the primary inferential quantities.- Selection frequency magnitudes may be elevated because hyperparameters are fixed from full-N tuning while each subsample uses N/2. Relative ordering is preserved; no significance threshold is applied (purely descriptive).
- Tier 1 fold-ensemble p-values (report_fold_ensemble_importance.csv) treat K fold-level coefficient estimates as independent observations in a one-sample t-test. Because adjacent folds share overlapping training data, the naive variance estimator is downward biased (Bengio and Grandvalet, 2004), producing anti-conservative p-values whose Type I error exceeds the nominal alpha by an algorithm-dependent amount. Tier 1 is designed as a liberal sensitivity screen; Tier 2 bootstrap CIs (report_fold_bootstrap_ci.csv) provide the confirmatory inference and are not affected by this bias.
- For classification, the Partial Ridge bootstrap CI method (Liu et al., 2020) is adapted from its original linear-model formulation to logistic regression via column scaling (selected features scaled by sqrt(n) with fixed L2 penalty C=1). This adaptation is not prescribed by Liu et al. (2020) and should be considered a project-specific extension. Percentile bootstrap CI coverage for this configuration (fold-wise-pooled, Partial-Ridge-refitted, elastic net) has not been directly benchmarked in the literature; however, percentile CIs in regularized settings tend toward conservative overcoverage (wider intervals), which is favorable for the pipeline's zero-crossing thresholding use case.
- LOO cross-validation disables
n_inner_repeats(repeated CV is undefined for LOO). cv_strategy: "group"with classification:StratifiedGroupKFoldcannot guarantee perfect class balance across folds, particularly with small or imbalanced groups. A warning is logged when a fold's test set is missing a class or has fewer than 2 instances of a class.- Group permutation with fewer than 20 groups: the null distribution's resolution is bounded by the factorial of group counts within each size stratum and can be coarse. Interpret permutation p-values with caution in this setting.