ModelArray: An R package for statistical analysis of fixel-wise data

Abstract

Diffusion MRI is the dominant non-invasive imaging method used to characterize white matter organization in health and disease. Increasingly, fiber-specific properties within a voxel are analyzed using fixels. While tools for conducting statistical analyses of fixel-wise data exist, currently available tools support only a limited number of statistical models. Here we introduce ModelArray, an R package for mass-univariate statistical analysis of fixel-wise data. At present, ModelArray supports linear models as well as generalized additive models (GAMs), which are particularly useful for studying nonlinear effects in lifespan data. In addition, ModelArray also aims for scalable analysis. With only several lines of code, even large fixel-wise datasets can be analyzed using a standard personal computer. Detailed memory profiling revealed that ModelArray required only limited memory even for large datasets. As an example, we applied ModelArray to fixel-wise data derived from diffusion images acquired as part of the Philadelphia Neurodevelopmental Cohort (n = 938). ModelArray revealed anticipated nonlinear developmental effects in white matter. Moving forward, ModelArray is supported by an open-source software development model that can incorporate additional statistical models and other imaging data types. Taken together, ModelArray provides a flexible and efficient platform for statistical analysis of fixel-wise data.

Keywords: Big data; Development; Fixel-based analysis; MRI; Software; Statistical analysis.

Fig. 1.

Schematic of ModelArray and its companion converter ConFixel. The original fixel-wise data (.mif files; see the example fixel-wise image) are first converted into an HDF5 file (.h5) using ConFixel (top of the left box). ModelArray provides easy access to fixel-wise data in the HDF5 file (“accessor”). When performing statistical analysis of each fixel (top of the right box), to reduce memory usage, only a limited block of fixel-wise data is read into the memory. Using the phenotypes of interest (e.g.,: age, sex; provided by a CSV file), ModelArray fits a statistical model and calculates statistical output for each fixel. After iterating across fixels, the result matrix is generated (bottom of the right box) and saved to the original HDF5 file on disk by ModelArray (“write”). Finally, ConFixel converts the result matrix in this HDF5 file into a list of .mif files ready to be viewed (bottom of the left box). The fixels in the fixel-wise image are colored by fixel orientation (red: left–right, green: anterior–posterior, blue: inferior–superior), and background image is the fiber orientation distribution (FOD) template.

Fig. 2.

Example R code for executing analysis using ModelArray. ModelArray functions are highlighted in green.

Fig. 3.

Memory required by ModelArray does not vary by sample size. The maximal memory required by a linear model executed using ModelArray.lm() was evaluated when analyzing a range of sample sizes. All models were performed with a parallelization factor of 4.

Fig. 4.

ModelArray is memory-efficient even under different parallelization configurations. Maximal memory usage for a linear model run using ModelArray.lm() was evaluated across a sample of n = 30 (A) and n = 938 (B) with varying numbers of CPU cores requested.

Fig. 5.

ModelArray allows the estimation of nonlinear effects. Fixel-wise GAM fitted with ModelArray.gam() revealed nonlinear FDC changes with age in childhood and adolescence (n = 938). The GAM also included sex and motion quantification as covariates. (A) Fixels whose FDC was significantly associated with age (p-value of s(age) < 1 × 10⁻¹⁵); fixels are colored by effect size of s(age). Background image is the FOD template. (B) GAM fit for FDC averaged in the 2D slice of the cluster in CC highlighted in panel A by a white arrow.