Corrections for multiple comparisons in voxel-based lesion-symptom mapping

1. Introduction

Identifying relationships between location of brain damage and cognitive deficits is a foundational method in cognitive neuroscience, tracing its history at least to the behavioral neurologists of the mid-19th century (e.g., Lichtheim, 1885). Those early studies were based on individual case studies and, as data accumulated, researchers used lesion overlays to identify the locations where damage consistently produced deficits of interest. Recent advances in neuroimaging technology have allowed much finer-grained analyses at the level of individual voxels (Bates et al., 2003; Rorden and Karnath, 2004). In voxel-based lesion-symptom mapping (VLSM), an association between deficit severity and lesion status (lesioned vs. not lesioned) is tested in each voxel, producing a statistical map of the strength of relationship between lesion status and deficit. However, this map is the result of individual tests across tens or even hundreds of thousands of voxels.

The large number of tests involved in analysis of neuroimaging data requires some kind of statistical correction for multiple comparisons. Several strategies have been proposed, often by adaptation from analysis of functional neuroimaging (e.g., fMRI). One standard strategy of correction for multple comparisons is to control voxel-level family-wise error rate (FWER), which is the probability of making one or more false positive (Type 1) errors among the entire set of tests. The Bonferroni correction is a classic FWER correction method, though it is generally considered overly conservative for neuroimaging data. An alternative, non-parametric, approach to FWER correction is to use permutations of the observed data to build a null distribution of test statistics and compare the observed test statistic against that null distribution to determine the likelihood of observing the result if the null hypothesis were true. Because they are based on permuting the real data, this approach has the important advantage of not making assumptions about the distributions of scores or test statistics – assumptions that are likely to be violated by skewed distributions of behavioral deficit (symptom) scores and by the spatial contiguity inherent to stroke lesions. Building a null distribution based on permutations of real data offers a rather literal way to compute p-values: the p-value is literally the probability of observing a particular outcome if there were no relationship between the behavioral scores and lesion patterns (i.e., random permutations). A null distribution of the test statistic can be built based on permutations of real data and used to reject voxels where the true analysis does not sufficiently differ (e.g., p > 0.01) from the permutation-based null distribution to warrant rejecting the null hypothesis (e.g., Kimberg et al., 2007; Rorden et al., 2007).

Controlling the probability of making one or more false positive errors is based on the idea that each test is critically related to the researcher’s interpretation or inference, thus, a single false positive could potentially undermine the inference and needs to be controlled. Voxel-level FWER does not align with VLSM interpretation, which never depends on a single voxel. The misalignment between standard FWER correction and VLSM interpretation makes standard voxel-level FWER correction unnecessarily conservative: if no inferences are made based on a single voxel, then a single false positive voxel cannot be responsible for an invalid inference about lesion-symptom relations.

The general approach of permutation-based correction for multiple comparisons can be implemented in many different ways, depending on what aspect of the results is to be controlled. Instead of controlling the rate of a single false positive voxel, permutation-based FWER can be used to set a minimum cluster size, thus controlling the rate of a single false positive cluster of voxels. Setting a minimum cluster size is a common “clean-up” step in neuroimaging data analysis; the addition of a principled strategy for selecting the minimum cluster size is the critical component that turns this into a statistical correction method. The permutation-based strategy is to set a voxel-wise cluster-forming threshold (e.g., p < 0.0001), then use permutations to determine the null distribution of cluster sizes that pass this threshold, and use that distribution to set a minimum cluster size (e.g., Nichols and Holmes, 2002). This approach is appealing because stroke lesions are inherently contiguous and VLSM interpretation tends to focus on clusters.

Another technique for controlling the rate of false positives is False Discovery Rate (FDR), which quantifies the proportion of above-threshold results that can be expected to be false positives (Genovese et al., 2002). That is, at FDR threshold q = 0.05, 5% of above-threshold voxels are expected to be false positives, which is likely to be substantially more than one voxel but not likely to affect interpretation of the overall patten (for a clear description see Bennett et al., 2009). FDR is widely used for analysis of functional neuroimaging data and VLSM, however, we have encountered informal criticism that FDR is inappropriate for VLSM. FDR is certainly less conservative than FWER (see also Rorden et al., 2007), but that is true by design – FDR is designed to allow a small percentage of false positive voxels, whereas FWER aims to make it unlikely that there is even a single false positive voxel – and we are not aware of any published analysis showing that FDR incorrectly quantifies the rate of false positive voxels in VLSM.

Correction for multiple comparisons is an attempt to manage variability, but it cannot remove all of the noise and leave all of the signal. Either some noise will get left behind or some of the signal will be removed. That is, there is an inherent trade-off between false positives and false negatives; incorrectly generalizing a result and overlooking a generalization that is warranted. By convention, data analysis requires setting a threshold to identify results that warrant rejection of the null hypothesis. There is a substantial price associated with adopting the conservative position that the probability of even a single false positive voxel needs to be controlled: VLSM analysis is based on a single data point per participant (each participant only has one lesion and only one deficit profile) and sample sizes are often limited by the practical challenges of recruiting and testing large numbers of participants with the targeted neurogenic deficits. This price is further exacerbated by publication bias: studies that meet the statistical threshold may be published, but studies that fall short are relegated to the “file drawer”, leaving a biased scientific literature. Publication bias also encourages various forms of “p-hacking” or “researcher degrees of freedom”, in which researchers try alternative analysis strategies (excluding certain “outlier” participants, transforming scores, etc.) until they find one that surpasses the statistical threshold. The result is a report that appears to use rigorous statistical methods, but the actual rate of false positives far exceeds the nominal p-value (e.g., Simmons et al., 2011; Nosek et al., 2012; Gelman and Loken, 2014). In addition to statistical soundness, the analytical strategy should allow researchers to transparently report their observations and the strength of the evidence that supports their conclusions.

The present study investigated two permutation-based methods of correcting for multiple comparisons in VLSM. The next section describes our investigation of using permutations to determine a minimum cluster size. Our analyses found that this approach produces consistent spill-over into neighboring regions (i.e., the identified region extends well beyong the boundaries of the true lesion-symptom relation), making it not well-suited to identifying brain-behavior relationships. Although some spill-over is a necessary consequence of high spatial correlations among neighboring voxels inherent to stroke lesion data, cluster-based correction was substantially more susceptible to this problem than standard voxel-level FWER correction. The subsequent section describes a generalization of the permutation-based FWER correction approach that captures some of the inferential advantages of FDR and cluster-based correction without making parameteric assumptions about the data. This approach makes it possible to balance control of false postives against risk of false negatives, and to transparently report results in a way that allows others to evaluate the evidence. We also compare this approach with the parametric FDR method and describe conditions under which FDR may produce misleading results. The final section of this report summarizes our findings and conclusions, and discusses future directions.

2. Minimum Cluster Size

Using permutations to determine a minimum cluster size proceeds as follows: (1) permute behavioral data and conduct VLSM analysis, (2) apply a pre-set cluster-forming threshold for each voxel (e.g., p < 0.0001), (3) compute size of largest supra-threshold voxel cluster, (4) repeat steps 1–3 many times to build up a null distribution of supra-threshold cluster sizes (e.g., Nichols and Holmes, 2002). This null distribution is the distribution of largest cluster sizes that are observed when there is no relationship between deficit scores and lesion location. Clusters from the original (true) VLSM analysis that are larger than 95% of the null distribution of cluster sizes are taken to reflect true lesion-symptom associations (for examples of application to VLSM see Pillay et al., 2014, 2017; Binder et al., 2016; Mirman et al., 2015a). This strategy involves two separate thresholds: the first is a pre-set voxel-level p-threshold; the second is a permutation-based cluster size threshold. This strategy has two intuitively appealing properties. First, strokes and other neurological disorders tend to produce spatially contiguous lesions, resulting in high spatial correlations between the lesion status of neighboring voxels. Using permutation to determine a null distribution of cluster sizes intuitively controls for this spatial correlation and produces a minimum cluster size threshold that should not be observed by chance. Second, it is typical for interpretation of VLSM (and other neuroimaging) results to focus on clusters, so correcting at the cluster level (rather than the voxel level) makes this statistical strategy more closely aligned with the interpretion strategy. In the following analyses we examine this permutation-based cluster-size correction strategy for detecting true lesion-symptom relations.

3. Continuous Permutation-Based FWER

The standard permutation-based FWER correction method proceeds as follows: (1) permute behavioral data and conduct VLSM analysis, (2) identify the maximal test statistic (typically, the most extreme t-value), (3) repeat steps 1 and 2 many times to build a null distribution of maximal t-values, (4) compute the n-th percentile of that null distribution to determine a threshold for the test statistic, which corresponds to n% of the permutations having 0 voxels that exceed this threshold (Rorden et al., 2007). A typical value of n is 95, which produces a FWER-corrected p < 0.05: less than 5% of the permutations had even a single voxel that exceeded this t-threshold.

This approach controls the rate of single-voxel false positives, but it can be generalized to multi-voxel false positives by focusing on the v-th most extreme test statistic⁴. The standard strategy is the special case when v = 1, thus using the most extreme voxel-wise test statistic from each permutation, and controlling the rate of 1 false positive voxel. If, for example, v = 10, one would similarly use the 10th most extreme voxel-wise test statistic from each permutation, and control the rate of up to 10 false positive voxels. An example is shown in Figure 4 where the left panel shows the sorted t-values from the first 10 permutations and the right panel shows permutation-based t-value distributions and 95% thresholds at v = 1 (red), v = 10 (green), v = 100 (blue), and v = 1000 (purple). A single set of permutations produced the set of possible 95% thresholds at different critical voxel (v) values, reflecting the expected number of false positive voxels at the corresponding t threshold. Not surprisingly, that t threshold decreases as the v value increases, but the relationship remains the same as the standard FWER case: the v value specifies the maximum number of voxels that exceeded the corresponding t-threshold in 95% of the permutations. That is, a reasonable upper bound on number of false positive voxels at that t threshold. We refer to this extension of the standard permutation-based FWER as continuous permutation-based FWER because it uses the same permutation-based FWER strategy but allows values of v > 1.

Since interpretation of VLSM results typically relies on a large set of voxels, that interpretation is unlikely to be affected by a small (but > 1) number of false positive voxels. This is not to say that, for example, v = 10 should be adopted instead of v = 1; rather, this generalization of FWER allows investigators to assess the strength of their evidence in a flexible way and to report that assessment in a way that allows readers (and reviewers) to evaluate the claims.

Continuous permutation-based FWER is somewhat similar to the false discovery rate (FDR) approach in that both allow multiple false positive voxels and quantify that rate of false positives. However, the two methods differ in two important ways. First, FDR is designed to control the proportion of supra-threshold voxels that are expected to be false positives (this proportion is usually reported as the q-value), whereas the FWER approach quantifies the number of possible false positive voxels, which may be a high or low proportion of supra-threshold voxels. Second, continuous FWER is permutation-based, which means that (unlike FDR) it makes no assumptions about distributions of data or test statistics. The latter property is important because lesion data may violate the assumptions of FDR severely enough to make FDR unreliable for VLSM. Here we report results from application of this continuous permutation-based FWER approach in several contexts. In addition, we used the permutation data to evaluate whether the nominal FDR q-value correctly quantifies the proportion of supra-threshold voxels that are expected to be false positives.

4. General Discussion and Conclusions

Permutation-based FWER correction is the current “gold standard” correction for multiple comparisons in VLSM. The standard permutation-based FWER strategy is to build a null distribution using only most extreme voxel-wise test statistic from each permutation. The main weakness of this approach is that it aims to control the occurrence of even a single false positive voxel, which is not the scale at which VLSM results are interpreted. We described an extension - continuous permutation-based FWER - which better aligns with VLSM interpetation and allows researchers a more flexible balance between false positives and false negatives. Continuous FWER uses the v-th most extreme voxel-wise test statistic, so the standard approach is the special case where v = 1, but other values v > 1 may be used as appropriate for a particular data set and hypothesis. This provides a principled way for researchers to flexibily set the upper limit of how many false positive voxels are allowed and to transparently report this limit along with their results, so readers can also evaluate the evidence. Since single voxels rarely (if ever) affect the interpretation of VLSM results, this flexibility lets researchers align their statistical method with their interpretations of the results.

In addition to continuous FWER, we examined two other methods of correction, but those results were not encouraging. Using permutations to set a minimum cluster size tended to produce clusters that extended well beyond the correct region. This was partly due to a true correlation between damage in adjacent regions and damage in the correct region (i.e., spatial coherence); however, since the other correction methods seemed less susceptible to this spill-over problem, spatial auto-correlation does not appear to be the full explanation. Instead, we suspect this spill-over occurred because weak or noisy effects in adjacent voxels were incorporated into true clusters, with the unfortunate consequence of blurring the boundary of the true symptom-related region. Cluster-based correction appeared to be effective at controlling the occurrence of false positive clusters, but the spill-over effect poses a problem for identifying brain-behavior relationships. The spill-over effect may have been further exacerbated by the strong relationship between our simulated deficit scores (percent damage) and their corresponding brain regions. The false discovery rate (FDR) approach is inferentially similar to continuous FWER: it aims to quantify the rate of false positive voxels. FDR performed quite well for larger samples with strong lesion-symptom relations, but consistently underestimated the rate of false positive voxels when the sample sizes were smaller and in a real data case (where the lesion-symptom relation is likely to be weaker). This suggests that researchers should be wary of using FDR in conventional VLSM analyses.

There is an inherent trade-off between false positives and false negatives: striving to eliminate false positives will necessarily result in missing many true effects, but generalizing from every observation will necessarily produce some incorrect inferences. Setting arbitrary thresholds of statistical significance makes evidence appear more dichotmous than it really is; statistical thresholds encourage binary thinking in which an effect is either significant or non-existent. This dichotomy is further exacerbated by publication bias because weaker, not statistically significant results are simply not published. This reifies the sense that effects that do not pass the significance threshold are non-existent, leading to a biased scientific literature and undermining evidence accumulation. Balancing false positives and false negatives is particularly challenging in VLSM, where participant recruitment and testing is difficult and relatively expensive, and samples are generally large relative to typical research in neuropsychology and cognitive neuroscience. A typical VLSM research project might require a long period of expensive data collection to reach a reasonable sample size of, say, N=50. If analysis of that data set produced results just short of the standard FWER v = 1 statistical threshold, the researchers would be left with an unpublishable result. Substantially increasing the sample size is likely to be impractical (and perhaps impossible) as it would require another long, labor-intensive, and expensive data collection effort. Addressing this challenge requires a statistical correction method that allows researchers to flexibly balance false positives and false negatives and to report how they struck that balance in a transparent fashion so that readers can interpret the evidence. The v threshold plays this role in the continuous permutation-based FWER correction method: v is the expected upper limit of false positive voxels, which can be adjusted to suit the researchers’ hypotheses and reported for readers to use in their evaluation. Making the full statistical maps easily available to other researchers through a repository would further support evidence accumulation through re-analysis and meta-analysis.

We have deliberately avoided recommending a specific v threshold to be used in continuous FWER correction, or how many supra-threshold voxels (effective q) should be considered “significant”. Such thresholds are fundamentally arbitrary – there is nothing qualitatively different between v = 10 and v = 11, just as there is nothing qualitatively different between p = 0.04 and p = 0.06. Setting significance thresholds contributes to mis-interpretation of these continuous statistics and results that do not pass the threshold tend not to be published, which creates an incentive for researchers to use (undisclosed) variations in data pre-processing and analysis methods to achieve a “significant” result (for more discussion see Amrhein et al., 2017). Thus, in lieu of recommending specific choices of hard thresholds, we offer guidelines for effective use of continuous FWER. First, the v threshold should match the neural specificity of the hypothesis under investigation. Neurally precise hypotheses (e.g., anterior vs. posterior portion of the inferior frontal gyrus) may require a small v threshold; in contrast, a high v threshold may be sufficient for evaluating a neuroanatomically broader hypothesis (e.g., role of anterior temporal lobe vs. the temporo-parietal cortex). It may also be useful to consider the overall number of voxels that are entered into the analysis, especially for VLSM analyses that are constrained to specific regions of interest. Second, the number of supra-threshold voxels should be substantially larger than v (the effective q, which should be fairly small, may serve as a useful summary statistic). For example, observing 101 supra-threshold voxels at v = 100 would not license strong conclusions because almost that many could be false positives. That is, the conclusions drawn by investigators should reflect both the number of supra-threshold voxels relative to the v threshold (effective q) and the neural precision of their hypothesis. These general guidelines may not be sufficient for selection of a single appropriate v threshold, so it may be useful to test and report results at multiple v thresholds as a way to assess the robustness of the effects. Continuous FWER is an extension of standard permutation-based FWER that enables the researcher to explore and report the strength of the effects by testing different v thresholds. Further investigation by independent groups will help to refine and set community standards for how v thresholds should be set and explored. We have made available a basic implemenation⁵ of continuous FWER and one that is integrated into a user-friendly lesion-symptom mapping package⁶, both of which make it easy to test multiple different v thresholds in order to provide a more complete understanding of the strength of the evidence.

It may be possible to further improve the correction methods described here by considering unique patches (where lesion status is perfectly correlated across participants) rather than individual voxels (Kimberg et al., 2007). Our (unsystematic) comparisons suggest that, in a typical data set, the number of unique patches may be an order of magnitude smaller than the number of voxels. Such a vast reduction in the number of (redundant) tests would have substantial consequences for these methods, as well as for the processing time required to compute them. More generally, including voxel neighborhood information as part of voxel-level corrections may lead to even more effective algorithms. In addition, it may be possible to develop algorithms to automatically establish an optimal v-threshold by using spatial correlations in the data (e.g., number of patches, neighborhood correlations) as an indication of the amount of continuous FWER correction required for the specific dataset under investigation. Finally, recent development of multivariate lesion-symptom mapping methods (Zhang et al., 2014; Pustina et al., 2017), which evaluate lesion-symptom relationships across all voxels simultaneously, provide a better method for studying brain-behavior relationships and mitigate the need for multiple comparisons correction. Such methods are not in wide use yet, but they offer a promising alternative approach.

In addition to optimizing corrections for multiple comparisons, getting the most out of lesion-symptom studies requires optimizing the test statistics themselves. Recent analysis show that correcting for overall lesion volume and only testing voxels with sufficient lesion involvement is important for avoiding mis-localization (Sperber and Karnath, 2017). For functional neuroimaging data, threshold-free cluster enhancement (TFCE; Smith and Nichols, 2009) has been shown to improve voxel-wise signal-to-noise ratio by incorporating local support infromation from nearby voxels. This method considers all possible cluster-forming thresholds, so it can caputre support from small strong clusters as well as weaker diffuse clusters. Once the voxel-wise TFCE statistic has been calculated, it still requires correction for multiple comparisons. Because they operate at different steps in the analysis process, calculating voxel-wise test statistics using TFCE, then correcting those test statistics for multiple comparisons using continuous FWER correction may jointly enhance power in VLSM analysis⁷.

VLSM is important for basic and translational human neuroscience. Analysis of lesion-symptom relations has been at the core of cognitive neuroscience research since the mid-19th century and remains critical to the field. Lesions that produce chronic deficits in a particular domain or task provide the strongest evidence that the damaged neural structures were critical for that domain or task. This method is an important complement to functional neuroimaging in neurologically-intact populations, but the differences in data collection challenges create somewhat different statistical demands. Cognitive neuroscience has tremendous potential for stimulating development of new and improved diagnosis, treatment, rehabilitation, and education strategies. That potential cannot be realized without testing the affected populations. VLSM offers a unique opportunity for research that answers fundamental questions about the neural basis of cognition while addressing the real-world problem of understanding neurogenic cognitive deficits. Robust, flexible, and transparent statistical methods play an important role in maximizing the impact of VLSM research.