Advanced Regional MVPA Analysis in rMVPA

Introduction

Regional MVPA evaluates prediction performance within predefined brain regions. This vignette walks through the complete workflow using rMVPA:

We’ll generate synthetic data, define ROIs with a region mask, build the MVPA model with cross-validation, run the analysis, and examine the results. The implementation follows the approach in regional.R, mvpa_model.R, and dataset.R.

Data Generation and Preparation

We begin by generating a synthetic volumetric dataset using the gen_sample_dataset() function. This function creates a 4D array (with spatial dimensions and multiple observations), along with a binary mask and an associated design for cross-validation.

library(rMVPA)
library(neuroim2)
# Generate a synthetic dataset with dimensions 6x6x6, 80 observations, divided into 4 blocks
data_out <- rMVPA::gen_sample_dataset(D = c(6,6,6), nobs = 80, blocks = 4, nlevels = 2)
print(data_out)

## $dataset
## 
##  █▀▀ MVPA Dataset ▀▀█ 
## 
## ├─ Training Data 
## │  ├─ Dimensions:  6 × 6 × 6 × 80 observations 
## │  └─ Type:  DenseNeuroVec 
## ├─ Test Data 
## │  └─  None 
## └─ Mask Information 
##    ├─ Areas:  1 : 120 
##    └─ Active voxels/vertices:  120 
## 
## 
## $design
## 
##  █▀▀ MVPA Design ▀▀█ 
## 
## ├─ Training Data 
## │  ├─ Observations:  80 
## │  ├─ Response Type:  Factor
## │  ├─ Levels:  a, b 
## │  └─ Class Distribution:  a: 40, b: 40 
## ├─ Test Data 
## │  └─  None 
## └─ Structure 
##    ├─ Blocking:  Present
##    ├─ Number of Blocks:  4 
##    ├─ Mean Block Size:  20  (SD:  0 ) 
##    └─ Split Groups:  None

The returned list contains:

dataset: an MVPA dataset object with training data and a binary mask.
design: an MVPA design object specifying the response variable and block structure.

Creating a Region Mask

For regional analysis, we need to define ROIs. Here, we create a region mask by randomly assigning each active voxel in the binary mask a region label (from 1 to 3). This simulates a scenario where the brain is partitioned into three regions of interest.

# Extract the binary mask from the dataset
mask <- data_out$dataset$mask
nvox <- sum(mask)

# Create a regional mask: assign each voxel a random region number (1 to 3)
set.seed(123)  # for reproducibility
region_mask <- neuroim2::NeuroVol(sample(1:3, size = nvox, replace = TRUE), neuroim2::space(mask), indices = which(mask > 0))
table(region_mask)

## region_mask
##  0  1  2  3 
## 96 36 44 40

Setting Up the MVPA Model

Next, we create an MVPA model to evaluate a classification task. This involves:

Constructing an MVPA dataset using mvpa_dataset().
Specifying the design (including the block variable and response) via mvpa_design().
Defining the model with mvpa_model() using a chosen classifier and cross-validation strategy.

# Create MVPA dataset object from the generated training data and mask
 dset <- mvpa_dataset(data_out$dataset$train_data, mask = data_out$dataset$mask)

# Build cross-validation structure using block information from the design
cval <- blocked_cross_validation(data_out$design$block_var)

# Load a classification model; here we use "sda" (Shrinkage Discriminant Analysis)
mod <- load_model("sda")
tune_grid <- data.frame(lambda = 0.01, diagonal = FALSE)

# Create the MVPA model object
mvpa_mod <- mvpa_model(mod, dataset = dset, design = data_out$design, crossval = cval, tune_grid = tune_grid)
print(mvpa_mod)

## mvpa_model object. 
## model:  sda 
## model type:  classification 
## tune_reps:  15 
## tune_grid:  
##   lambda diagonal
## 1   0.01    FALSE
## 
##  █▀▀ Blocked Cross-Validation ▀▀█ 
## 
## ├─ Dataset Information 
## │  ├─ Observations:  80 
## │  └─ Number of Folds:  4 
## └─ Block Information 
##    ├─ Total Blocks:  4 
##    ├─ Mean Block Size:  20  (SD:  0 ) 
##    └─ Block Sizes:  1: 20, 2: 20, 3: 20, 4: 20 
## 
## 
##  █▀▀ MVPA Dataset ▀▀█ 
## 
## ├─ Training Data 
## │  ├─ Dimensions:  6 × 6 × 6 × 80 observations 
## │  └─ Type:  DenseNeuroVec 
## ├─ Test Data 
## │  └─  None 
## └─ Mask Information 
##    ├─ Areas:  1 : 120 
##    └─ Active voxels/vertices:  120 
## 
## 
##  █▀▀ MVPA Design ▀▀█ 
## 
## ├─ Training Data 
## │  ├─ Observations:  80 
## │  ├─ Response Type:  Factor
## │  ├─ Levels:  a, b 
## │  └─ Class Distribution:  a: 40, b: 40 
## ├─ Test Data 
## │  └─  None 
## └─ Structure 
##    ├─ Blocking:  Present
##    ├─ Number of Blocks:  4 
##    ├─ Mean Block Size:  20  (SD:  0 ) 
##    └─ Split Groups:  None

The mvpa_model() function, as defined in mvpa_model.R, packages all necessary parameters including cross-validation and performance computation.

Running the Regional Analysis

The regional analysis is executed with the run_regional() function, which:

Prepares the ROIs by extracting voxel indices from the region mask.
Iterates over the regions, applying the MVPA model to each ROI.
Compiles performance metrics and prediction tables.

# Run the regional analysis on the defined region mask
regional_results <- run_regional(mvpa_mod, region_mask)

The output is a regional_mvpa_result object containing:

performance_table: Cross-validated performance metrics per region.
prediction_table: Detailed predictions for each ROI and trial.
vol_results: Volumetric maps representing performance distributed across the brain.

Examining the Results

We can inspect the performance table to evaluate model accuracy in each region.

# Display performance metrics for each region
print(regional_results$performance_table)

## # A tibble: 3 × 3
##   roinum Accuracy     AUC
##    <int>    <dbl>   <dbl>
## 1      1    0.475 -0.0194
## 2      2    0.375 -0.114 
## 3      3    0.612  0.160

For a more detailed view, the prediction table shows trial-by-trial predictions:

# Display the prediction table
print(regional_results$prediction_table)

## # A tibble: 240 × 8
## # Rowwise: 
##    .rownum roinum observed pobserved predicted correct prob_a    prob_b
##      <int>  <int> <fct>        <dbl> <chr>     <lgl>    <dbl>     <dbl>
##  1       1      1 b           0.966  b         TRUE    0.0337 0.966    
##  2       2      1 a           0.663  a         TRUE    0.663  0.337    
##  3       3      1 a           0.599  a         TRUE    0.599  0.401    
##  4       4      1 b           0.952  b         TRUE    0.0481 0.952    
##  5       5      1 b           0.0166 a         FALSE   0.983  0.0166   
##  6       6      1 a           0.996  a         TRUE    0.996  0.00354  
##  7       7      1 b           0.507  b         TRUE    0.493  0.507    
##  8       8      1 a           0.0589 b         FALSE   0.0589 0.941    
##  9       9      1 a           1.00   a         TRUE    1.00   0.000411 
## 10      10      1 a           1.00   a         TRUE    1.00   0.0000256
## # ℹ 230 more rows

Volumetric results (vol_results) can be further visualized with neuroimaging tools to determine spatial patterns of performance.

Under the Hood: How It Works

The run_regional() function internally calls prep_regional() (from regional.R) to process the region mask, and then uses mvpa_iterate() to apply the MVPA model across each ROI. Functions such as combine_regional_results() and combine_prediction_tables() merge the individual regional outputs into a comprehensive result.

This modular design, laid out in mvpa_model.R and dataset.R, ensures that:

Data integrity is maintained.
Cross-validation is properly applied.
Results are aggregated for clear interpretation at the regional level.

Summary

This vignette showed you how to generate synthetic neuroimaging data, define ROIs with region masks, set up MVPA models with cross-validation, run analyses across ROIs, and interpret the performance metrics. You now have the tools to conduct regional MVPA analyses on your own neuroimaging data.

For further details, please refer to the source files:

regional.R for regional analysis methods.
mvpa_model.R for MVPA model creation and result formatting.
dataset.R for dataset generation routines.

Happy analyzing!

Bradley Buchsbaum

2025-03-20