A sensitivity analysis of preprocessing pipelines: Toward a solution for multiverse analyses Open Access

Modern neuroimaging techniques have provided unique opportunities to measure complex signaling pathways in the living human brain with the goal of identifying reliable biomarkers of disease states and treatment outcomes (Poldrack et al., 2020). Data arising from state-of-the-art neuroimaging techniques are, however, often contaminated with noise confounds such as motion-related artefacts, affecting both the spatial and temporal correlation structure of the data (Nørgaard, Ozenne, et al., 2019). Carefully designed preprocessing steps have been developed to remove unwanted noise sources, but in the absence of a “ground truth” it remains a major challenge to evaluate the impact of preprocessing choices on subsequent statistical analyses and results. Over time, preprocessing pipelines (i.e., a set of preprocessing steps) have become more complex and flexible, and this increase in researcher degrees of freedom (termed multiverse analyses) has consistently been shown to affect the outcomes of neuroimaging studies (Botvinik-Nezer et al., 2020; Carp, 2012; Nørgaard, Ganz, et al., 2019). The most common approach in the neuroimaging field is, to date, to use a single pipeline and not consider the heterogeneity of preprocessing choices. This approach not only makes abstraction of the multitude of possible results more difficult, but is likely also sub-optimal because the best pipeline is more often than not, study, population, or even subject dependent (Churchill et al., 2015; Nørgaard, Ganz, et al., 2019). More concerning, neuroscientists might be tempted to “tune” the pipeline in order to obtain the most satisfying results. This will generally lead to spurious and non-reproducible results since the variability induced by the choice of pipeline is now conditioned on the desired outcome. However, since it is neither realistic nor optimal to move toward a single unified preprocessing pipeline, there is an urgent need for a statistical framework allowing to explore results among many preprocessing pipelines in a principled way.

The aim of this work is, thus, to provide a statistical framework that can aggregate the evidence from multiverse analyses to produce conclusions robust to the choice of the pipeline. More specifically, the present paper proposes a statistical sensitivity analysis providing:

• (i)

  visualizations of the heterogeneity of several preprocessing pipelines

• (ii)

  estimation of a global effect across all preprocessing pipelines

• (iii)

  quantification of the proportion of pipelines with evidence for an effect

• (iv)

  a statistical framework for testing hypotheses across pipelines such as “no effect across all pipelines”, “at least one pipeline with no effect”

In this work, as an exemplary case, we use Positron Emission Tomography (PET) neuroimaging data to carry out a real-world multiverse analysis. However, our framework is generic and may be used in any multiverse scenario. The corresponding software is available as an R package called LMMstar freely accessible on CRAN, and our Github repository contains the code and data to reproduce all simulations and figures. Finally, we also provide a step-by-step tutorial on GitHub so other researchers can carry out these analyses using their own data.

We used two different data sources in our analyses: in sillico and real data. For the in sillico data, different noise structures were chosen to reflect how different configurations of pipelines may interact with the distribution of the data, and the sample size was varied to encompass small to larger scale clinical studies. In the real data analysis, we mimicked how real neuroimaging studies compare an intervention to a reference measurement. Pipelines were selected independently of the intervention data using the healthy/placebo arm of Frokjaer et al. (2015) using results from Nørgaard, Ganz, et al. (2019), investigating the impact of preprocessing pipelines. The proposed sensitivity analysis was illustrated on the intervention arm of the study where the follow-up value was compared to a reference value taken from a normative serotonergic atlas (Beliveau et al., 2017).

To simulate in sillico data, we considered the simple case of a single brain measurement (⁠$R=1$⁠), with a single binary exposure (⁠$P=1$⁠) following a Bernoulli distribution with parameter $π=0.5$ (i.e., two balanced groups) and no covariates (⁠$C=0$⁠). Latent $Y$ values for the brain measurement were simulated using a normal distribution with variance 1 and mean $β$ times the exposure $X$ values, where $β=0$ (null hypothesis) or $β=0.5$ (alternative hypothesis). The observed $Y$ was simulated for each pipeline, adding pipeline specific noise to the latent $Y$⁠. This noise was simulated using a multivariate normal distribution with mean 0 and variance $∑scenario 1$⁠, $∑scenario 2$⁠, and $∑scenario 3$ depending on the scenario (see Supplementary Material A for details). In scenario 1, we simulated many pipelines (⁠$J=20$⁠) with correlated homoscedastic noise; in scenario 2, a few pipelines (⁠$J=6$⁠) with uncorrelated heteroscedastic noise; and in scenario 3, many pipelines (⁠$J=20$⁠) with correlated heteroscedastic noise. Scenario 2 and 3 included one pipeline with high signal to noise ratio (SNR), that is, low variance, and many pipelines with low SNR, that is, high variance. Additional simulations with non-normally distributed noise (⁠$∑scenario 4$ and $∑scenario 5$⁠) can be found in Supplementary Material E, investigating the impact of outliers (⁠$∑scenario 4$⁠) and skewed data (⁠$∑scenario 5$⁠). The sample size was varied from $n=10$ to $n=500$ in each group, such that the smallest sample size was well below the number of parameters (⁠$2×J$ mean parameters, $J2$ variance-covariance parameters) in scenario 1 and 3 and the asymptotic regime was reached for the largest sample size. We generated 10,000 datasets per scenario and sample size—this provides sufficient precision about the mean, standard deviation, and rejection rate to neglect the Monte Carlo uncertainty. These data will be used to assess the large sample size properties of the procedure (bias, variance, type 1 error control) in finite samples.

To illustrate the use of the proposed sensitivity analysis on real data, we utilize neuroimaging results from a placebo-controlled, double-blinded, clinical study (Frokjaer et al., 2015). The study was registered and approved by the ethics committee for the capital region of Copenhagen (protocol-ID: H-2-2010-108) and registered as a clinical trial: www.clinicaltrials.gov under the trial ID NCT02661789. All subjects provided written informed consent prior to participation, in accordance with The Declaration of Helsinki II. The aim of the study was to assess the association between the emergence of depressive symptoms and change in cerebral serotonin transporter (SERT) availability following a hormonal treatment (⁠$P=1$⁠). Raw data is available from the CIMBI database (https://nru.dk/index.php?view=category&id=128) upon request under GDPR law, and derived data is available on Github. It consists of structural Magnetic Resonance Imaging (MRI) and Positron Emission Tomography (PET) imaging data for $60$ healthy females who underwent a baseline scan, received either Placebo (⁠$n=30$⁠) or a GnRHa implant intervention (⁠$n=30$⁠), and participated in a follow-up scan. SERT availability estimates were modeled and extracted for each subject for 28 subcortical and cortical regions (see Nørgaard, Ganz, et al., 2019 for details), and averaged across hemispheres, producing a final sample of $R=14$ regions per subject and pipeline. These regions (amygdala, thalamus, putamen, caudate, anterior cingulate cortex, hippocampus, orbital frontal cortex, superior frontal cortex, occipital cortex, superior temporal gyrus, insula, inferior temporal gyrus, parietal cortex, and entorhinal cortex) were chosen because they cover the entire brain, and many are target regions in published serotonin transporter (SERT) PET studies. No covariates were considered (⁠$C=0$⁠).

Five preprocessing steps were used to analyze the data and estimate the SERT availability (outcome measure). These steps include motion correction (with/without), co-registration (4 options), delineation of volumes of interest (3 options), partial volume correction (4 options), and kinetic modeling for quantification of SERT availability (MRTM, SRTM, Non-invasive Logan and MRTM2). More information about the preprocessing choices and software can be found in (Nørgaard, Ganz, et al., 2019). The combination of individual preprocessing steps leads to a number of $J=2×3×43=384$ possible combinations.

We now introduce notations that will be used to describe the proposed sensitivity analysis (see Supplementary Material F for a summary). We are interested in relating $R$ brain measurements (⁠$Y=(Y1,…,YR)$⁠) to $P$ exposures or treatments (⁠$X=(X1,…,XP)$⁠) accounting for $C$ covariates (⁠$W=(W1,…,WC)$⁠). We consider a set of $J$ pipelines used to preprocess the neuroimaging data. For a given pipeline $j∈{1,…,J}$ we fit a statistical model with parameters $θj$ that relates $Y$ processed by pipeline $j$⁠, $X$⁠, and $W$⁠. We then obtain from this model an estimate $ψ^j$ of the effect of interest (denoted $ψ$⁠). In our real-life example, we use, for each pipeline, a paired t-test to compare the observed change in SERT availability to an atlas value so for $j∈{1,…,J}$⁠, $θ^j=(ψ^j,σ^j2)$ where $ψ^j$ is the empirical mean and $σj2$ the empirical variance of the change in SERT availability (processed with pipeline $j$⁠) between baseline and follow-up.

We make the following working assumptions: first, the observed data $(Oi)i∈{1,…,n}=(yi,xi,wi)i∈{1,…,n}$ correspond to independent and identically distributed replicates of $(Y,X,W)$⁠. Second, we have chosen a set of reasonable pipelines, meaning that the estimated effects $ψ^1,…,ψ^J$ found in the follow-up statistical analysis will converge to the right value $ψ$ as the sample size increases. This set can include pipelines distorting the signal $Y$ (e.g., adding a fixed value) if that has no consequence, asymptotically, on the estimated effect (the mean change is not biased as the added value cancels out when substracting the baseline value to the follow-up value). Finally, when considering asymptotic results we will consider a fixed number $J$ of pipelines and let the sample size $n$ increase to infinity.

To be able to draw conclusions across pipelines, we not only need the result of each pipeline but also some information about how they relate. If all pipelines were equally reliable and equally different, we would weight each pipeline equally. If there exists one independent pipeline and a block of correlated pipelines, all equally reliable, then we would weight the independent pipeline more compared with other pipelines. By treating pipelines as black boxes, we can investigate their relation in terms of each observation’s influence on the estimated effects across pipelines, $ψ^=(ψ^1,…,ψ^J)$⁠. This relation is fully characterized by the joint distribution of the effects. Once estimated, we can extract summaries of this distribution, for example, an average value, and carry out statistical tests, for example testing the compatibility between the observations and the joint distribution that would have been observed under a specific hypothesis.

The joint distribution could be obtained using a multivariate model, for example, modeling data from all pipelines at once using a mixed model. Because of the complexity of the dependency among pipelines, this is, however, rarely feasible with the available sample size. Instead, and this matches common practice, we perform the same analysis separately for each pipeline and obtain a vector of estimated associations $ψ^$ with their standard errors $σψ^=(σψ^1,…,σψ^J)$⁠. Using tools from the semi-parametric theory (see Kennedy, 2017 and Tsiatis, 2006 for more details), we can approximate the influence of each observation on the estimate by a random variable called the influence function, denoted $φψ^j$ for pipeline $j$⁠, and satisfying:

where $op(1)$ denotes a residual term that convergences toward zero in probability as the sample size tends to infinity. As shown in Supplementary Material B, $φψ^j$ has a simple expression when $ψ^j$ is the empirical mean or an element of a maximum likelihood (ML) estimator. Since this decomposition applies to all pipelines, we get from the multivariate central limit theorem that the joint distribution of the estimates is asymptotically multivariate normal. It has mean $ψ$ and its variance-covariance, denoted $∑ψ^$⁠, is the same as the one of $φ^ψ^=(φψ^1,…,φψ^J)$ divided by $n$⁠. Note that with limited number of observations, typically when $n<J$⁠, the estimated variance-covariance $Σ^Ψ^$ based on the estimated influence function is not guaranteed to be positive definite.

The global null hypothesis “no effect across all pipelines” can be tested using a max-test approach: more extreme realizations would correspond to larger values of the maximum statistic $t^max=max(|t^1|,..,|t^J|)$ where $|.|$ denotes the absolute value and $t^j=ψ^jσψ^j$⁠. A p-value may, therefore, be computed by integrating the joint density under the null hypothesis outside of the domain $D(t^max)=[−t^max, t^max]⊗J$ (see Supplementary Fig. B in Supplementary Material C). Here, we use the notation $[a,b]⊗J=∏j=1J[a,b]$ that represents the Cartesian product between $J$ intervals $[a,b]$⁠. The value $tc$⁠, such that the integral outside $D(tc)$ equals $α$⁠, provides a critical threshold for the estimated test statistics $(|t^1|,..,|t^J|)$⁠. This threshold can also be used to derive confidence intervals.

The null hypothesis “at least one pipeline with no effect” is an intersection union test. As such, it can be rejected if and only if all the un-adjusted p-values relative to each pipeline are below $α$ or equivalently if the largest un-adjusted p-value is below $α$⁠.

The proportion of pipelines where there is evidence for an effect $η$ can be estimated as $η^𝟙=1J∑j=1J𝟙|t^j|<tc$ where $%.$ denotes the indicator function. One drawback with this non-parametric estimator is that it is a non-smooth function of $t^j$⁠, making the associated uncertainty difficult to evaluate. Instead, one can use a parametric approach, assuming normally distributed test statistics:

which is a smooth (but complex) function of the model parameters. Here, $Φ$ refers to the cumulative distribution function of a standard normal distribution. The uncertainty about the estimator can, therefore, be derived using a non-parametric bootstrap or a delta method where $Var(η^Φ)=∂η^Φ∂Θ∑Θ^∂η^Φ⊺∂Θ$⁠, and $Θ=(θj)j∈{1,…,J}$ is the set of parameters of the statistical model across pipelines.

In section 3.1, it was stated that the estimated associations are, asymptotically, normally distributed. They can, therefore, be summarized by their expectation and variance-covariance matrix (i.e., standard errors and correlation matrix). In this work, we suggest two graphical displays to visualize the heterogeneity of the results across pipelines:

• A heatmap of the estimated correlation among estimates, obtained by converting $Σ^Ψ^$ into a correlation matrix.¹

• A forest plot displaying $(ψ^1,…,ψ^J)$ and $(σ^ψ^1,…,σ^ψ^J)$ through the confidence intervals, possibly using the previously established re-ordering.

Several methods can be used for estimating a global effect across pipelines. A naive method would be to compute the mean of the estimated associations:

This estimator will, however, not be efficient if some pipelines lead to more precise estimates, that is, the elements in $(σψ^1,…,σψ^J)$ are not equal. Intuitively, we would like to pool the estimates with weights inversely proportional to the standard errors such that we put more weight on precise estimates:

However, we also need to take into account the correlation between estimates. Indeed, perfectly correlated estimates should weight as if there was only one estimate. To do so, we can use the following GLS estimator of the global effect:

where $1$ is a column vector filled with ones and $∑j=1JwjGLS=1$⁠. Indeed, equation 1 can be shown to be equivalent to performing a spectral decomposition of $Σ^Ψ^$⁠, and using the eigenvectors to combine estimates into independent components that can be pooled according to weights proportional to the eigenvalues (Supplementary Material D1). This is used when $∑ψ^$ is singular to compute $Ψ^GLS$⁠, by restricting the spectral decomposition to the eigenvalues above a given threshold (⁠$∈=10−10$⁠). In the simple case where $R=1$⁠, $P=1$⁠, $C=0$⁠, brain measurements are jointly normally distributed, $X$ is binary, and there are no missing data, the GLS estimator can be shown to be asymptotically efficient (Supplementary Material D2).

Cox et al. (2006) studied a similar estimator when $J=2$ and found that, under unequal variance, the global estimate can be outside of the range of the (pipeline specific) estimates. As a remedy, we propose a constrained GLS estimator, denoted $Ψ^constrained GLS$⁠, which, in addition to have weights that sum up to 1, constrains the weight of each estimate to be at most 1 in absolute value:

where $κ$ is chosen to satisfy the constraint (Supplementary Material D3). These constrained weights can be seen as a regularization of the GLS weights (first term) toward the weights for the mean (second term) with regularization coefficient $1κ+maxj∈{1, …,J} |wjGLS|$$=J−1J maxj∈{1, …,J} |wjGLS|−1$⁠. This regularization coefficient is 1 when the largest weight is 1 (⁠$Ψ^constrained GLS=Ψ^GLS$⁠) and decreases toward 0 as the largest weight increases toward infinity (⁠$Ψ^constrained GLS→Ψ^average$⁠).

Note that if some of the pipelines may induce some bias, the previous approaches will propagate this bias and therefore be unsatisfactory. Systematic differences between pipelines can be investigated by comparing the estimates between pipelines, for example, $ψ^j−ψ^j′$ and using $Σ^Ψ^$ to obtain the corresponding uncertainty: $Var[ψ^j−ψ^j′]=Var[ψ^j]+Var[ψ^j′]−2ℂov(ψ^j,ψ^j′)$⁠.

We compare the performance of the proposed estimators when using unbiased pipelines on simulated data. For each dataset $∑scenario 1$⁠, $∑scenario 2$⁠, and $∑scenario 3$⁠, we computed the four previously described estimators of the global effect, $Ψ^average$⁠, $Ψ^pool-se$⁠, $Ψ^GLS$⁠, and $Ψ^constrained GLS$⁠. P-values for $Ψ^pool-se$⁠, $Ψ^GLS$⁠, and $Ψ^constrained GLS$ were computed, neglecting the uncertainty of the weights (⁠$wjse$⁠, $wjGLS$⁠, $wjconstrained GLS$⁠).

Based on $∑scenario 1$⁠, $∑scenario 2$⁠, and $∑scenario 3$⁠, we can compute how $Ψ^average$⁠, $Ψ^pool-se$⁠, and $Ψ^GLS$ (or $Ψ^constrained GLS$⁠) would weight the results from each pipeline if the variance-covariance matrix of the pipeline estimates was known (Fig. 1). In scenario 1, $Ψ^average$ and $Ψ^pool-se$ would provide equal weight to all pipelines with a weight of 5%, while $Ψ^GLS$ would weight by 1.88% the correlated pipelines and by 14.37% the uncorrelated pipelines (in this paragraph, all weights are rounded for readability). In scenario 2, $Ψ^average$ equally weights all pipelines by 16.67% while $Ψ^pool-se$ and $Ψ^GLS$ would use the following weights: 18.75%, 52.51%, 10.94%, 7.72%, 5.97%, and 4.10%, favoring the high SNR pipeline. In scenario 3, $Ψ^average$ would equally weight all pipelines by 5%, $Ψ^pool-se$ would weight each correlated pipeline by 5.17% and the remaining pipelines by 14.48% (high SNR pipeline), 3.02%, 2.13%, 1.65%, and 1.13%, while $Ψ^GLS$ would weight by 1.65% the correlated pipelines and by 48.61% (high SNR pipeline), 10.13%, 7.15%, 5.52%, and 3.80% the remaining pipelines.

The average estimate of the global effect estimators (⁠$Ψ^average$⁠, $Ψ^pool-se$⁠, $Ψ^GLS$⁠, $Ψ^constrained GLS$⁠) was close to the true value for all sample sizes and all scenarios, typically $±0.005$⁠, and the largest discrepancy was 0.015 (⁠$Ψ^GLS$⁠, n=10 per group, scenario 3, alternative hypothesis).

Under the null, the estimated proportion of pipelines with evidence for an effect ranged between 5.35% and 6.79% based on $η^Φ$ and between 0.50% and 1.71% based on $η^%$⁠. $η^Φ$ showed no obvious trend with sample size, whereas $η^%$ showed a small decrease with the sample size (upper panel of Supplementary Fig. C). Under the alternative, the proportion of pipelines dramatically increased with the sample size for both estimators (upper panel of Fig. 3): typically from 10% or below when $n=10$ per group to 99% in scenario 1, 65% in scenario 2, and around 80% in scenario 3 with $n=500$ per group. The “usual” estimator $η^%$ exhibited a smaller value at low sample size but a steeper increase compared to the proposed estimator $η^Φ$⁠.

Similar results were obtained for $β=0$ and $β=0.5$⁠, so we will only discuss the former case (upper panels of Fig. 2). In scenario 1, $Ψ^average$ and $Ψ^pool-se$ showed similar variability, whereas the variability of $Ψ^GLS$ was higher in very small samples (standard deviation: +156% for $n=10$⁠, +11% for $n=25$⁠) and lower for larger sample sizes (e.g., -11% for $n=500$⁠). In scenario 2, $Ψ^pool-se$ and $Ψ^GLS$ showed similar variability (slightly higher for $Ψ^GLS$ in low sample sizes and slightly lower afterward), both smaller than the variability of $Ψ^average$⁠, for example, -9.93% for $n=10$ and -24.86% for $n=500$ for $Ψ^GLS$⁠. Results in scenario 3 were similar to scenario 1 up to a larger decrease in variability in large samples for $Ψ^GLS$ (-27.75% for n=500). The variability of $Ψ^constrained GLS$ was similar to one of $Ψ^GLS$ except in small samples where it was smaller and closer to the one of $Ψ^average$⁠.

The proposed estimator $η^Φ$ had a lower variability compared to the usual estimator $η^1$ in all scenarios and across sample sizes under the alternative hypothesis (except in scenario 1, $n=500$⁠): the reduction in standard deviation ranged from 5.67% to 36.94% (lower panels of Fig. 3).

Across all scenarios, the type 1 errors for $Ψ^average$ and $Ψ^pool-se$ were well controlled except in very small samples where small deviations from the nominal level were observed (maximum of 6.31% for $Ψ^average$ and 7.36% for $Ψ^pool-se$⁠, lower panels of Fig. 2). When neglecting the uncertainty about the weights, the type 1 error control of $Ψ^GLS$ and $Ψ^constrained GLS$ was only controlled for large samples, that is, for $n=500$⁠; large type 1 error rates were found in very small samples (e.g., 22.02% for $Ψ^constrained GLS$ with $n=10$⁠).

While the previous section evaluated the performance of the sensitivity analysis assuming ground-truth data generated from a multivariate normal distribution, this assumption is often violated in real neuroimaging data (especially with patient data). To address this issue, we expanded the simulation scenarios to include non-normal data (scenario 4 and 5), allowing us to understand how well the sensitivity analysis works, how these data distributions interact with sample size, and what additional constraints are necessary to govern the usage of the sensitivity analyses. More specifically, we created two additional simulations largely corresponding to scenario 3 (many pipelines with correlated heteroscedastic noise), but with different noise distributions. Scenario 4 used a noise distribution following a Student’s t-distribution with 5 degrees of freedom, resulting in a distribution with heavy tails to reflect the inclusion of outliers. Scenario 5 used a noise distribution following a half-normal distribution, resulting in a heavily skewed non-negative distribution, with most of the density centered around 0. Figure 4 highlights the differences between scenario 3, 4, and 5, across different noise levels. Additional simulation details for scenario 4 and 5 can be found in the Supplementary Material E.

Notably, the results for these additional scenarios were similar to the multivariate normal cases in terms of bias, variability, and rejection rate (Supplementary Fig. E). However, while these additional scenarios may not fully represent real-world data, they show the robustness of the proposed sensitivity analysis to assumptions about the residual distribution. Intuitively, they benefit from the law of large numbers, that is, for large enough sample sizes they will behave as expected. Worse low sample size performances with non-Gaussian noise distributions is expected, which is in line with some of the simulation results (e.g., bias in scenario 5, Supplementary Fig. D). The difference between the Gaussian and non-Gaussian case is not striking though, and factors such as sample size seem to be driving the behavior of the estimators.

We illustrate the statistical sensitivity analysis with real data described in Section 2.1 in which the SERT availability was assessed after a drug intervention and compared to normative values. We start by studying the behavior of the four statistical estimators (pooled GLS, pooled constrained GLS, pooled average and pooled SE) to estimate a common effect across pipelines for the given null hypothesis for a single subcortical region - amygdala - first, reported as a forest plot in Figure 5 (left panel). The dashed vertical line in Figure 5 represents the normative value, and all horizontal error bars represent the estimated effect (mean and 95% CI) for a given pipeline and estimator. Across a reasonable set of preprocessing pipelines, three of the eight selected reject the null hypothesis (as indicated by the non-overlapping CI with the normative value) with estimated percent differences between groups ranging between -9% (pipeline 6) and 2.5% (pipeline 1). The pooled constrained GLS, pooled SE and pooled average all fail to reject a common effect across pipelines, whereas the GLS estimator rejects the null hypothesis. However, when inspecting the pipeline weights for the GLS estimator for these data, it assigned a very high weight to four pipelines (i.e., weight above 1 in absolute values), leading to an unreliable estimate as illustrated by the large standard deviation found in the simulation study for low sample sizes (Fig. 2). The constrained GLS estimator did not exhibit this problem and had weights between -0.79 and 1. We also visualize the effect of the four global effect estimators across the nine cortical regions. This is shown in Figure 6. Here, we can observe a similar pattern, where the pooled GLS estimator exhibited the largest deviations from the expected value based on the SERT normative values, whereas the pooled constrained GLS estimator is close to the pooled SE and pooled average.

Figure 5 (right panel) shows a heatmap for the estimated correlations across preprocessing pipelines, ranging from 0.68 to 1. Pipelines 1-4 show a very high correlation with each other (only varying the registration choices), whereas pipelines 5, 6, and 8 are equally correlated with each other (different kinetic models), and pipeline 7 is somewhat in between (no motion correction). The heatmap captures important differences in the correlation structure between pipelines, suggesting that not all pipelines perform similarly with moderate levels of unexplained variance. Pipelines 4 to 8 exhibit numerically smaller correlation compared to pipelines 1 to 4 (and similar variance), which explains why GLS estimators produce estimates with lower numerical values compared to pooled estimators ignoring the correlation, as GLS estimators assign more weight to the last four pipelines.

Notably, no correction for multiple comparisons was carried out across regions and pipelines at this point. The rationale for not including this (as should otherwise always be carried out) is that we wanted to make our analysis as comparable as possible to the Neuroimaging Analysis Replication and Prediction Study (Botvinik-Nezer et al., 2020), where each participating institution analyzed the data using their own established pipeline and tested only a single region in a hypothesis-driven fashion. Future work should address how adequate correction for multiple comparisons should be carried out across multiple estimators and pipelines.

The proposed statistical framework leverages information about the joint distribution of the test statistics across pipelines to provide insights about the impact of preprocessing strategies. In particular, we introduce a new estimator of the proportion of pipelines with evidence for an effect (⁠$η^Φ$⁠) and a new estimator of the effects across pipelines (⁠$Ψ^constrained GLS$⁠), able to pool correlated and heteroschedastics estimates while exhibiting reasonable small sample size properties. This information can complement a standard forest plot, as exemplified in Figure 5 where we added a heatmap of the correlation between pipelines.

The simulation results show that while $Ψ^GLS$ outperforms standard pooling techniques in large samples (⁠$Ψ^average$ or $Ψ^pool-se$⁠), it is not suited for the commonly observed cases in neuroimaging studies of smaller N (number of subjects) relative to the large P (number of pipelines). The proposed alternative estimator (⁠$Ψ^constrained GLS$⁠) can be seen as regularizing the GLS estimator toward the empirical average in small samples. Even though this considerably improves the small sample performance of the estimator, other regularization approaches may lead to further gain. An alternative approach could be to regularize the estimated variance-covariance matrix between pipeline-specific estimates, for example, using graphical lasso. However, this is challenging since the correlation structure among pipelines is typically complex and non-sparse.

The proposed estimator of the proportion of pipelines (⁠$η^Φ$⁠) showed lower variability in the simulation studies which should benefit statistical power. It is also intuitively appealing due to less sensitivity to the threshold chosen for statistical significance compared to averaging the number of pipelines for which the p-value was below, say, 0.05.

In the real-world application, we observed that all estimators could be readily applied, and three of them performed as expected based on the results from individual pipelines. Only the $Ψ^GLS$ estimator performed differently and was the only one rejecting the null hypothesis hinting at an effect across pipelines. This is, though, due to the estimator not being able to fit the weights properly in small sample sizes. Since neuroimaging studies in general and especially in PET are rarely beyond sample sizes of $n>50$ (Szucs & Ioannidis, 2020), other estimators should be used. We recommend using $Ψ^pool-se$ or $Ψ^constrained GLS$⁠: the latter when pipelines are not similarly related and the sample size is moderate to large, otherwise the former.

The framework that we are proposing is not without limitations. Fast and accurate estimation of statistical uncertainty is still being investigated. One attempt, neglecting the uncertainty of the weights, is numerically fast but unreliable in small samples or with a large number of pipelines for $Ψ^GLS$ and $Ψ^constrained GLS$⁠. For the moment, resampling methods such as non-parametric bootstrap or permutation test, while more computationally demanding, are recommended.

The underlying assumption of combining results across pipelines in our analysis is that all pipelines are unbiased and positively correlated (i.e., they all estimate the “same” effect). This can, however, only sometimes be guaranteed. For example, in the well-known NARPS analysis (Botvinik-Nezer et al., 2020), the effect detected by different groups was, for a few cases, in the opposite direction (which is not typical of multiverse analyses). Alternative approaches, for example, STAPLE (Dafflon et al., 2022), could be used to reduce the bias of the pooled estimators by assuming a majority of unbiased pipelines or identifying clusters of pipelines and pooling pipeline-specific estimates within clusters. Finally, our proposed framework currently only considers a single region at a time and not multiple regions of interest jointly. Previous work has shown that optimal preprocessing strategies for single regions exist (Nørgaard, Ganz, et al., 2019), but also that multiple regions of interest are correlated within-subject (Eklund et al., 2016; Ganz et al., 2020). Therefore, future work is still needed to further jointly investigate and summarize the complex and nested dependence structure across multiple regions, and how it interacts with multiple preprocessing strategies.

In this work, we have developed a statistical sensitivity analysis that can quantify the impact of different preprocessing choices on subsequent statistical analyses. As has been reported in previous studies, we observe that the influence of preprocessing pipelines on subsequent statistical analysis can be quite large. Hence, we provide tools for statistical analyses that can aggregate multiple analyses of the same data. We introduced a new pooling estimator across pipelines (⁠$Ψ^constrained GLS$⁠) but also discuss how to estimate the proportion of pipelines with evidence for an effect and test hypotheses across pipelines, such as “no effect across all pipelines” or “at least one pipeline with no effect”. Code to reproduce all simulations and figures, including a step-by-step tutorial, is available at https://github.com/openneuropet/multiverse_tools/tree/main.

The authors ensure that the work submitted is original and all authors are credited appropriately and have agreed to the journal’s open access, ethical, and data-sharing policies. All clinical data used in this study was ethically acquired following relevant ethical and legal guidelines. There are no conflicts of interest with respect to the work presented in this article.

The script to re-create the in silico data is available on GitHub (https://github.com/openneuropet/multiverse_tools/tree/main/sensitivity_analysis-simulation). The data for the real-world example originate from a clinical study by Frokjaer et al. (2015). The study was registered and approved by the ethics committee for the capital region of Copenhagen (protocol-ID: H-2-2010-108) and registered as a clinical trial: www.clinicaltrials.gov under the trial ID NCT02661789. All subjects provided written informed consent prior to participation, in accordance with The Declaration of Helsinki II. Data are available from the CIMBI database (https://nru.dk/index.php?view=category&id=128) upon request.

Code implementing the proposed methods is available on CRAN in the LMMstar package (https://cran.r-project.org/web/packages/LMMstar/index.html). Code for the simulation study and the application is available on GitHub (https://github.com/openneuropet/multiverse_tools/tree/main/sensitivity_analysis-simulation).

All authors participated in the conceptualization and first drafting of the manuscript. M.N. performed the data curation. B.O. developed the methodology and wrote the software. M.N., C.P., and M.G.B. participated in the validation and formal analysis. B.O. and M.G.B. added visualizations. All authors read and approved the final manuscript.

We thankfully acknowledge funding through the BRAIN Initiative grant #1R24MH120004-01A1 and the Novo Nordisk Foundation grant #N20OC0063277, as well as the Lundbeck foundation grant BrainDrugs #R279-2018-1145.

There are no competing interests with respect to the presented work.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00523.