BLMM: Parallelised computing for big linear mixed models

Abstract

Within neuroimaging large-scale, shared datasets are becoming increasingly commonplace, challenging existing tools both in terms of overall scale and complexity of the study designs. As sample sizes grow, researchers are presented with new opportunities to detect and account for grouping factors and covariance structure present in large experimental designs. In particular, standard linear model methods cannot account for the covariance and grouping structures present in large datasets, and the existing linear mixed models (LMM) tools are neither scalable nor exploit the computational speed-ups afforded by vectorisation of computations over voxels. Further, nearly all existing tools for imaging (fixed or mixed effect) do not account for variability in the patterns of missing data near cortical boundaries and the edge of the brain, and instead omit any voxels with any missing data. Yet in the large-n setting, such a voxel-wise deletion missing data strategy leads to severe shrinkage of the final analysis mask. To counter these issues, we describe the "Big" Linear Mixed Models (BLMM) toolbox, an efficient Python package for large-scale fMRI LMM analyses. BLMM is designed for use on high performance computing clusters and utilizes a Fisher Scoring procedure made possible by derivations for the LMM Fisher information matrix and score vectors derived in our previous work, Maullin-Sapey and Nichols (2021).

Fig. 1

Activity diagram detailing the BLMM pipeline. The boundary of the BLMM code is indicated by the gray outline. The start and end nodes of the pipeline are represented by the black circle and nested black and white circles, respectively. Decision nodes are represented by diamonds and parallel stages of computation are represented with vertical bars. Also included are dotted lines indicating the distinct “stages” of the BLMM pipeline. Of particular note are the image-wise and voxel-wise batching stages of the pipeline, in which computation is parallelised across $B_I$ groups of images and $B_v$ groups of voxels, respectively.

Algorithm 1

Product Form Computation Pseudocode.

Algorithm 2

Full Simplified Fisher Scoring Pseudocode.

Fig. 2

A visual representation of the pipeline employed to generate the simulated data of Section 3.1. The first box depicts the model which was used for data generation, notably highlighting that $ϵ$ and $b$ varied across space. The second box details the smoothing process, with the $\otimes$ symbol representing convolution in this instance. The third box details the masking stage, with $\odot$ representing the Hadamard (element-wise) product.

Fig. 3

Observed serial computation times for each experimental design, displayed as a function of the number of observations, $n$. Displayed are the SCT in kiloseconds for BLMM (dashed) and lmer (dotted).

Fig. 4

First row: The MNI152 2mm anatomical template, for reference. Second row: ${\chi }^2$ statistics for comparison of model 1 and model 2, displayed on the square root scale; outlined in black are regions where evidence was found that inclusion of a random subject intercept significantly affected (at the $5\%$ Bonferroni-significance level) the results of the analysis. Third row: ${\chi }^2$ statistics for comparison of model 2 and model 3, displayed on the square root scale; outlined regions indicate where the inclusion of a random subject slope significantly affected the results of the analysis. Fourth row: Effect estimates for the “Faces $>$ Shapes” contrast, derived from model 2; for this row, voxels demarcated are those Bonferroni-significant at the $5\%$ level. Fifth row: The fixed effects standard deviation ($\sigma$), derived from model 2. Sixth row: The standard deviation of the subject-level random intercept ($\sigma \sqrt{d}$ where $d$ is the only non-zero element of $D$ in model 2).