An empirical comparison of univariate versus multivariate methods for the analysis of brain–behavior mapping

Lesion symptom mapping (LSM) tools are used on brain injury data to identify the neural structures critical for a given behavior or symptom. Univariate lesion symptom mapping (ULSM) methods provide statistical comparisons of behavioral test scores in patients with and without a lesion on a voxel by voxel basis. More recently, multivariate lesion symptom mapping (MLSM) methods have been developed that consider the effects of all lesioned voxels in one model simultaneously. In the current study, we provide a much‐needed systematic comparison of several ULSM and MLSM methods, using both synthetic and real data to identify the potential strengths and weaknesses of both approaches. We tested the spatial precision of each LSM method for both single and dual (network type) anatomical target simulations across anatomical target location, sample size, noise level, and lesion smoothing. Additionally, we performed false positive simulations to identify the characteristics associated with each method's spurious findings. Simulations showed no clear superiority of either ULSM or MLSM methods overall, but rather highlighted specific advantages of different methods. No single method produced a thresholded LSM map that exclusively delineated brain regions associated with the target behavior. Thus, different LSM methods are indicated, depending on the particular study design, specific hypotheses, and sample size. Overall, we recommend the use of both ULSM and MLSM methods in tandem to enhance confidence in the results: Brain foci identified as significant across both types of methods are unlikely to be spurious and can be confidently reported as robust results.


Abstract
Lesion symptom mapping (LSM) tools are used on brain injury data to identify the neural structures critical for a given behavior or symptom. Univariate lesion symptom mapping (ULSM) methods provide statistical comparisons of behavioral test scores in patients with and without a lesion on a voxel by voxel basis. More recently, multivariate lesion symptom mapping (MLSM) methods have been developed that consider the effects of all lesioned voxels in one model simultaneously. In the current study, we provide a much-needed systematic comparison of several ULSM and MLSM methods, using both synthetic and real data to identify the potential strengths and weaknesses of both approaches. We tested the spatial precision of each LSM method for both single and dual (network type) anatomical target simulations across anatomical target location, sample size, noise level, and lesion smoothing. Additionally, we performed false positive simulations to identify the characteristics associated with each method's spurious findings. Simulations showed no clear superiority of either ULSM or MLSM methods overall, but rather highlighted specific advantages of different methods. No single method produced a thresholded LSM map that exclusively delineated brain regions associated with the target behavior. Thus, different LSM methods are indicated, depending on the particular study design, specific hypotheses, and sample size. Overall, we recommend the use of both ULSM and MLSM methods in tandem to enhance confidence in the results: Brain foci identified as significant across both types of methods are unlikely to be spurious and can be confidently reported as robust results.

K E Y W O R D S
aphasia, brain-behavior relationships, language, lesion symptom mapping, multivariate, stroke, VLSM Maria V. Ivanova and Timothy J. Herron contributed equally to this study.

| Lesion symptom mapping
Throughout the 19th and much of the 20th century, systematic clinical observations of neurologic patients along with postmortem autopsy remained the main method for establishing brain correlates of cognitive functioning (Damasio & Damasio, 1989;Dronkers, Ivanova, & Baldo, 2017;Luria, 1980). The advent of modern neuroimaging methods in the 1970s greatly enhanced the ability to determine neural foundations of cognition, as the actual lesion site could be Voxels that show significant differences for a particular behavior or symptom are inferred to be critical for the behavior under examination. ULSM methods complement functional neuroimaging studies in healthy participants, by testing the necessity of particular brain areas for a particular behavior, thereby demonstrating the crucial causal link in brain-behavior relationships (Bates et al., 2003;Karnath, Sperber, & Rorden, 2018;Rorden, Karnath, & Bonilha, 2007;Vaidya, Pujara, Petrides, Murray, & Fellows, 2019). Contemporary ULSM methods provide a fundamental shift in broadening our understanding of brain-behavior relationships, both confirming (Baldo, Arevalo, Patterson, & Dronkers, 2013) and challenging previously held beliefs about key neural structures for different cognitive functions (Baldo et al., 2018;Dronkers, Wilkins, Van Valin, Redfern, & Jaeger, 2004;Ivanova et al., 2018;Mirman, Chen, et al., 2015).
More recently, new multivariate lesion symptom mapping (MLSM) methods have been developed as an alternative to ULSM.
The principal difference between ULSM and MLSM methods is that MLSM considers the entirety of all lesion patterns in one model simultaneously. This is in contrast to the parallel, independent analysis of lesion patterns on a voxel by voxel basis performed with ULSM models. While some papers have argued that MLSM methods should be superior to ULSM methods (DeMarco & Turkeltaub, 2018;Mah, Husain, Rees, & Nachev, 2014;Pustina, Avants, Faseyitan, Medaglia, & Coslett, 2018;Zhang, Kimberg, Coslett, Schwartz, & Wang, 2014), many of these arguments have been presented theoretically without rigorously comparing the two approaches (see also .
Below, we focus on two critical properties of LSM analyses that are foundational to the validity of the method: spatial accuracy and ability to detect networks. With regard to these two aspects, we first review issues that impact brain-behavior inferences made with both ULSM and MLSM methods and then report a comprehensive empirical evaluation of these two different approaches to lesionbehavior mapping.
1.2 | Spatial accuracy 1.2.1 | Issues affecting spatial accuracy The original ULSM method was superior in terms of its spatial accuracy to lesion overlays and lesion subtraction analyses as it provided a quantifiable statistical approach to capturing the continuous nature of behavioral data in relation to lesion site. However, recently multiple concerns over the spatial accuracy of the resulting LSM maps have been raised.
One issue that affects both ULSM and MLSM methods is that lesion distributions in stroke (the most frequently studied etiology with LSM techniques) are influenced by the vascular anatomy and are thus nonrandomly distributed in the brain with certain areas being more likely to be lesioned than others (Mah et al., 2014;Phan, Donnan, Wright, & Reutens, 2005;Sperber & Karnath, 2016;Xu, Jha, & Nachev, 2018). This nonrandom distribution of lesions impacts LSM analyses in the following ways. First, it limits analysis of certain brain areas that are rarely affected in stroke (e.g., the temporal pole).
Second, neighboring voxels have a higher probability of being lesioned together, as strokes never affect just one voxel. The ULSM approach is potentially susceptible to this spatial autocorrelation, because it assumes independence of lesioned voxels throughout the brain, as thousands of independent tests (nonparametric, t tests, or linear regressions) are carried out serially in the affected voxels. While independence of tests is not an assumption of MLSM methods per se (since only one multivariate model incorporating all the lesion patterns is tested), lack of sufficient spatial distinction is an issue. In other words, if two voxels are always either damaged together or always spared, it is not possible to differentiate their unique contribution to the observed deficits with any LSM method. A third related concern that affects both ULSM and MLSM methods is differential statistical power across voxels/regions of the brain. For example, a voxel in which 50% of patients have a lesion has more power than a voxel where only 10% of patients have a lesion . In general, given both the nonrandom nature of lesions and the inability to predict their specific pattern in a given study, it is hard to estimate statistical power in advance for any LSM method. Generally, one can only perform post hoc power analyses to determine the amount of power in different brain regions. Even with very large samples, statistical power can be low everywhere in the brain, for example, if lesions are small and voxels are only affected by a small proportion of all lesions. Finally, stroke lesions are also typically larger than the functional anatomical targets that LSM analyses attempt to uncover, thereby limiting the spatial resolution of the analysis.
Cumulative effects of nonrandom lesion distribution, autocorrelation across voxels, and differential power distribution can potentially lead to distortion in spatial localization of critical regions. Significant clusters are often "diverted" toward the most frequently damaged regions, which are regularly impacted together with the true correlates of a cognitive function (Inoue, Madhyastha, Rudrauf, Mehta, & Grabowski, 2014;Xu et al., 2018) and potentially along the brain's vasculature (Mah et al., 2014;. The open question is, to what degree do these lesion-anatomical biases impact different ULSM and MLSM methods (Sperber, 2020), and how are they ameliorated by sample size, method choices, additional corrections, and interpretation? Currently, very few studies have investigated these biases systematically and compared them across different LSM methods.

| Empirical studies investigating spatial accuracy of LSM methods
One of the original papers to raise awareness about spatial distortion in LSM by Mah et al. (2014) suggested that ULSM analyses mislocalized foci by an average of 16 mm. However, their model did not include lesion volume as a covariate in their analysis. The importance of using lesion volume as a nuisance covariate in LSM has been a standard recommendation for several years (Baldo, Wilson, & Dronkers, 2012;De Haan & Karnath, 2018;DeMarco & Turkeltaub, 2018;Price, Hope, & Seghier, 2017;Sperber & Karnath, 2017). In addition, Mah et al. used a minimum lesion load per voxel of <1% in their ULSM analyses, which is far below the standard recommendation of 5-10% (Baldo et al., 2012). Moreover, the displacement maps using synthetic data in Mah et al. showed single voxels (i.e., a single voxel leading to a specific deficit), which is an oversimplified and exclusively theoretical case that does not occur naturally. Furthermore, when damage to an anatomical region was used as a synthetic behavioral score in their study, the score was binarized rather than continuous, likely further reducing spatial resolution. Finally, Mah et al. did not explore spatial bias for MLSM, so within their study, it was not possible to directly compare the spatial displacement between ULSM and MLSM methods.
In another simulation study critiquing accuracy of ULSM (Inoue et al., 2014), lesion volume was included as a covariate, but the authors again used binarized synthetic behavioral scores (a deficit was indicated when 20% of voxels in the target parcel were damaged) and did not apply a minimal lesion load threshold. Also, the results in this study were predominantly analyzed with false discovery rate (FDR)based thresholding. This method of correction for multiple comparisons has been discontinued for some time in the ULSM literature, as it frequently leads to an increase in false positives (Baldo et al., 2012;Kimberg, Coslett, & Schwartz, 2007;Mirman et al., 2018). Also, accuracy of mapping was not systematically explored across sample sizes.
Finally, the lesion data for this study came from highly heterogenous etiologies (stroke, traumatic brain injury, encephalitis), contrary to standard recommendations for any LSM study (De Haan & Karnath, 2018).
Most prominently, Sperber and Karnath (2017) empirically demonstrated that ensuring a sufficiently large minimal lesion load threshold as well as including lesion volume as a covariate have significant additive effects on improving spatial precision of results, although not entirely removing spatial bias. In their study, spatial displacement was calculated for single voxels, with displacement of larger clusters expected to be smaller. Furthermore, the lesion volume correction to enhance accuracy of localization has been strongly recommended for at least some MLSM approaches (DeMarco & Turkeltaub, 2018), again highlighting that MLSM methods are not immune to these types of spatial biases. In another simulation study, Sperber, Wiesen, Goldenberg, and Karnath (2019) and  showed that a common support vector regression (SVR)-based MLSM method was also susceptible to mislocalization along the brain's vasculature, even after applying a correction for lesion volume, and that this displacement error was actually higher than that observed for a ULSM method. However, since displacement was determined for single voxels in a single axial slice, these spatial biases require further exploration to fully understand their impact on LSM results with real behavioral data.
The most comprehensive simulation study to date by Pustina et al. (2018) showed that even one of the most advanced MLSM algorithms, sparse canonical correlation analysis for neuroimaging (SCCAN), exhibited spatial bias in the results. Here, the superiority of SCCAN using synthetic data was consistently demonstrated, but only when compared to the univariate analyses with inappropriate FDRbased thresholding. Modern ULSM methods instead use a conservative, permutation-based familywise error rate (FWER) correction, a nonparametric resampling approach to significance testing, which sets the overall probability rate of false positives across all of the results, while making almost no assumptions about the underlying data distributions (Hayasaka & Nichols, 2003;Nichols & Holmes, 2001).
Permutation-based FWER provides the most stringent and robust form of correction for multiple comparisons, providing an optimal balance between false positives and false negatives (see Kimberg et al., 2007;Mirman et al., 2018). Accordingly, in the same paper, the ULSM results obtained with this more appropriate thresholding using permutation-based and Bonferoni FWER corrections, were comparable to SCCAN results across a number of spatial indices . Moreover, Pustina et al. (2018) did not include lesion size as a covariate in the ULSM analysis, running counter to standard recommendations for ULSM and potentially biasing the comparison (Baldo et al., 2012;Sperber & Karnath, 2017). Furthermore, limited spatial metrics were used as measures of accurate mapping in comparing LSM methods, and most of these metrics produced similar levels of performance for all methods tested. For example, while the dice index (measure of overlap between two regions) was shown to be significantly higher for SCCAN compared to a nonparametric Brunner-Munzel version of ULSM, values were very low in nearly all cases with every method (predominantly <0.5 and often <0.2) rendering the statistical advantage uninformative. Also, results of statistical comparisons across different sample sizes for other spatial metrics were not provided (see Figure 4, p. 161, Pustina et al., 2018).
In short, the degree to which spatial bias affects ULSM versus MLSM methods has not yet been systematically and rigorously tested across a wide range of LSM methods with a wide range of metrics of spatial accuracy, while implementing best practices such as lesion volume control, minimum lesion load threshold, and proper correction for multiple comparisons.

| Detection of networks
Another issue in LSM is the ability to detect complex relationships and functional dependencies in the data (i.e., networks). Given that most complex cognitive functions are supported by a number of regions working together in a coordinated fashion, it is pivotal that LSM methods are able to uncover multiple regions underpinning the target behavior.
Some papers have argued that MLSM should be better than ULSM at detecting multifocal relationships between lesion location and behavioral deficits (i.e., when damage to multiple areas leads to a specific behavioral deficit), as MSLM takes into account all the voxels simultaneously in a single model (DeMarco & Turkeltaub, 2018;Mah et al., 2014;Pustina et al., 2018;Sperber, 2020;Zhang et al., 2014). However, the empirical evidence favoring the superior ability of MLSM methods to identify networks remains inconclusive. This is in part due to the strong regularization (e.g., sparse vs. dense solutions) and additional assumptions (e.g., restriction on possible locations of solutions) required in order to solve a single, massively underdetermined, multivariate system of equations (typically with 1 patient per 100 or 1,000 voxels). Mah et al. (2014) claimed that MLSM resulted in higher sensitivity and specificity compared to ULSM in detecting a two-parcel fragile network (when the synthetic score was based on the maximal lesion load among a set of anatomical regions). However, as described above, synthetic behavioral scores were binarized for their analysis, the statistical threshold used was not specified, and there was no quantification of the differences in spatial bias between ULSM and MLSM. Pustina et al. (2018) showed an advantage of MLSM over ULSM in detecting an extended network ("AND" rule; when the synthetic score was based on the average lesion load among a set of anatomical regions) consisting of three parcels. However, there was no significant advantage of the MLSM over ULSM in detecting other types of two-and three-parcel networks when a proper FWER correction was included.
Findings with real behavioral data also remain mixed. For instance, in one study MLSM methods were able to detect a brain network underlying apraxia of pantomime, while ULSM could not (Sperber, Wiesen, Goldenberg, & Karnath, 2019). However, a number of studies have repeatedly shown that ULSM methods are able to detect spatially distinct regions in a network (Akinina et al., 2019;Baldo et al., 2018;Gajardo-Vidal et al., 2018;Mirman, Chen, et al., 2015). With real behavioral data, however, there is no ground truth with which to compare the results, so it remains possible that the analyses should have uncovered even more relevant regions.
Thus, further empirical evidence is needed to show what specific measurable advantage MLSM has over ULSM in detecting multifocal behavioral determinants.

| Aims of the current study
To summarize, there are a number of theoretical concerns for both ULSM and MLSM methods. Some of these concerns raised originally with respect to ULSM (e.g., spatial bias and autocorrelation, differential statistical power) are actually concerns for both ULSM and MLSM methods and require further elucidation with respect to both approaches. Moreover, efficient controls already exist for both ULSM and MLSM methods that can be implemented to minimize the biasing effect of lesion physiology (e.g., lesion size correction, minimum lesion load threshold). The theoretical concerns about ULSM methods being less able to detect networks of brain regions (as opposed to a single target region) have not been systematically confirmed. Properties of new LSM models require further delineation and comparative evaluation in order to assess the mapping power and accuracy under varying conditions. Further, certain factors that can potentially impact accuracy of analysis such as lesion mask smoothing and behavioral noise levels have not been properly addressed in previous papers comparing different LSM methods. To date, the comparisons of ULSM and MLSM in the literature have been limited and when they are contrasted, a sub-standard version of ULSM is often implemented without proper correction, leading to an unfavorable impression of ULSM (Inoue et al., 2014;Pustina et al., 2018;Zhang et al., 2014). Finally, neither ULSM nor MLSM methods have been properly explored with respect to the incidence of false positive results.
The current paper aimed to address these gaps in the LSM literature and provide a comprehensive appraisal of several versions of ULSM and MLSM methods with a large stroke lesion dataset, using both synthetic and real behavioral data, across a range of relevant parameters. Synthetic data were used to test the spatial accuracy of different LSM methods: (a) ULSM with five different permutationbased thresholding approaches; (b) MLSM with voxel-level lesion data using two different approaches; and (c) MLSM with dimensionreduced lesion data using three different strategies for feature reduction. Obtained results were compared across different anatomical target locations, sample sizes, noise levels, lesion mask smoothing values, types of networks, and false positive simulations. We used a number of different distance-and overlap-based spatial metrics as indices of mapping accuracy. We also compared performance of these LSM methods using real behavioral data (language scores) with multiple demographic and sampling covariates, along with subsampling to check the stability and agreement across methods. Our goal was to provide the first comprehensive comparison of ULSM and MLSM methods, in order to afford guidance on selecting the most appropriate LSM method(s) for a particular study with a specific lesion dataset.

| Participants
For the simulation analyses, lesion masks from 340 chronic left hemisphere stroke patients were obtained from two different sources: our Northern California stroke dataset (n = 209, NorCal) and the Moss Rehabilitation stroke dataset provided with the open-source LESYMAP software (n = 131, LESYMAP, Pustina et al., 2018). Synthetic behavioral scores were based on lesion load to different cortical areas (described further below).
For the analysis of real behavioral data, we analyzed language data and lesion masks from a subset of patients in the NorCal database (n = 168; 36 female) who completed behavioral testing and met the following inclusion criteria: History of a single left hemisphere stroke (including both embolic and hemorrhagic etiologies), premorbidly right-handed (based on the Edinburgh Handedness Inventory), native English speaker (English by age 5), minimum high school or equivalent education (i.e., 12 years), in the chronic stage of recovery (at least 12 months poststroke) at the time of behavioral testing, no other neurologic or severe psychiatric history (e.g., Parkinson's, dementia, schizophrenia), and no substance abuse history. The mean age of this subset of patients was 61.0 years (range 31-86, SD = 11.2), mean education was 14.9 (range 12-22, SD = 2.4), and mean months poststroke was 51.4 (range 12-271, SD = 54.0). All patients were administered the Western Aphasia Battery (WAB, Kertesz, 1982Kertesz, , 2007, which classified 47 patients with anomic aphasia, 45 with Broca's aphasia, 6 with conduction aphasia, 4 with global aphasia, 1 with transcortical motor aphasia, 3 with transcortical sensory aphasia, 14 with Wernicke's aphasia, and 48 patients who scored within normal limits (i.e., overall WAB language score of ≥93.8 points out of 100). This latter group included patients with very mild aphasic symptoms, such as mild word-finding difficulty.

| Behavioral data
Data for the LSM analyses with real behavioral scores were derived from 168 patients in the NorCal stroke dataset. Patients were tested on the WAB (Kertesz, 1982(Kertesz, , 2007, which consists of several subtests measuring a wide range of speech and language functions. Here, we analyzed the most reliable and least-confounded speech-language scores on the WAB, which index three distinct language domains: speech fluency, single-word auditory comprehension, and verbal repetition. All patients signed consent forms and were tested in accordance with the Helsinki Declaration.

| Imaging and lesion reconstructions
Real lesion masks were obtained from two different sources as detailed previously. All lesion masks were reconstructed from MRI or CT data acquired during the chronic phase of stroke (at least 2 months poststroke). Detailed information about data acquisition, lesion reconstruction, and normalization procedures for the NorCal and LESYMAP datasets can be found in Baldo et al. (2013) and Pustina et al. (2018), respectively. The lesion masks were converted to standard MNI space with a 2 mm isovoxel resolution. The overlay of patients' lesions from the two different databases is shown in

| LSM methods
In the current study, we compared five ULSM and eight MLSM methods, using both synthetic and real behavioral data. All ULSM and MLSM methods discussed in this paper were implemented in a downloadable MATLAB script (freely available at https://www.nitrc.org/ projects/clsm/) that was based on the original VLSM software developed by Stephen Wilson (Bates et al., 2003;Wilson et al., 2010).

| ULSM methods
All ULSM variants were conducted using a linear regression with a voxel lesion value as the dependent variable (lesioned or not), the continuous behavioral data as the independent variable, and lesion volume as a covariate (Wilson et al., 2010;Baldo, Ivanova, Herron, Wilson, & Dronkers, in press). Linear regression was chosen because it is flexible, popular, and applicable to multiple data types (categorical or continuous) with multiple covariates. Also, the synthetic behavioral data contained linear effects, similar to other recent simulation studies (e.g., Pustina et al., 2018). Both the synthetic and real behavioral datasets were well controlled with respect to outliers; all the z-scores were within 3 SD of the mean.
To correct for multiple comparisons, we used permutation-based thresholding, which is a nonparametric approach to FWER correction that randomly permutes behavioral scores and which records the relevant t-value or maximal cluster size for each permutation with the final threshold set at p = .05 (Kimberg et al., 2007;Mirman et al., 2018). Unlike some previous methodological LSM papers, we did not evaluate the performance of the ULSM method with FDR-correction, given that it is a fundamentally inappropriate correction for lesion data (Baldo et al., 2012;Mirman et al., 2018). Finally, the use of permutation-based significance testing with all of the ULSM methods protects against inflated false positives that can accompany non-Gaussian noise that may appear in the real data (all synthetic behavioral noise was Gaussian).
ULSM results were generated with five different nonparametric FWER thresholding approaches that are commonly used in contemporary ULSM studies: 1. Maximum statistical t-value (ULSM T-max). This statistic corresponds to the most standard and conservative version of the permutation-based FWER-correction.
3. Cluster-size thresholding with a fixed voxel-wise threshold of

| MLSM methods
MLSM methods included multivariate regression methods with voxellevel lesion data (two approaches) and with dimension-reduced lesion data (three approaches). For all MLSM methods, lesion volume was regressed out of both the behavioral and the lesion variables. While other corrections for lesion volume are possible, such as direct total lesion volume control (Zhang et al., 2014), we opted for the most conservative option here based on recommendations from the only study to date that systematically tested different lesion volume corrections (DeMarco & Turkeltaub, 2018). We normalized both lesion and behavioral data to SD of 1 and centered the behavioral data (to mean 0), as is customary for multivariate methods to optimize regression estimation. Finally, permutation testing identical to that used with ULSM methods above (Kimberg et al., 2007;Mirman et al., 2018), was also applied to the voxel-wise feature weights obtained for the MSLM methods. The resultant values were used to threshold and identify significant voxels in the maps at p < .05, using the maximum voxel-wise feature weight value obtained (as in the ULSM T-max). From here on, we use the term LSM statistical values to refer to both the actual voxel-wise statistics for ULSM methods and voxel-wise feature weights for MLSM methods that are ascertained following permutation thresholding and presented in the resulting LSM map output.

MLSM methods with voxel-level lesion data.
SVR differs from ordinary multivariate regression in a number of ways, as it computes a solution based on many voxels' lesion status given the relatively small number of patients. First, it incorporates two regularization hyperparameters that help the regression model keep the model's parameter values small (to avoid overfitting), while at the same time controlling the model's prediction accuracy by partially ignoring small fitting errors. Second, SVR incorporates a radial basis kernel function that implicitly projects the lesion data into a highdimensional space in order to help model fitting succeed, in part by allowing some nonlinear effects to be incorporated into the model. SVR has been used in several previous LSM studies (Ghaleh et al., 2017;Griffis, Nenert, Allendorfer, & Szaflarski, 2017;Zhang et al., 2014). We used the SVR routine encoded as part of the SVR- The second voxel-level MLSM regression method was partial least squares (PLS), which jointly extracts dual behavior and lesion factors that maximize the variance between behavior and lesion locations in a single step. PLS algorithms (and closely related canonical correlation algorithms) have been developed extensively in bioinformatics for use in genetics where there is a similar "wide" data structure: there are far more genes/voxels to be considered in a regression solution than there are subjects providing such data (Boulesteix & Strimmer, 2006). PLS has also previously been used in LSM (Phan et al., 2010). However, we used a basic version of PLS regression based on the singular value decomposition (SVD) function (Abdi, 2010;Krishnan, Williams, Randal, & Abdi, 2011) because it is known to be both a fast and reliable regression technique over wide data sets containing highly correlated variables, which is important given the number of simulations (using permutation testing) that were run. Although this version of PLS regression is known to produce "dense" solutions (Mehmood, Liland, Snipen, & Saebø, 2012), resulting in (overly) large clusters, another important consideration for including it here as an exemplar of this class of algorithms is that it is easy to generalize basic PLS regression to integrate multiple target behaviors simultaneously (Abdi, 2010); an inviting prospect for investigating behavioral test batteries used to assess patient populations.
MLSM methods with dimension-reduced lesion data.
Three different types of MLSM data reduction methods were tested in the current study. These methods reduce the spatial dimensionality of the lesion data first without considering behavior. Lesion status of thousands of voxels is reduced to a number of spatial lesion components that is fewer than the number of patients in the analysis. This results in a more tractable system of equations to solve, and the components' estimated weights are then transformed back into spatial maps, relating brain areas to the behavior being investigated. In other words, lesion status in all voxels for each patient are replaced with sums of weighted lesion components for that patient. SVD ("svd", MatLab v.7) (Ramsey et al., 2017) identifies an ordered set of orthogonal spatial components, each a weighted mixture of all voxels inside the lesion mask, that explain as much lesion variance across all voxels with as few of the ordered components as possible. These components can be linearly combined to reconstruct each patient's lesion mask. Our preliminary trials found that most individual patient lesion masks were well-reconstructed (median dice >0.98; mean dice >0.85) when 90% of the cumulative variance in the SVD diagonal matrix was accounted for. Thus, in this data-reduction version of MLSM, we used the number of components (approximately equal to half the number of patients) required to explain 90% variance in order to speed up LSM computations.
The second data reduction method we used was a logistic principal component analysis (LPCA) (Schein, Saul, & Ungar, 2003), which iteratively identifies a set of ordered spatial components whose lesion incidence maps are orthogonal under the logistic function (Siegel et al., 2016). Given the better fit of method to lesion data type, trials with lesion masks showed that approximately the first 40 LPCA components were able to reconstruct individual patient lesion masks quite well (median dice >0.99; mean dice >0.9). Thus, the number of components we used for the LPCA data reduction was the number of patients capped at 40 components in order to speed up LSM estimates.
The third data reduction method was an independent component analysis (ICA) (FastICA v2.5 as used by Hyvärinen & Oja, 2000) which is a generalization of PCA. The latter method is widely employed in fMRI studies for clustering brain regions, and although it has not been typically used in lesion analysis studies, we included it here for exploratory purposes.
ICA in this context estimates independent linear mixtures of lesion incidence voxel data that are the most non-Gaussian sources found within the data. We used the default cubic function as the fixed-point nonlinearity for finding components under coarse iterations first and then used a hyperbolic tangent function (default setting) for fine iterations in order to reflect the bounded nature of lesion data. For ICA, we used the same number of components as with LPCA (maximum of 40), in order to see if ICA can outperform LPCA given its usefulness in other spatial dimension reduction applications in neuroimaging (Calhoun, Liu, & Adali, 2009).
After application of the three data reduction approaches described above, an elastic net linear regression ("glmnet" package; Qian, Hastie, Friedman, Tibshirani, & Simon, 2013;Tibshirani et al., 2010) was performed with the target behavior (real or simulated) as the dependent variable along with the data reduced spatial lesion components and the lesion size covariate. While data reduction is not required for MLSM, it is a necessary step for implementation of elastic net regression. Elastic net regression is only appropriate for relatively low-dimensional datasets, because for high-dimensional data (original voxel-level lesion data), it will select too many voxels unrelated to behavior (Gilhodes et al., 2020).
We used two different elastic net regressions to see if either is superior in producing accurate or reliable LSM maps: one near to a pure lasso (which we will call "L1") case (95% L1 penalty mixed with a 5% L2 penalty) and one with the opposite mixture (95% L2 and 5% L1), a ridge ("L2") case. The elastic net regressions use cross-validation to solve for the penalty hyperparameter that best fits the data. Thus, in total six variants of MLSM methods with dimension-reduced lesion data were tested: SVD-L1, SVD-L2, LPCA-L1, LPCA-L2, ICA-L1, and ICA-L2. For technical details on implementation of the MLSM methods described above see Appendix B.

| Simulations with synthetic behavioral data
Three sets of simulations were performed with all ULSM and MLSM methods using synthetic behavioral scores and real lesion masks: single anatomic target, dual anatomic targets, and zero anatomic target (i.e., false positive simulation). In each simulation, we varied several factors (described below for each simulation) in a fully crossed manner in order to systematically compare effect sizes and significance across the different ULSM and MLSM methods. For each simulation analysis, the specified number of lesion masks was randomly selected from one of the two datasets (i.e., not mixing NorCal and LESYMAP masks together). For all analyses, we only included voxels in which at least five patients had lesions, and which had statistical power ≥.1 at p < .01 (Hsieh, Bloch, & Larsen, 1998). We also performed behavioral value outlier scrubbing at jzj ≥ 3.0, which was particularly important at higher behavioral noise levels.
Synthetic behavioral scores (also called "artificial" or "fake" in the literature) for the single and dual anatomical target simulations were derived from the lesion load to the target anatomical parcels (or ROIs).
For simple single target simulations, the synthetic behavioral score was calculated as the fraction of the target anatomical parcel that was spared (see next section for more details), that is the synthetic score was directly proportional to the lesion load of that anatomical parcel.
Use of synthetic behavioral scores allows one to determine how well the different LSM methods are able to localize behavior, since we know the ground truth as to exactly which region in the brain it should localize to (i.e., the target anatomical parcel) . To create target anatomical parcels, we used ROI masks of gray matter areas in the left middle cerebral artery region from FSL's version of the Harvard-Oxford (H-O) atlas, thresholded at 50% incidence. We used 16 such parcels that had 5% or greater lesioned area within at least 25% of the lesion masks. To create a set of smaller parcels, each of these 16 parcels was divided into two sections along the axis of maximal spatial extent. Two of the subdivided parcels failed to intersect the lesion mask sufficiently according to the above criteria, rendering a total of 30 smaller parcels (see Figure 2). We specifically chose larger ROIs (similar to Mirman et al., 2018), because small ROIs are unlikely to be accurately identified in patients who generally have much larger lesions than focal fMRI activation areas. Further, even assuming that fMRI properly delineates the size and location of specific functional areas, in the chronic stage of stroke recovery, these functional areas are likely to be altered by neural reorganization (Kiran, Meier, & Johnson, 2019;Stefaniak, Halai, & Ralph, 2019).

| Single anatomical target simulations
In the single anatomical target simulations, the synthetic behavioral data for each patient were calculated as one minus the lesion load (fraction of the target anatomical parcel covered by the patient's lesion), before noise was added. Accordingly, a score of 1 indicated "perfect performance"-complete sparing of the target parcel by the lesion, and a score of 0 indicated "complete impairment"-full coverage of the target parcel by the lesion. For the simulations, we varied the following factors in a fully crossed manner (obtaining all possible combinations of these factors): • Number of patients (n = 32, 48, 64, 80, 96, 112, 128). Additional simulations were run with n = 144-208 patients from the NorCal dataset for descriptive purposes only. Patients were randomly selected from one dataset on each run. • Size of parcels (16 larger or 30 smaller as described above).
• Anatomical target parcel (see Figure 2 for list).

| Dual anatomical target simulations
In the dual anatomical target simulations, two spatially distinct target parcels were used to simulate a minimal brain "network." We considered three different types of dual-target networks: F I G U R E 2 Target anatomical parcels (n = 30) used to generate synthetic behavioral scores • Redundant network in which the behavioral score is reduced only if both target parcels are lesioned (one minus the minimum lesion load of the two parcels is used to generate the synthetic behavioral score).
• Fragile network in which the behavioral score is reduced if either target is lesioned (one minus the maximal lesion load of the two parcels is used to generate the synthetic behavioral score; corresponds to what Pustina et al., 2018 called the "OR" rule for generating multiregion simulations and approximates the partial injury problem described in Rorden, Fridriksson, & Karnath, 2009).
• Extended network in which the behavioral score is reduced proportionally to the overall damage to the two regions, which is similar to the single target simulation except the parcel is divided into two spatially separate components (one minus the average lesion load of the two parcels is used to generate the synthetic behavioral score; corresponds to what Pustina et al., 2018 called the "AND" rule).
For the dual-target simulations, we only analyzed results with the larger 16 parcels, moderate behavioral noise level (0.36 SD), and lesion smoothing at 4 mm FWHM. We tested all 120 pairwise combinations of the target parcels for each type of network. The number of patients was varied systematically from n = 64 to 208. We did not use n = 32 or n = 48, because preliminary results showed a lack of power with this sample size. We only used the NorCal dataset for this analysis for consistency across these simulations.

| Measures of LSM success in simulation analyses
As a proxy for statistical power in both single and dual-target simulations, we calculated the percentage of trials that yielded any significant (above-threshold) LSM statistical values (referred to here as the probability of obtaining a positive result).
To gauge the accuracy of an LSM method in our single and dual anatomical target simulations, we calculated two types of accuracy measures: distance and overlap. For each measure, the target anatomical parcel, whose lesion load was used to generate the specific synthetic behavioral score, was compared to the LSM output map (LSM Overlap-based measures included overlap and weighted overlap metrics between target anatomical parcel(s) and LSM output map.
First, we used the dice coefficient, which is a simple measure of overlap between two binary maps, with 1 being perfect overlap between the two maps and with zero voxels falling outside of the overlapping area, and 0 being no overlap between the two (calculated as proportion of area 2B to the combined areas of the anatomical target and LSM cluster [A + 2B + C] in Figure 4). In addition, we also looked at the proportion of false negatives (part of the anatomical target not covered by statistically significant LSM values-A) and false positives (statistically significant LSM values falling outside the anatomical target-C) relative to the combined areas of the anatomical target and LSM cluster (A + 2B + C) to better understand what is driving specific dice values (see Figure 4). While the dice coefficient remains a standard means of measuring mapping accuracy, there are two clear limitations to this measure. First, the dice coefficient is greatly dependent on the size of the parcels being compared, with larger parcels generating a larger dice index compared to smaller regions, even when the relative overlap is smaller (Pustina et al., 2016). Second, the actual statistical values of the LSM output are not taken into account in the calculation of the dice index. To account for this latter limitation, we also computed a one-sided Kuiper (OSK) distribution difference (Rubin, 1969) between statistically significant LSM values inside the target anatomical parcel(s) versus outside. This measure compares the LSM statistical values outside the target (C) to those inside the target (B) and it also assigns zero values to target areas not reaching threshold in the LSM output (A) (see Figure 4). The idea is that we want to reward finding LSM hotspots inside the anatomical target(s) while ignoring lower values in LSM clusters that are outside the target. OSK values range from −1 to 1 with 1 representing that the entire target is covered with LSM statistical values higher than all those outside the target; −1 representing that the lowest LSM values (or none at all) are inside the target, and 0 representing no difference in LSM statistic distributions inside versus outside the target anatomical parcel.

| LSM analysis with real behavioral scores
In addition to the simulation analysis with synthetic behavioral scores, we also compared ULSM and MLSM methods using real behavioral scores. These data were generated from a subset of patients in the NorCal dataset (n = 168) who met specific inclusion criteria (described above in the participants section) and who were tested on speech fluency, single-word auditory comprehension, and repetition subtests from the WAB (Kertesz, 1982(Kertesz, , 2007. Since the aim of the current paper was to evaluate performance of different LSM methods with real behavioral data under typical conditions, we used a number of standard covariates. So in these analyses, in addition to lesion volume, we also covaried for age, education, gender, months poststroke (log-transformed), test site (referring to one of the two locations where behavioral data were collected), image resolution (referring to the resolution of the original scan: low or high), and overall aphasia severity (WAB AQ minus the target subscore). The last covariate allowed us to account for overall aphasia severity, while focusing on the specific language function. Lesion smoothing was done at 4 mm FWHM. The minimum power per voxel was 0.25 at p < .001 for a d 0 of 1.
We also performed 100 repetitions of the same 3 analyses but with 75% subsampling (n = 126) of the original full cohort in order to determine the stability of the solutions for each given LSM method.
We considered the original LSM map generated with the full cohort (n = 168) to be the "target parcel" in this case and correspondingly analyzed the stability of the results obtained on subsequent runs (n = 126) with the same metrics we used to analyze the single target simulation results. We then tested whether there were significant differences between the different LSM methods with respect to the stability of the results produced.

| Statistical analyses
Simulations with synthetic behavioral data were performed in a fully Analysis of important differences between each simulation's factors was established using mixed between/within analysis of variance (ANOVA) applied to each fully crossed dataset type where we treated anatomical target parcel location as the random factor.
We used high-speed ANOVA software (CLEAVE, nitrc.og.projects/ cleave) to compute partial omega squared) effect sizes (Olejnik & Algina, 2004) for factors and first-order interactions. Typically, weak, moderate, and strong partial omega effect sizes are taken to be 0.2, 0.5, and 0.8, respectively (Keppel, 1991). We used partial omega squared cutoff of 0.35 to restrict ourselves to reporting moderate or strong effects (unless otherwise noted). Finally, we note that the distance, dice, and OSK values often did not conform to a Gaussian distribution as required by ANOVA. In these instances, we used a Kumaraswamy distribution (Jones, 2009)  Additionally, we also looked at the different distance-based spatial accuracy metrics as factors and explored whether they yielded similar accuracy estimates or not. In the results, we concentrate specifically on effect sizes rather than significance testing, given that large-scale simulations produce so much data, that even tiny differences can be significant (Kirk, 2007;Schmidt & Hunter, 1995;Stang, Poole, & Kuss, 2010).
The same ANOVA analyses were performed on the subsample analysis of the real behavioral data by using the full analysis LSM maps as targets for their respective LSM variant (e.g., SVR target for SVR LSM analysis). The only addition was that we included the maximum statistic location of the target ROI as another centroid to anchor three more distances per LSM map, because these target ROIs (the n = 168 maps) already have statistics.  Table 1

| Spatial accuracy using distance-based metrics
The spatial accuracy of the different LSM methods using distance-based metrics was evaluated (see Table 1 for mean displacement error values). There was a main effect of LSM Method on accuracy (ω 2 p = .84). ULSM methods demonstrated numerically higher mean accuracy across all output map locations and target locations, with the average displacement error ranging from 6 to 7.5 mm.
Average displacement error across MLSM methods (excepting SVR) ranged from 8 to 13 mm. Across sample sizes, the conservative ULSM methods (T-max and T-nu = 125) and SVR produced the most accurate maps. Figure 5 shows the spatial displacement error based on mean Max LSM to Closest target (these two metrics in combination provide the highest accuracy estimates; more on this below) of the different LSM methods at different sample sizes. As can be seen in  Tables with the mean simulation values for these analyses and further information on different distance-based metrics can be found in Appendix A (Tables A1-A3).
An additional post hoc analysis was run with smaller-sized anatomical targets (n = 30) and the highest behavioral noise level (0.71 SD), with and without lesion size as a covariate. Using lesion size as a covariate resulted in higher spatial accuracy of LSM output maps across all sample sizes and all LSM methods (ω 2 p = .87, mean displacement of 8.6 mm vs. 12.1 mm without it). The effect of covarying for lesion size varied as factor of LSM Method (ω 2 p = .43) with MLSM methods showing more improvement in spatial accuracy (3-8 mm) relative to ULSM methods (2-3 mm) with the lesion size covariate (see Appendix A, Table A4). In contrast to the dice coefficient values, OSK distribution values were substantially higher and showed higher variability across factor levels, rendering it more useful than the dice coefficient for evaluation of spatial accuracy (see Table 1). The OSK distribution statistic greatly

| Dual anatomical target simulations
Dual-target simulations were run on three network types: redundant, fragile, and extended. LSM Method, sample size and type of network all had a strong effect on the probability of obtaining a positive result F I G U R E 8 Dice index as a function of lesion symptom mapping (LSM) method and sample size (top panel) or behavioral noise level (bottom panel) for single target simulations (ω 2 p ranged from .80 to .86; see Table 2 for mean simulation values). A sample size of 64 or larger was required to achieve a probability of 80% for all three network types. Two-way interactions between these factors demonstrated medium to large effect size (ω 2 p ranged from .36 to .50) (see Appendix A, Tables A5 and A6 for relevant data). The Redundant network was the hardest to detect for all LSM methods, requiring a much larger sample size to obtain probability levels similar to that of the fragile and extended networks.  Table A7). Patterns in the proportion of false negatives and false negatives across sample sizes were similar to those observed in single target simulations (see Appendix A, Table A8). See also Figure 7 for LSM output maps for the three network types for planum temporale and pars triangularis regions across different LSM methods.

| False positive simulations
For the false positive simulations, SVR and ULSM T-max performed best, yielding the smallest number of clusters with the smallest number of voxels in those 5% of trials when a false positive solution was produced (see Table 3

| LSM analysis of real behavioral data
We next compared results across the different LSM methods, using real behavioral data collected from 168 patients from the NorCal dataset. The behavioral data consisted of WAB language subscores F I G U R E 9 Proportion of false negatives (part of the anatomical target not covered by statistically significant lesion symptom mapping [LSM] values) shown in the top panel and false positives (statistically significant LSM values falling outside the anatomical target) shown in the bottom panel. Proportions are calculated relative to the combined areas of the anatomical target and LSM cluster as a function of LSM method and sample size for single target simulations for speech fluency, single-word auditory comprehension, and verbal repetition (Kertesz, 1982(Kertesz, , 2007. Results showed distinct areas identified for different behavioral scores, such as frontal regions for speech fluency and posterior temporal regions for single word comprehension (see Figure 13 for LSM maps across methods). However, as can be seen in Figure 13, substantial variability in the identified regions was observed between different LSM methods. MLSM methods tended to identify more regions, including regions not typically associated with the behavior (e.g., frontal regions for single-word comprehension). In contrast, ULSM maps typically clustered around the maximum statistic location more traditionally associated with the behavior (e.g., posterior temporal cortex for single-word comprehension). The methods that provided more sparse solutions (smaller LSM clusters) showed better differentiation between different behaviors under examination (e.g., conservative ULSM methods and ICA).
Since there is no ground truth in the case of real behavioral data, we evaluated accuracy in terms of the stability of the obtained solutions using a random subsampling approach for n = 126 (using the full sample as the target map). The LSM maps produced in the subsample analysis varied between themselves and often clearly differed from the full sample map (see again Figure 13). Stability varied substantially across LSM techniques, as shown by spatial metrics averaged across the 100 subsample analyses and by the average correlation between the full LSM output and each of the subsample analyses (see Table 4). The methods that provided more dense solutions (larger LSM clusters) tended to show higher stability of results (e.g., liberal ULSM methods and PLS). Moderately conservative ULSM methods showed higher stability of results compared to MLSM as indicated by lower average displacement and higher OSK statistic, as well as higher correlations.
Finally, when evaluating the stability of the subsampling analyses substantial variability across different distance-based metrics was observed (see Table 5). The maximum LSM statistic of each subset LSM map usually landed on or very near the full target LSM map (0.3 mm), but often not very close to the full map's maximum value (20.6 mm). Since the location of the maximum statistic is most commonly reported in LSM papers, we also compared the stability of the  Table 4). Again, moderately conservative ULSM methods showed higher stability (range: 0.6-0.8 mm) relative to MLSM methods (range: 0.8-2.6 mm) with the exception of PLS.

| DISCUSSION
In the current study, we conducted the first comprehensive, empirical comparison of several univariate and multivariate LSM methods. The methods were evaluated with respect to several characteristics, including the probability of obtaining positive results, spatial accuracy, and robustness. Both synthetic and real datasets were employed across a range of relevant parameters, including sample size, lesion smoothing, noise level, and type of network. A number of different metrics were used to evaluate spatial accuracy. Cumulatively, our analyses indicated that both ULSM and MLSM methods can be equally robust in locating brain behavior relationships, depending on the design of the study, the research question being asked, and with proper spatial metrics.
The results provide crucial insights into the accuracy of different LSM methods and their susceptibility to artifact, providing a first of its kind data-driven navigational guide for users of LSM analyses. Based on these findings, we suggest that both ULSM and MLSM methods should be used in tandem, given an adequate sample size for each method, to increase confidence in the results of an LSM analysis.

| Comparison of ULSM and MLSM methods
Real findings will be consistent across different classes of LSM methods, but spurious findings will be inconsistent.
With respect to LSM performance on the real behavioral data, moderately conservative ULSM (T-nu = 125 and T-0.001) and specific MLSM methods (ICA, SVD) identified largely similar regions. This is again in line with previous empirical studies that employed both ULSM and MLSM methods to investigate brain-behavior relationships and yielded highly coherent results between the two approaches (e.g., Fridriksson et al., 2018;Thye & Mirman, 2018). However, some MLSM methods also identified brain regions that are not typically associated with the behavior under examination. Also, stability (robustness) of the obtained solutions was generally higher for ULSM, relative to MLSM methods. On the whole, methods that provided more dense solutions tended to show higher stability of results (e.g., PLS, ULSM T A B L E 3 Mean values of above-threshold (significant) cluster characteristics in false positive simulations for the 5% of trials for each LSM method that yielded false positive results   T-0.01). However, these methods showed less spatial discrimination with results obtained for different language functions yielding more similar LSM output maps relative to the sparse solutions. The opposite pattern was observed for methods providing sparse solutions are based on statistical, but not causal inference, they do not demonstrate a notable advantage over ULSM methods in overcoming these limitations (Sperber, 2020). Moreover, in the current study, data reduction methods in particular often showed pronounced spatial distortion. It is very likely that these methods can enhance vasculaturebased effects as the brain's vasculature greatly drives lesion structure similarities and thus determines the dimension-reduced spatial components that these methods use. While no single method was able to fully overcome lesion-anatomical biases and produce a thresholded LSM map that even approached an exclusive delineation of brain regions associated with the target behavior, under specific conditions certain methods performed better than others. ULSM methods demonstrated particularly good performance in terms of accuracy and robustness for identifying single targets. Here, MLSM methods clearly lagged in power, accuracy, susceptibility to behavioral noise, and stability. Although, not all MLSM methods were equally good at identifying networks, in general, they performed better than ULSM in uncovering various functional dependencies in the lesion data (for a related argument, see Sperber, 2020). This underscores the notion that no single LSM method can provide an ultimate solution for establishing brain correlates of cognitive functioning.

| Impact of other parameters on LSM accuracy
In addition to the comparison of different LSM methods, we also investigated the impact of a number of other parameters on mapping accuracy across methods, including anatomical target location and size, mask smoothing, behavioral noise level, sample size, and lesion volume. With respect to target location, we found variable accuracy across spatial locations, with especially poor performance in cortical locations on the edge of the lesion masks (areas of lower power). The size of target parcels (small or large) did not have a substantial impact on the LSM results for distance-based metrics, with accuracy being only slightly higher for the smaller parcels in the single anatomical target simulations. This implies that the current findings with distancebased metrics are applicable to LSM analyses irrespective of the size of the hypothetical target parcel(s). The size of the target parcels did have an impact on dice scores, with larger parcels showing higher dice coefficient values. These findings reflect a previous interpretation that the dice index is highly biased by the size of the chosen target parcels  and is not the optimal measure for interpretation of LSM simulation results.
We also compared the effects of lesion mask smoothing at 0 (no smoothing) versus 4 mm Gaussian FWHM. The 4 mm smoothing led to improved accuracy of localization across all LSM methods. We suspect that with real data, smoothing should be even more beneficial given that in simulations with synthetic behavioral scores there was effectively no lesion delineation error, anatomical normalization error, or anatomical size/shape mismatch across patients as there ordinarily is (typically 5 mm in extent) in stroke patients (Crinion et al., 2007;Rorden, Bonilha, Fridriksson, Bender, & Karnath, 2012).
Behavioral noise level was also manipulated in the current study, in order to investigate its impact on accuracy of localization. More behavioral noise in small sample sizes (n < 100), led to a reduction in the probability of obtaining positive results, particularly for MLSM methods. At larger sample sizes (n > 100), behavioral noise did not considerably diminish this probability. Behavioral noise also negatively impacted the accuracy of mapping, as evidenced by reduced OSK values, but had no noticeable impact on distance-based measures of spatial accuracy.
With respect to sample size, the current findings showed increased probability of obtaining significant results and accuracy with increasing numbers of patients, not surprisingly. It is important to note, however, that the benefit of increased sample sizes on spatial accuracy plateaued at about 130 participants for all LSM methods in single target simulations. LSM methods do not appear to be spatially unbiased statistical methods, and stroke lesions are not entirely random (Sperber & Karnath, 2016), so beyond a certain sample size the increase in coverage and in variability of lesion patterns is minimal leading to muted effects on accuracy.
Finally, including lesion volume as a covariate significantly improved spatial accuracy across all LSM methods. The importance of including lesion volume as a covariate has been shown in previous studies for ULSM (Sperber & Karnath, 2017) and some MLSM methods (DeMarco & Turkeltaub, 2018;Zhang et al., 2014). Our current findings further reinforce the stipulation that a lesion volume covariate should be included in all types of LSM analyses.

| Different metrics of LSM accuracy
Once sufficient power is achieved in LSM analysis, the greatest chal-  The final and very critical recommendation arising from our current study concerns which metrics should be reported in LSM studies.
As we have empirically shown here, LSM results are always spatially distorted to a certain degree-no LSM method provided output maps that were perfectly spatially accurate-so it is important to report localization using the most robust metrics that take lyses. In the current study, we used linear regression rather than logistic regression (appropriate when the dependent variable is bounded) because linear regression is much faster to compute, is more familiar to users, and is standard in previous studies (Akinina et al., 2019;Baldo et al., 2013Baldo et al., , 2018Ivanova et al., 2018;Mirman, Chen, et al., 2015). While linear regression is sufficient to model the perfectly linear relationships in synthetic data, real behavioral data might require a different approach such as the use of nonparametric measures (Rorden et al., 2007).
Another significant limitation of the current study is the lack of a sparse version of PLS (Phan et al., 2010) or similar methods, like SCCAN . We did not implement SCCAN in the current study because it is highly computationally intensive as one needs to determine several hyperparameters in order to optimize performance for a given hypothesis/target configuration.

| CONCLUSIONS
We hope that the current systematic approach to evaluating the prob-  ULSM methods demonstrated the greatest probability to detect significant voxels. MLSM methods required more participants to achieve comparable probability levels ( Figure A1, top panel). Sample size (ω 2 p = .79) and behavioral noise level (ω 2 p = .74) strongly impacted the probability of obtaining positive results across all LSM methods: larger sample size and less behavioral noise led to a higher probability of detecting significant voxels. Also, first-order interactions between these factors demonstrated a large effect size (ω 2 p ranged from .61 to .73). At higher behavioral noise levels, a larger sample size was required to achieve similar probability levels ( Figure A1, middle panel).
Probability of obtaining positive results was disproportionately reduced for MLSM methods relative to ULSM at higher levels of behavioral noise, with the exception of SVR ( Figure A1, bottom panel). data into a smaller (number of subjects × umber of patterns) that records how well each subject's lesion matches each spatial pattern component common to the group of subjects extracted by the data reduction procedure.

Impact of different distance-based metrics
a. An SVD is performed on the matrix of lesion data (number of subjects × number of voxels). A number of least significant (according to the SVD eigenvalues) components/patterns are eliminated according to a specified formula based on the sample size and variance remaining. This additional data reduction speeds up processing, eliminating spatial patterns where few subjects have lesions overlapping.
b. An LPCA (Schein et al., 2003) is an iterative method that similarly reduces bounded data into components/patterns, and again only the LPCA components/patterns capturing the most variance are kept for the next step.
c. ICA (Hyvärinen & Oja, 2000) is computed is similarly computed and the highest variance spatial patterns, and subject lesion overlap with each pattern are kept for the next step.
4. Elastic net regression: An underdetermined linear regression (the number of subjects is generally only two to three times the number of spatial patterns being produced by the data reduction step) is implemented using the elastic net method.
a. First, the subject overlap coefficient matrix (number of subjects × number of patterns) is normalized to N (0,1) in each subject. Similarly, the target behavioral vector is normalized to N (0,1).
b. An elastic net regression (Tibshirani et al., 2010) is performed using 10-fold cross validation to compute the optimal coefficient for the penalty term to control data noise from overly influencing the solution. The elastic net type can be chosen to be either a near-ridge type (L2) or a near LASSO type (L1) that can affect the pattern-sparsity of the results. The solution obtained is a vector of values representing how well each data reduction-produced voxel pattern correlates with the target behavioral variable. The cross-validation also produces an estimate of the variance of the target behavior explained by the subject voxel/pattern data.
c. The regression solution is then converted, using the saved spatial voxel patterns computed in the data reduction step, into (arbitrary scale) feature weights for all voxels.
5. Permutation test thresholding: The procedure above produces a feature weight at each voxel that was included in the analysis after preprocessing Step 1.d. Further thresholding of these data reduction MLSM maps is obtained by performing permutation testing as indicted in the main text by randomly permuting the behavioral data (but not lesion size) and repeating Steps 2-4 and computing maps to produce a maximum feature weight across each such permuted data map that can be used to threshold the true feature weight map at a familywise p < .05.

MLSM methods with voxel-level lesion data
The following describes the specific computational steps used in PLS method as implemented in the open software at https://www.nitrc.
org/projects/clsm used in the current study.
1. Preprocessing: Same as in the Data Reduction method section above.
2. Covariate elimination: Same as in the Data Reduction method section above.
3. PLS regression: An SVD-based PLS regression (Abdi, 2010) is used to jointly produce a voxel lesion solution pattern that explains (linearly) the highest amount of target behavioral variance. This is distinct from the data reduction methods above that produce the most common spatial lesion pattern of the subject population. PLS regression was designed to handle highly correlated variables such as voxel lesion data; however, it often produces dense solutions.
The amount of variance explained by PLS is computed using predicted residual errors within a 10-fold cross validation.
4. Permutation test thresholding: Similar to the Data Reduction method permutation section above, except redoing PLS Steps 2 and 3 above on permuted behavioral data to compute the maximal PLS voxel solution value distribution that provides the familywise p < .05 thresholding.