Advancing statistical approaches, from data acquisition to inference testing

Spinal Muscular Atrophy (SMA) is a genetic disease that causes progressive dysfunction and death of spinal motor neurons, leading to motor deficits ranging from lower limb weakness (type 4) to severe muscle weakness with respiratory failure (type 1). Recent experiments in mice indicate that SMA motor deficits are due to motor neuron death and decreased firing rates in surviving motor neurons due to a maladaptive response to a loss in the excitatory Ia sensory synapses [1]. Epidural spinal cord stimulation (SCS) can selectively activate Ia sensory fibers; thus we hypothesize that targeted stimulation of Ia afferents via epidural SCS would increase inputs to the motor neurons, resulting in increased firing ability and improved leg functions through long-term stimulation effects (Figure 1A) [2-4]. To test the efficacy of our SCS therapy we quantified long term changes in motor neuron functions by performing functional magnetic resonance imaging (fMRI) of the lumbar spinal cord during active and passive mobilization of the knee joint pre- and post- SCS therapy (Figure 1B). Spinal cord fMRI is a rapidly growing field, but the lumbar spine has largely been ignored. Therefore, we leverage recently developed cervical spinal cord fMRI techniques to create a robust lumbar spine acquisition and processing paradigm, which can be applied to any clinical population [5-7]. 

Recent work has leveraged massive datasets and advanced image harmonization algorithms to construct normative models of imaging-derived phenotypes (IDPs)[1,2]. These brain chart models, which can produce centile or z-scores to benchmark individuals' morphology within a population, are often fit on magnetic resonance imaging data collected across hundreds of scanners. One popular method for harmonizing these data is ComBat, which preserves the effects of specified covariates on the IDPs' means. However, evidence suggests that biological factors, such as sex, also impact an IDP's variance across a population[3]. These scale effects, which directly impact centile and z-score distributions, are not preserved by current harmonization methods. Thus, harmonization may induce error in centile and z-scores, particularly when factors that impact scale are distributed unequally across sites.

Here, we propose a new method in the ComBat family of harmonization tools, ComBatLS, that preserves biological variance in IDPs' location and scale. We tested ComBatLS's ability to preserve variation in scale and its impacts on centile and z-scores by harmonizing across sex-imbalanced artificial "sites" in data from the UK Biobank. 

Spatial Bayesian models are a powerful way to account for spatial dependencies in fMRI analysis [3]. While massive univariate modeling treats each voxel or vertex as a separate entity, spatial Bayesian models place a multivariate prior distribution on the underlying maps of activation, which encodes the spatial dependence and implicitly smoothes the activation estimates. This avoids ad-hoc data smoothing, which induces spatially dependent noise and can lead to false positive clusters [2]. In addition, spatial Bayesian models can use the joint posterior distribution across brain locations to identify areas of activation. This dramatically increases power to detect activations and facilitates the use of meaningful minimum effect sizes, even in individual-level analysis [5]. The open-source BayesfMRI R package provides a user-friendly interface for spatial Bayesian models for task fMRI analysis.

Importantly, the spatial Bayesian models implemented in BayesfMRI are surface-based and subcortical parcel-constrained. These grayordinates-based models leverage spatial dependencies in a neurobiologically appropriate way and avoid the mixing of signals known to occur with whole-brain volumetric smoothing or spatial modeling [1]. 

Researchers using task-fMRI data have access to a wide range of analysis tools to model brain activity. This diversity of analytical approaches has been shown to have substantial effects on neuroimaging results (Botvinik-Nezer et al., 2020; Bowring et al., 2018; Carp, 2012; Glatard et al., 2015). Combined with selective reporting, this analytical flexibility can lead to an inflated rate of false positives and contributes to the irreproducibility of neuroimaging findings (Poldrack et al., 2017). Multiverse analyses are a way to systematically explore and integrate pipeline variation on a given dataset. We focus on the setting where multiple statistic maps are produced as an output of a set of analyses. Meta-analysis is a natural approach to extract consensus inferences from these maps, yet the traditional assumption of independence amongst input datasets does not hold. In this work we consider a suite of methods to conduct meta-analysis in the multiverse setting, accounting for inter-pipeline dependence among the results. 

Cluster-size inference is one of the most popular approaches in neuroimaging analysis but has been criticized for its dependence on an arbitrary (but important) cluster-forming threshold. Threshold-Free Cluster Enhancement (TFCE) was proposed in Smith (2009) with the aim of 1) decreasing the arbitrariness of the cluster-forming threshold and 2) being more sensitive to brain activation. Here we will show that TFCE does not achieve these aims. First, TFCE introduces more free parameters into the model than cluster-size inference, making it more dependent on researchers' choices. Second, the large clusters of TFCE, often mistakenly interpreted voxel-wise, actually hamper interpretation as the spatial specificity paradox increases. We will show, using recent advances in multiple testing, that cluster-size inference is a more robust way of doing neuroimaging analyses. 

The brain expresses activity in complex spatiotemporal patterns, reflecting cytoarchitectural and genetic influences that possess specific spatial properties. These brain patterns, also known as brain maps, frequently have high smoothness and spatial organization, i.e., spatial autocorrelation (SA), reflecting its central position in modern neuroimaging analyses [1–6]. In regimes of high SA, correlation between two brain maps can be spuriously elevated leading to false positive associations. An appropriate null hypothesis test to exclude false positives requires surrogate brain maps that preserve SA. Here we introduce "eigenstrapping", a technique for generating null hypotheses for maps possessing SA. This method uses geometric eigenmodes derived from various surfaces to produce surrogate brain maps that preserve SA. We show that these surrogate maps appropriately represent the null distribution and control false positives for cortical maps with SA, providing a versatile approach for investigations of cortical and subcortical topography.