Searchlight Analysis

## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display

## Warning: 'rgl.init' failed, will use the null device.
## See '?rgl.useNULL' for ways to avoid this warning.

Simulated data

We first create a small volumetric dataset for demonstration. The call below generates a 4D image with 6×6×6 spatial voxels and 80 observations along the time dimension. Observations are grouped into 4 blocks (20 trials each), and the response factor y has two levels. The return value contains both the data (mvpa_dataset) and the design (mvpa_design).

dataset <- gen_sample_dataset(D=c(6,6,6), nobs = 80, blocks=4, nlevels=2)
print(dataset)

## $dataset
## 
##  MVPA Dataset 
## 
## - Training Data 
##   - Dimensions:  6 x 6 x 6 x 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

Cross‑validation

Because trials are organized into runs (blocks), we use a blocked cross‑validation scheme: train on k−1 blocks and test on the held‑out block. This leave‑one‑group‑out strategy respects temporal correlations within runs and yields a more realistic estimate of generalization.

block <- dataset$design$block_var
crossval <- blocked_cross_validation(block)
crossval

## 
##  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

Model specification

We use the shrinkage discriminant analysis model (sda_notune), which estimates its shrinkage parameter from the training folds. See the sda package for details.

sda_model <- load_model("sda_notune") 
model <- mvpa_model(model=sda_model, dataset=dataset$dataset, design=dataset$design, crossval=crossval)
model

## mvpa_model object. 
## model:  sda_notune 
## model type:  classification 
## 
##  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 x 6 x 6 x 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

Standard searchlight

run_searchlight() returns image volumes with performance metrics at each sphere center. For two‑class problems we report cross‑validated accuracy and AUC (centered at 0 by subtracting 0.5). The radius is in millimeters; method = "standard" evaluates one model per sphere centered on each voxel.

result <- run_searchlight(model, radius=4, method="standard")
result

Randomized searchlight

The randomized variant samples non‑overlapping spheres and propagates each sphere’s score to all voxels it covers. Repeating this for niter exhaustive passes yields, at each voxel, the average performance across the spheres that included it. This can improve spatial localization and often reduces computation, since the number of models scales with nvoxels / radius × niter rather than the total number of voxels.

result <- run_searchlight(model, radius=4, method="randomized", niter=8)
result

Using different classifiers

Any classifier in the model registry can be used in a searchlight. You can also register your own with register_mvpa_model(). For example, to run a linear SVM:

svm_model <- load_model("svmLinear") 
model <- mvpa_model(model=svm_model, dataset=dataset$dataset, design=dataset$design, crossval=crossval)
result_svm <- run_searchlight(model, radius=4, method="randomized", niter=2)

For random forests, you can either fix mtry or tune it over a small grid. Tuning quickly becomes expensive in whole‑brain analyses, so sda_notune is often a good default when compute is limited.

if (requireNamespace("randomForest", quietly = TRUE)) {
  rf_model <- load_model("rf")
  model <- mvpa_model(model = rf_model, dataset = dataset$dataset, design = dataset$design,
                      crossval = crossval, tune_grid = data.frame(mtry = 2))
  result_rf <- run_searchlight(model, radius = 4, method = "randomized", niter = 2)
  result_rf
} else {
  message("Package 'randomForest' not installed; skipping RF example.")
}

if (requireNamespace("randomForest", quietly = TRUE)) {
  grid <- data.frame(mtry = c(2, 4, 6, 8))
  model2 <- mvpa_model(model = rf_model, dataset = dataset$dataset, design = dataset$design,
                       crossval = crossval, tune_grid = grid, tune_reps = 2)
  result_rf_tuned <- run_searchlight(model2, radius = 6, method = "randomized", niter = 1)
  result_rf_tuned
} else {
  message("Package 'randomForest' not installed; skipping tuned RF example.")
}

Searchlight Analysis

Bradley Buchsbaum

2025-09-28