Lesion to Symptom Mapping

Lesymap uses univariate and multivariate methods to map functional regions of the brain that, when lesioned, cause specific cognitive deficits. It requires is a set of binary lesion maps (nifti files) in template space and the vector of corresponding behavioral scores. Note, lesions must be already registered in template space, use the built-in function registerLesionToTemplate or other ANTs tools to register lesions. Lesymap will check that lesion maps are in the same space before running. For traditional mass-univariate analyses (i.e., BMfast, ttest, etc.), voxels with identical lesion patterns are grouped together in unique patches. Patch-based analysis decreases the number of multiple comparisons and speeds up the analyses. Multivariate analysis are performed using an optimized version of sparse canonical correlations (SCCAN, method='sccan')or support vector regression (method='svr').

lesymap(lesions.list, behavior, mask = NA, patchinfo = NA,
  method = "sccan", correctByLesSize = "none",
  multipleComparison = "fdr", pThreshold = 0.05, flipSign = F,
  minSubjectPerVoxel = "10%", nperm = 1000, saveDir = NA,
  binaryCheck = TRUE, noPatch = FALSE, showInfo = TRUE, ...)

Arguments

lesions.list	list of antsImages, or a vector of filenames, or a single antsImage with 4 dimensions.
behavior	vector of behavioral scores or filename pointing to a file with a single column of numbers.
mask	(default=NA) binary image to select the area where analysis will be performed. If not provided will be computed automatically by thresholding the average lesion map at minSubjectPerVoxel.
patchinfo	(default=NA) an object obtained with the getUniqueLesionPatches function. Useful to set if your run repetitive analysese and want toavoid the computation of patches each time. You can also pass the patchinfo object obtained from a previous analysis.
method	what analysis method to use to run, one of 'BM', 'BMfast', 'ttest', 'welch', 'regres', 'regresfast', 'regresPerm', 'sccan' (default) or 'svr',. BM - Brunner-Munzel non parametric test, also called the Generalized Wilcoxon Test. The BM test is the same test used in the npm/Mricron software (see Rorden (2007)). This method is slow, use 'BMfast" for a compiled faster analysis. BMfast - ultrafast Brunner-Munzel with compiled code. BMfast can be combined with multipleComparison='FWERperm' to perform permutation based thresholding in a short time. ttest - Regular single tailed t-test. Variances of groups are assumed to be equal at each voxel, which may not be true. This is the test used in the voxbo software. Relies on t.test function in R. It is assumed that 0 voxels are healthy, i.e., higher behavioral scores. See the alternative parameter for inverted cases. (see Bates (2003)). welch - t-test that does not assume equal variance between groups. Relies on t.test function in R. regres - linear model between voxel values and behavior. Uses the lm function in R. This is equivalent to a t-test, but is useful when voxel values are continuous. This method is R-based and slow, for faster analysis and to add covariates use the "regresfast" method compiled in LESYMAP. regresfast - fast linear regressions with compiled code. This method allows setting covariates. If covariates are specified the effect of each voxel will be estimated with the formula: behavior ~ voxel + covar1 + covar2 + ... The effect of covariates at each voxel is established with the Freedman-Lane method (see Winkler (2014)). This LESYMAP method allows multiple comparison correction with permutation based methods "FWERperm" and "clusterPerm". regresPerm - linear model between voxel values and behavior. The p-value of each individual voxel is established by permuting voxel values. The lmPerm package is used for this purpose. Note, these permutations do not correct for multiple comparisons, they only establish voxel-wise p-values, a correction for multiple comparisons is still required. chisq - chi-square test between voxel values and behavior. The method is used when behavioral scores are binary (i.e. presence of absence of deficit). Relies on the chisq.test R function. By default this method corrects individual voxel p-values with the Yates method (the same approach offered in the Voxbo software). chisqPerm - chi-square tests. P-values are established through permutation tests instead of regular statistics. Relies on the chisq.test R function. sccan - sparse canonical correlations. Multivariate method that considers all voxels at once by searching for voxel weights that collectively explain behavioral variance. For our purposes, this method can be considered a sparse regression technice. By default, lesymap will run a lengthy procedure to determine the optimal sparseness value (how extensive the results should be). You can set optimizeSparseness=FALSE if you want to skip this optimization. The search for optimal sparsness provides a cross-validated correlation measure that shows how well the sparseness value can predict new patients. If this predictive correlation is below significance (i.e., below pThreshold), the entire solution will be dropped and LESYMAP will return an empty statistical map. If the predictive correlation is significant, LESYMAP will return a statistical image with normalized weights between -1 and 1. The raw SCCAN weights image is returned in rawWeights.img as well. LESYMAP scales and centers both lesion and behavior data before running SCCAN (hardcoded in lsm_sccan). See more details in Pustina (2018) svr - support vector regression. Multivariate method that considers all voxels at once by searching for voxel weights. To establish p-values, a number of permutations are needed as set with SVR.nperm. The current implementation of SVR in LESYMAP is not parallelized and takes many hours to finish all permutations. The SVR method is initially described in Zhang (2014), LESYMAP uses a contribution by the Tuebingen group.
correctByLesSize	whether to correct for lesion size in the analysis. Options are "none", "voxel", "behavior", "both": "none": (default) no correction "voxel": divide voxel values by 1/sqrt(lesionsize). This is the method used in Mirman (2015) and Zhang (2014). This correction works only with 'regres' methods. Two sample comparisons (t-tests and Brunner-Munzel) use binary voxels and will ignore this correction. "behavior": residualize behavioral scores by removing the effect of lesion size. This works on all methods, but is more agressive on results. "both": both voxel and behavior residualized.
multipleComparison	(default='fdr') method to adjust p-values. Standard methods include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr". (see p.adjust ) Permutation methods include: "FWERperm" (permutation based family-wise threshold) is enabled with methods 'BMfast' and 'regresfast'. In this case, many analysis are run with permuted behavioral scores, and the peak score is recorded each time (see Winkler 2014). The optimal threshold is established at 95th percentile of this distribution (or whatever pThreshold you choose). You can choose to use as reference another voxel lower in the ranks by specifying another `v` value (i.e., lesymap(..., v=10) will record the 10th highest voxel). "clusterPerm" (permutation based cluster correction) is enabled for 'regresfast'. It records the maximal cluster size from many random permutations of the behavior score and sets a cluster threshold based on that distribution. You must select pThreshold (voxel-wise, default=0.05) and clusterPermThreshold (cluster-wise, default 0.05) to achieve optimal thresholding with this approach.
pThreshold	(default=0.05) threshold statistics at this p-value (after corrections or permutations)
flipSign	logical (default=FALSE), invert the sign in the statistics image.
minSubjectPerVoxel	(default='10%') remove voxels/patches with lesions in less than X subjects. Value can be speficifed as percentage ('10%') or exact number of subjects (10).
nperm	(default=1000) number of permutations to run when necessary. This is used mostly for univariate analyses, while multivariate methods have their own permutation arguments. Check the documentation of each method to know more.
saveDir	(default=NA) save results in the specified folder.
binaryCheck	logical (default=FALSE), make sure the lesion matrix is 0/1. This will help if lesion maps are drawn in MRIcron or other software which label lesioned voxel with value 255.
noPatch	logical (default=FALSE), if True avoids using patch information and will analyze all voxels. It will take longer and results will be worse due to more multiple comparison corrections. This argument is ignored when performing multivariate analyses, SCCAN or SVR, for which all voxels are always used.
showInfo	logical (default=TRUE), display time-stamped info messages
...	arguments that will be passed down to other functions (i.e., sparsness=0.045)

Value

The following objects are typically found in the returned list:

• stat.img - statistical map

• rawWeights.img - (optional) raw SCCAN weights

• pval.img - (optional) p-values map

• zmap.img - (optional) zscore map

• mask.img - mask used for the analyses

• average.img - map of all lesions averaged, produced only if no mask is defined.

• callinfo - list of details of how you called lesymap

• outputLog - terminal output in a character variable

• perm.vector - (optional) the values obtained from each permutation

• perm.clusterThreshold - (optional) threshold computed for cluster thresholding

• perm.FWERthresh - (optional) threshold computed for FWERperm thresholding

• patchinfo - list of variables describing patch information:

  • patchimg - antsImage with the patch number each voxels belongs to

  • patchimg.samples - antsImage mask with a single voxel per patch

  • patchimg.size - antsImage with the patch size at each voxel

  • patchimg.mask - the mask within which the function will look for patches

  • npatches - number of unique patches in the image

  • nvoxels - total number of lesioned voxels in mask

  • patchvoxels - vector of voxel count for each patch

  • patchvolumes - vector of volume size for each patch

  • patchmatrix - the lesional matrix, ready for use in analyses. Matrix has size NxP (N=number of subjects, P=number of patches)

Details

Several other parameters can be specified to lesymap() which will be passed to other called fuctions. Here are some examples:

permuteNthreshold - (default=9) for Brunner-Munzel tests in method='BMfast' or method='BM'. Voxels lesioned in less than this number of subjects will undergo permutation-based p-value estimation. Useful because the BM test is not valid when comparing groups with N < 9. Issue described in: Medina (2010)

clusterPermThreshold - threshold used to find the optimal cluster size when using multipleComparison='clusterPerm'.

alternative - (default='greater') for two sample tests (ttests and BM). By default LESYMAP computes single tailed p-values assuming that non-lesioned 0 voxels have higher behavioral scores. You can specify the opposite relationship with alternative='less' or compute two tailed p-values with alternative='two.sided'.

covariates - (default=NA) enabled for method = 'regresfast'. This will allow to model the effect of each voxel in the context of other covariates, i.e., formula "behavior ~ voxel + covar1 + covar2 + ...".
I,.e., lesymap(lesions,behavior, method='regresfast', covariates=cbind(lesionsize, age)).
If you choose permutation based thresholding with covariates, lesymap will use the Freedman-Lane method for extracting the unique effect of each voxel (see Winkler 2014, Freedman 1983)

template - antsImage or filename used for plotting the results if a saving directory is specified (see saveDir)

v - (default=1) which voxel to record for permutation based thresholding. Normally the peak voxel is used (1), but other voxels can be recorded. See Mirman (2017) for this approach.

Examples

lesydata = file.path(find.package('LESYMAP'),'extdata')
filenames = Sys.glob(file.path(lesydata, 'lesions', 'Subject*.nii.gz'))
behavior = Sys.glob(file.path(lesydata, 'behavior', 'behavior.txt'))
template = antsImageRead(
 Sys.glob(file.path(lesydata, 'template', 'ch2.nii.gz')))
lsm = lesymap(filenames, behavior, method = 'BMfast')
#> 07:56:04 Running LESYMAP 0.0.0.9221
#> 07:56:04 Checking a few things...
#> 07:56:04 Loading behavioral data... 131 scores found.
#> 07:56:04 Filenames as input, checking lesion values on 1st image...
#> 07:56:04 Searching voxels lesioned in >= 10% subjects... 326828 found
#> 07:56:35 Computing unique patches...
#> 07:56:51 Found 195102 patches in 326828 voxels - 1.7 times more voxels
#> 07:56:51 Using existing lesion matrix...  131x195102
#> 07:56:51 Checking matrix values are binary 0/1...
#> 07:56:51 Running analysis: BMfast ... 
#> 07:56:53 Correcting p-values: fdr ...
#> 07:56:53 Preparing images...
#> 07:56:54 Logging call details...
#> 07:56:54 Done! 49.7 secs
plot(template, lsm$stat.img, window.overlay = range(lsm$stat.img))
#> NULL

if (FALSE) {
# Same analysis with SCCAN
lsm = lesymap(filenames, behavior, method = 'sccan',
sparseness=0.045, validateSparseness=FALSE)
plot(template, lsm$stat.img, window.overlay = range(lsm$stat.img))
save.lesymap(lsm, saveDir='/home/dp/Desktop/SCCANresults')
}