Voxel or Vertex-Based GAMLSS

VBGAMLSS fits Generalized Additive Models for Location, Scale and Shape voxel-wise or vertex-wise, designed for normative modelling in neuroimaging. VBGAMLSS is based on gamlss2, so it uses the same formula syntax.

📦 Requirements

ANTsR
doFuture
gamlss
gamlss2
itertools
pbmcapply
progressr
tibble
qs2
future.mirai

---

🔧 Installation

Install gamlss2 and ANTsR

# devtools is required to install from GitHub
install.packages("devtools")

# gamlss2 from R-universe and ANTsR from GitHub
install.packages(
  "gamlss2", 
  repos = c("https://gamlss-dev.R-universe.dev", 
            "https://cloud.R-project.org")
)

# qs2: quick serialization 2
remotes::install_cran("qs2", type = "source", configure.args = "--with-simd=AVX2")

# Note for vertex based data ANTsR is not mandatory
devtools::install_github("ANTsX/ANTsR")

Install VBGAMLSS

devtools::install_github("tmspvn/VBGAMLSS", dependencies = TRUE)

Note: VBGAMLSS is heavily based on gamlss2 which is under development so changes are common. If you encounter issues, please let me know.

---

🎭 Quick Start

library(VBGAMLSS)

# Paths to example data
nsubj <- 258
img_controls <- "~/controls.nii.gz"  # 90x90x90x258
img_patients <- "~/patients.nii.gz"  # 90x90x90x258
mask <- "~/mask.nii.gz"      # 90x90x90


# Covariates
covs_controls <- data.frame(
  x0  = 1:nsubj,
  x1 = as.factor(rbinom(nsubj, 1, 0.6)),
)

covs_patients <- data.frame(
  x0  = rnorm(nsubj),
  x1 = rnorm(nsubj) * rnorm(nsubj),
)


# Convert to 2D subject × voxel/vertex
imageframe_controls <- images2matrix(img_controls, mask)
imageframe_patients <- images2matrix(img_patients, mask)


# Fit the voxel-wise GAMLSS models. Note: formulas follow gamlss2 syntax, 
#  but the Y is the response variable placeholder is mandatory. 
vbgamlss_model <- vbgamlss(
  imageframe_controls,
  g.formula = Y ~ pb(x0) + x1 | x1,
  g.family  = NO,
  num_cores = 20,
  train.data = covs_controls,
  debug = TRUE
)

After few hours (less than 8h for 150k voxels on 40 cores)

# Save / load
save_model(vbgamlss_model, "~/vbgamlss.model/fitted_model")
models_loaded <- load_model("~/vbgamlss.model/fitted_model.vbgamlss")


# Predict & compute Z-scores
predictions <- predict.vbgamlss(models_loaded, newdata = covs_patients)
zscores <- zscore.vbgamlss(predictions, imageframe_patients)


# Map the z-score to nifti. `_subj-<ID>.zscore.nii.gz` is then appended. 
#  ID id the patient index of the imageframe. 
#  Note, this function computes z-scores from a `vbgamlss.predictions" object. 
#  It differes from the `map_zscores_from_map` function which computes the maps 
#  z-scores directly from coefficient maps.
map_zscores(
  zscores,
  mask = mask,
  filename = "~/patients_deviation_maps/deviation_map",
)


# Map the predicted μ,σ,ν,τ response distribution parameters to nifti. `_subj-<ID>_fam-<FAMILY>_par-<PARAM>.nii.gz` is then appended. FAMILY is the distribution family (e.g. SHASH), PARAM is the parameter (e.g. mu, sigma).
map_model_predictions(
  predictions,
  mask = mask,
  filename = "~/patients_parameter_maps/prediction", 
)

---

📚 Functions Overview

Core

Function	Description
vbgamlss	Fit GAMLSS voxel/vertex-wise with optional segmentation and parallel processing. Process the data in chunks.

Support & Utility

Function	Description
images2matrix	Convert 4D NIfTI/list of 3D images to subject × voxel matrix.
load_input_image	Load image/matrix from NIfTI, RDS, or data frame.
save_model	Save fitted vbgamlss models to file.
load_model	Load saved vbgamlss models from file.
predict.vbgamlss	Predict parameters or responses from fitted models.
zscore.vbgamlss	Compute per-voxel z-scores.
restore_family	For space efficiency reason, vbgamlss remove the distribution class from each submodel after fitting. use this function to restore it if needed (e.g. if you want to do custom operation on vbgamlss submodels).
vbapply	Apply a function to each voxel submodel in a vbgamlss object and return a vector. restore_family() is already included in each iteration before calling the custom function.

Mapping

Function	Description
map_model_coefficients	Save per-voxel model coefficients to NIfTI images (e.g. β coeff. of a regression model as maps).
map_model_predictions	Save predicted parameters to NIfTI images (i.e. save μ,σ,ν,τ distribution coeff. as a maps).
map_zscores	Save z-score maps to NIfTI images.
map_zscores_from_map	Compute and save z-score directly from coefficients maps (i.e. after map_model_predictions).

Cross-validation

Function	Description
vbgamlss.cv	Perform stratified k-fold cross-validation for voxel-wise models and summaries results.

---

📝 General notes

• Using an HPC is strongly encouraged. Fitting models voxel-wise is computationally and memory intensive and may take a long time.
• Please consider tweaking chunk_max_mb parameter of vbgamlss to fit your system's memory constraints (256Mb is a good starting point). The function will proceed chunking the imageframe and parallize a chunk at the time to avoid loading everyhting in memory at once.
• When performing model selection is advise to do it on a subset of voxels (regularly or randomly sampled) with a penalized likelihood (e.g. GAIC with k=2 or k=log(n)) method instead of CV for time reasons.
• Providing an already subsampled mask is simple and effective. vbgamlss's subsample parameter is DEPRECATED.
• vbgamlss can be forced to save each chunk and resume the fit from the last chunk in case of interruptions (e.g. time limits on HPC) by providing cache=T and cachedir. See ?vbgamlss for details.

---

🏗️ Work in progress

• vbgamlss.cv: ⚠️ implemented, but not fully tested. Undocumented.
• vbgamlss segmentation handling: ⚠️ implemented, but not fully tested.
• vbgamlss.model_selection: ❌ "It will not work on other HPC systems. Undocumented.

---

⚠ Known Issues

• It takes a long time to fit models voxel-wise.
• pbmcappyly runs sequentially on some systems (i.e. Windows).
• vbgamlss.model_selection currently does not generalize across all HPC systems

---

All Functions (click to expand)

Core

Function	Description
vbgamlss	Fit GAMLSS voxel/vertex-wise with optional segmentation and parallel processing.

Cross-validation

Function	Description
vbgamlss.cv	Perform stratified k-fold cross-validation for voxel-wise models and summarise deviance.
predictGD	Predict parameters/responses and compute global deviance per voxel.
testGD	Compute deviance, prediction error, and residuals for a voxel.
statGD	Summarise deviance and prediction error across voxels.
describe_stats	Return mean, SD, quantiles, min, max for a vector.
stratCVfolds	Generate stratified fold assignments from a factor.
getCVGD	Aggregate global deviance across CV folds.
getCVGD.pen	Aggregate penalised global deviance across CV folds.
getCVGD.all	Aggregate both penalised and unpenalised deviance.
akaike_weights	Compute AIC/Akaike model weights.

Mapping

Function	Description
map_model_coefficients	Save per-voxel model coefficients to NIfTI images (e.g. β coeff. of a regression model as map).
map_model_predictions	Save predicted parameters to NIfTI images (i.e. save μ,σ,ν,τ distribution coeff. as a map).
map_zscores	Save z-score maps to NIfTI images.

Model-selection

Function	Description
vbgamlss.model_selection	Run multi-model CV jobs on HPC (Slurm).
slurm_template	Return a Slurm job script template.
slurm_resources	Build Slurm resource parameter list.
slurm_registry	Create registry for job tracking.
sanity_check	Verify required files exist before running.
sbatch_jobs	Submit jobs to Slurm.
jobs_status	Query Slurm job statuses.
monitor_jobs	Monitor jobs until completion.
gather_jobs_outputs	Collect outputs from all jobs.

Support

Function	Description
save_model	Save fitted vbgamlss models to file.
load_model	Load saved vbgamlss models from file.
predict.vbgamlss	Predict parameters or responses from fitted models.
zscore.vbgamlss	Compute per-voxel z-scores.

Utilities

Function	Description
images2matrix	Convert 4D NIfTI/list of 3D images to subject × voxel matrix.
load_input_image	Load image/matrix from NIfTI, RDS, or data frame.
restore_family	For space efficiency reason, vbgamlss remove the distribution class from each submodel after fitting. use this function to restore it if needed (e.g. if you want to do custom operation on vbgamlss submodels).
vbapply	Apply a function to each voxel submodel in a vbgamlss object and return a vector. restore_family() is already included in each iteration before calling the custom function.
estimate_nchunks	Calculate chunks to fit memory constraints.
get_subsample_indices	Generate indices for subsampling.
combine_formulas_gamlss2	Merge formulas for mu, sigma, nu, tau models.
quite	Suppress console output of an expression.
rand_names	Generate random string IDs.
check_formula_LHS	Ensure formula LHS is Y.
TRY	Try-catch with error/warning logging.