Regional Analysis
Bradley Buchsbaum
2024-04-23
Regional_Analysis.RmdMVPA Analysis on Regions of Interest
If we have one or more pre-defined regions of interests (ROIs) we can
use rMVPA to run a clasisifation/regression model in each
region, yielding peformance measures for each region.
Generate a volumetric dataset with 100 observations and two classes
To generate a dataset we use the gen_sample_dataset
function. We are creating a 4-dimensional neuroimaging dataset, with
6-by-6-by-6 spatial dimensions and 80 observations in the 4th dimension.
These 80 observations are divided into 4 blocks, each consisting of 20
trials. The generated y variable is a factor
with 2 levels (‘a’ and ‘b’). the gen_sample_dataset
function creates a list with two elements: an mvpa_dataset
object (dataset) and an mvpa_design object
(design). The first contains information about the data
itself and the second contains information about the experimental
design.
dset <- gen_sample_dataset(D=c(6,6,6), nobs = 80, blocks=4, nlevels=2)
print(dset)## $dataset
## $train_data
## DenseNeuroVec
## Type : DenseNeuroVec
## Dimension : 6 6 6 80
## Spacing : 1 X 1 X 1
## Origin : 0 X 0 X 0
## Axes : Left-to-Right Posterior-to-Anterior Inferior-to-Superior
## Coordinate Transform : 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
##
## $test_data
## NULL
##
## $mask
## NeuroVol
## Type : LogicalNeuroVol
## Dimension : 6 6 6
## Spacing : 1 X 1 X 1
## Origin : 0 X 0 X 0
## Axes : Left-to-Right Posterior-to-Anterior Inferior-to-Superior
##
## attr(,"class")
## [1] "mvpa_image_dataset" "mvpa_dataset" "list"
##
## $design
## mvpa_design:
## training observations: 80
## training levels: a b
## no test observations.
## training response: Factor w/ 2 levels "a","b": 2 1 1 2 2 1 2 1 1 1 ...
## block var: int [1:80] 1 1 1 1 1 1 1 1 1 1 ...Now we generate an arbitrary set of ROIs by assigning each voxel in the input mask a value ranging from 1 to 3.
mask <- dset$dataset$mask
nvox <- sum(mask)
region_mask <- neuroim2::NeuroVol(sample(1:3,replace=TRUE,size=nvox), space(mask), indices=which(mask>0))
table(region_mask)## region_mask
## 0 1 2 3
## 123 30 33 30The region_mask volume now defines a set of three ROIs
where each ROI is defined by an integer value. Next, we can run a
classification/regression model within each ROI using the
run_regional function:
mod <- load_model("sda")
tune_grid <- data.frame(lambda=.01, diagonal=FALSE)
cval <- blocked_cross_validation(dset$design$block_var)
mspec <- mvpa_model(mod, dataset=dset$dataset, design=dset$design, tune_grid=tune_grid,crossval=cval)
res <- run_regional(mspec,region_mask)run_regional produces several output values, including
performance_table, which includes cross-validated
performance metrics for each ROI:
print(res$performance_table)## # A tibble: 3 × 3
## roinum Accuracy AUC
## <int> <dbl> <dbl>
## 1 1 0.488 -0.0300
## 2 2 0.525 0.0437
## 3 3 0.588 0.149To get fine-grained information, the individual predictions for each
ROI and trial are included in prediction_table:
print(res$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.000426 a FALSE 1.00 0.000426
## 2 2 1 a 0.00602 b FALSE 0.00602 0.994
## 3 3 1 a 0.988 a TRUE 0.988 0.0119
## 4 4 1 b 0.0146 a FALSE 0.985 0.0146
## 5 5 1 b 0.0231 a FALSE 0.977 0.0231
## 6 6 1 a 0.990 a TRUE 0.990 0.0103
## 7 7 1 b 0.997 b TRUE 0.00264 0.997
## 8 8 1 a 0.0454 b FALSE 0.0454 0.955
## 9 9 1 a 0.0113 b FALSE 0.0113 0.989
## 10 10 1 a 0.881 a TRUE 0.881 0.119
## # ℹ 230 more rows