What's new in v0.0.0.9201
Version 0.0.0.9201 focused mostly on improving the SCCAN method, and is described below. Other changes in more recent versions are described in the NEWS file.
In earlier LESYMAP versions, SCCAN weights were unidirectional, they could not capture both negative and positive relationships with behavior. While this approach works fine for many analyses, the unidirectional relationship does not always reflect lesion-to-symptom relationships. On one side, there is literature showing that partial lesion or preservation of some regions may exacerbate cognitive deficits. On the other side, researchers have observed dual sided effects in mass-univariate analyses.
To address this issue, we updated LESYMAP to allow both positive and negative voxel weights.
To allow a full capture of all brain-behavior relationships changes were made to:
- allow the search of both positive and negative voxel weights
- return a map with correctly oriented voxel weights with respect to behavior.
Tests show that the new approach works well both with old simulations and with new simulations of multiple brain areas with opposite relationships with behavior (see example below).
The stat.img map produced from SCCAN will now contain correctly oriented voxel weights with respect to behavior. Negative weights are inversely related to behavioral scores (1-lesion = lower behavioral score), and positive weights are directly related to behavioral scores (0-healthy = lower behavioral score). While stat.img contains normalized voxels weights between -1 and 1, raw voxel weights can be found in rawWeights.img.
Be default, LESYMAP will now search for positive and negative voxel weights. This breaks the consistency with the previous behavior.
There is a new argument called directionalSCCAN, simply set directionalSCCAN=FALSE to obtain results similar to the older versions of LESYMAP.
When running SCCAN, a positive sparseness (i.e., 0.02) forces all the voxels to have weights on one side of the spectrum, i.e., all positive or all negative. In SCCAN terminology, the weights are unsigned. To allow signed voxel weights, sparseness must be negative. Thus, we expanded the sparseness search to go from -0.9 to 0.9 in LESYMAP v0.0.0.9201 .
You may think that a negative sparseness doesn't force voxels to have mixed signs, it just allows them. Therefore the solution might still be good since all voxels can sit on one side of the spectrum if they want to. This is true in many cases, but in practice we have noticed that negative sparseness may cause instability. For example, in a small dataset (N=27) we noticed that negative sparseness yields CVcorr=0.33 (p=0.09) while positive sparsness yields CVcorr=0.437 (p=0.02). The low predictive accuracy with negative sparseness might have derived from the SCCAN iterations being unnecessarily confused from voxels with an unusual relationship with behavior. Indeed this dataset contained both positive and negative relationships with behavior (found also with univariate analyses), but only the negative weights were truly predictive of behavior, while the positive weights were not. Thus, it is important to allow the optimization algorithm to search on both signs of sparseness for the point at which the CV predictive accuracy is at its maximum.
LESYMAP uses the optimize() function in R to search for optimal sparseness. This function identifies the best solution if the distribution is unimodal, i.e., if there is a single peak of sparseness that yields maximum CV correlation. The new approach of searching positive and negative sparseness likely forms a bimodal distribution, i.e., there are two peaks, one on the positive side and one on the negative side. In these cases, there is a risk that optimization is stuck in a local minimum without finding the global minimum, which might be on the other side. Our initial checks show that, despite this issue, the optimization converges systematically on the correct side.
If you want to check this issue on your own data, run LESYMAP twice: one time with
lowerSparseness=-0.9; upperSparseness=-0.005
and one time with
lowerSparseness=0.005; upperSparseness=0.9.
Check if the side with the maximum CV correlation corresponds to what LESYMAP found automatically by itself. If you see discrepancies, let us know so we can consider integrating dual side search in LESYMAP. Currently, it doesn't seem worth to integrate such a lengthy procedure as a standard approach in LESYMAP.
Note, this is an experimental feature.
The standard approach of LESYMAP is to run SCCAN and remove voxels with weights <10% of the peak weight. However, SCCAN can do this approach by itself in a sort of "internal thresholding" after each iteration. This process is much faster than standard SCCAN, and often produces similar results. However, we have noticed that maxBased solutions may diverge considerably from the standard solution found in LESYMAP, so this feature is there just for experimental approaches, not to be used routinely. From a runtime perspective, maxBased SCCAN is 3-20 times faster than regular SCCAN. Note that the two approaches yield different optimized sparseness values, so maxBased cannot be used to accelerate the search for optimal sparseness either.
A reproducible example with two opposite brain-behavior relationships:
library(LESYMAP)
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'))
parcellation = antsImageRead(
Sys.glob(file.path(lesydata, 'template', 'Parcellation_403areas.nii.gz')))
lesions=imageFileNames2ImageList(filenames)
# simulate regions 140 and 75
simulatedRegions = c(140, 75)
set.seed(123456)
sims = simulateBehavior(lesions, errorWeight = 0.5,
parcellation,
label = simulatedRegions)
# invert one of the region simulations and combine both
# with an AND simulation
sims$behavLoad[,2] = -sims$behavLoad[,2]
simBehavior = rowMeans(sims$behavLoad)
# create a folder to save results
dir.create('/DATA/dorian/lsmLesySimulations')
setwd('/DATA/dorian/lsmLesySimulations')
A search for sparseness both positive and negative converges correctly on the negative side:
lsmLesyDualSided = lesymap(lesions, simBehavior,
method='sccan',
directionalSCCAN=T,
upperSparseness=0.9,
lowerSparseness=-0.9)
23:15:42 Running LESYMAP 0.0.0.9201
23:15:42 Checking a few things...
23:15:43 SCCAN method: ignoring patch, nperm, and multiple comparison...
23:15:43 Searching voxels lesioned >= 10% subjects...326828 found
23:15:43 noPatch true - Patches will not be used...
23:15:43 Computing lesion matrix... 131x326828
23:15:49 Running analysis: sccan ...
Searching for optimal sparseness:
lower/upper bound: -0.9 / 0.9
cvRepetitions: 3
nFolds: 4
sparsenessPenalty: 0.03
optim tolerance: 0.03
23:15:53 Checking sparseness -0.212 ... CV correlation 0.739 (0.785) (cost=0.267)
23:50:33 Checking sparseness 0.212 ... CV correlation 0.456 (0.587) (cost=0.550)
00:03:09 Checking sparseness -0.475 ... CV correlation 0.731 (0.767) (cost=0.284)
00:35:29 Checking sparseness -0.314 ... CV correlation 0.737 (0.778) (cost=0.273)
01:08:58 Checking sparseness -0.05 ... CV correlation 0.727 (0.789) (cost=0.274)
01:44:24 Checking sparseness -0.189 ... CV correlation 0.739 (0.786) (cost=0.266)
02:18:59 Checking sparseness -0.151 ... CV correlation 0.739 (0.788) (cost=0.265)
02:54:11 Checking sparseness -0.112 ... CV correlation 0.738 (0.791) (cost=0.265)
03:29:47 Checking sparseness -0.102 ... CV correlation 0.739 (0.791) (cost=0.264)
04:05:17 Checking sparseness -0.082 ... CV correlation 0.734 (0.791) (cost=0.268)
04:40:45 Checking sparseness -0.102 ... CV correlation 0.739 (0.791) (cost=0.264)
Found optimal sparsenes -0.102 (CV corr=0.739 p=0.00000000000000000000000717)
Calling SCCAN with:
Components: 1
Use ranks: 1
Sparseness: -0.102
Cluster threshold: 150
Smooth sigma: 0.4
Iterations: 20
maxBased: FALSE
directionalSCCAN: TRUE
Permutations: 0
05:19:42 Preparing images...
05:19:42 Logging call details...
05:19:42 Done! 6.1 hours
save.lesymap(lsmLesyDualSided, 'lsmLesyDualSided')
Both areas are correctly found, and the relationship with behavior is retained:
(blue=negative correlation with behavior, red=positive correlation with behavior)
(dark areas are the simulated regions)
Now the old SCCAN approach, searching only for positive sparseness:
lsmLesyPosSided = lesymap(lesions, simBehavior,
method='sccan',
directionalSCCAN=T,
upperSparseness=0.9,
lowerSparseness=0.005)
23:16:08 Running LESYMAP 0.0.0.9201
23:16:08 Checking a few things...
23:16:08 SCCAN method: ignoring patch, nperm, and multiple comparison...
23:16:08 Searching voxels lesioned >= 10% subjects...326828 found
23:16:09 noPatch true - Patches will not be used...
23:16:09 Computing lesion matrix... 131x326828
23:16:15 Running analysis: sccan ...
Searching for optimal sparseness:
lower/upper bound: 0.005 / 0.9
cvRepetitions: 3
nFolds: 4
sparsenessPenalty: 0.03
optim tolerance: 0.03
23:16:19 Checking sparseness 0.347 ... CV correlation 0.438 (0.566) (cost=0.573)
23:38:10 Checking sparseness 0.558 ... CV correlation 0.457 (0.562) (cost=0.560)
00:02:03 Checking sparseness 0.689 ... CV correlation 0.457 (0.556) (cost=0.564)
00:26:56 Checking sparseness 0.568 ... CV correlation 0.457 (0.562) (cost=0.560)
00:51:04 Checking sparseness 0.477 ... CV correlation 0.460 (0.563) (cost=0.554)
01:12:02 Checking sparseness 0.428 ... CV correlation 0.464 (0.556) (cost=0.549)
01:31:11 Checking sparseness 0.397 ... CV correlation 0.462 (0.557) (cost=0.550)
01:49:54 Checking sparseness 0.418 ... CV correlation 0.463 (0.556) (cost=0.549)
02:09:10 Checking sparseness 0.447 ... CV correlation 0.464 (0.556) (cost=0.550)
02:28:26 Checking sparseness 0.428 ... CV correlation 0.464 (0.556) (cost=0.549)
Found optimal sparsenes 0.428 (CV corr=0.464 p=0.0000000246)
Calling SCCAN with:
Components: 1
Use ranks: 1
Sparseness: 0.428
Cluster threshold: 150
Smooth sigma: 0.4
Iterations: 20
maxBased: FALSE
directionalSCCAN: TRUE
Permutations: 0
02:49:02 Preparing images...
02:49:02 Logging call details...
02:49:02 Done! 3.5 hours
save.lesymap(lsmLesyPosSided, 'lsmLesyPosSided')
Only the anterior region is found, the posterior region is missed completely:
Switching to maxBased SCCAN with positive sparseness:
lsmLesyPosSidedMaxbased = lesymap(lesions, simBehavior,
method='sccan',
maxBased=T,
directionalSCCAN=T,
upperSparseness=0.9,
lowerSparseness=0.005)
02:49:04 Running LESYMAP 0.0.0.9201
02:49:04 Checking a few things...
02:49:04 SCCAN method: ignoring patch, nperm, and multiple comparison...
02:49:04 Searching voxels lesioned >= 10% subjects...326828 found
02:49:04 noPatch true - Patches will not be used...
02:49:04 Computing lesion matrix... 131x326828
02:49:10 Running analysis: sccan ...
Searching for optimal sparseness:
lower/upper bound: 0.005 / 0.9
cvRepetitions: 3
nFolds: 4
sparsenessPenalty: 0.03
optim tolerance: 0.03
02:49:13 Checking sparseness 0.347 ... CV correlation 0.493 (0.605) (cost=0.517)
02:56:26 Checking sparseness 0.558 ... CV correlation 0.496 (0.600) (cost=0.520)
03:03:37 Checking sparseness 0.216 ... CV correlation 0.508 (0.603) (cost=0.498)
03:10:47 Checking sparseness 0.136 ... CV correlation 0.488 (0.602) (cost=0.516)
03:17:55 Checking sparseness 0.24 ... CV correlation 0.515 (0.603) (cost=0.493)
03:25:03 Checking sparseness 0.281 ... CV correlation 0.508 (0.605) (cost=0.501)
03:32:11 Checking sparseness 0.25 ... CV correlation 0.516 (0.603) (cost=0.492)
03:39:20 Checking sparseness 0.26 ... CV correlation 0.511 (0.606) (cost=0.497)
03:46:26 Checking sparseness 0.25 ... CV correlation 0.516 (0.603) (cost=0.492)
Found optimal sparsenes 0.25 (CV corr=0.516 p=0.000000000291)
Calling SCCAN with:
Components: 1
Use ranks: 1
Sparseness: 0.25
Cluster threshold: 150
Smooth sigma: 0.4
Iterations: 20
maxBased: TRUE
directionalSCCAN: TRUE
Permutations: 0
03:54:20 Preparing images...
03:54:21 Logging call details...
03:54:21 Done! 1.1 hours
save.lesymap(lsmLesyPosSidedMaxbased, 'lsmLesyPosSidedMaxbased')
This time the anterior region is missed and the posterior is found.
Conclusion: unsigned voxel weights may miss complex brain-behavior relationships. This is resolved using directional SCCAN to obtained signed voxel weights.
Note:
Although we use the term SCCAN, the approach currently maps a single behavioral score, therefore can be considered a sparse regression approach.