Meta-analytic algorithms

[tsalo]

Here's a list of the methods we plan to support and their current statuses in my branch:

• Kernel-based coordinate-based meta-analyses:
  • ALE
  • SCALE
  • MKDA Density Analysis
  • MKDA Chi2 Analysis (Neurosynth speed-optimized version)
  • MKDA Chi2 Analysis with empirical null
  • KDA
• Model-based coordinate-based meta-analyses:
  • BHICP Available implementation only in C++
  • HPGRF/BHPGM Available implementation only in C++
  • SBLFR
  • SBR Cannot currently be run on a full brain. No source code available.
• Image-based meta-analyses:
  • Fisher's
  • Stouffer's FFX
  • Stouffer's RFX with theoretical null
  • Stouffer's RFX with empirical null
  • Weighted Stouffer's
  • FFX GLM
  • RFX GLM with theoretical null
  • RFX GLM with empirical null
  • MFX GLM

[tyarkoni]

Nice! Sorry I haven't had time to keep up with all the changes; will try to get to it later this week.

[tyarkoni]

Also looping @nicholst in in case he can think of other methods that are relatively straightforward to implement and are missing from the list. Tom, NiMARE is the better-named successor to the "PyCIBMA" package we proposed in our grant. @tsalo is a grad student working with Angie who's taking the lead on developing this.

[tsalo]

No pressure. I just wanted to set up a task list for milestone purposes (and also to surreptitiously ask for help on a couple of these algorithms).

I can get SCALE working soon. Just need to convert my code from pyale to work with NiMARE. I can also figure out the MKDA Chi2 empirical null from Dr. Wager's documentation. No idea how KDA works, but I'll start looking into it. The ALE and MKDA code that I do have working is probably a lot slower than it could be, because I'm masking and unmasking the data with nilearn much more than I probably need to.

I am stuck on the empirical nulls for Stouffer's RFX and RFX GLM. These aren't implemented in the IBMA MATLAB toolbox yet, so I did what I could from text in Maumet & Nichols (2016). I'm fairly certain I got something wrong.

Oh, and I forgot to add the model-based CBMA methods to the task list. Adding those now.

[tsalo]

There is a MATLAB package for SBLFR, so I should be able to translate that to Python with some time. The implementations available for BHICP and HPGRF/BHPGM are in C++, which is very cool, but unfortunately means I won't be able to translate them to Python any time soon. SBR, according to Samartsidis, Montagna, Nichols, & Johnson (2017), cannot currently be run on a full brain, so I'm going to skip it for now.

It seems like the best way to run an MFX GLM is through FLAME, so I was thinking of writing a short Nipype workflow for the NiMARE implementation.

[tsalo]

SCALE and MKDA Chi2 analysis with FWE correction seem to be working as of c237ff4, although when I tried to test them on the full Neurosynth database my computer ran out of application memory... They're also pretty slow in general, so there's still a fair amount of optimization to be done, I'm sure.

[tsalo]

I think I fixed the IBMA empirical nulls in 96aa5d2. Hopefully KDA is now working as well.

[tsalo]

I ran all of the meta-analyses on the 21 pain studies dataset. The results, thresholded at p<0.05 with a variety of corrections, are shown in the following image. The MFX GLM used FSL's cluster, with a cluster-defining threshold of p<0.01. RFX GLM, Z MFX, FFX GLM, Fisher's, Stouffer's, and Weighted Stouffer's all use Bonferroni correction. Contrast Permutation and Z Permutation use FDR correction because Bonferroni correction doesn't really work for the empirical null distributions. The MKDA and KDA density analyses use voxel-level FWE correction, and ALE and SCALE both use uncorrected thresholds. SCALE can't use multiple comparisons correction, but I chose to only plot the uncorrected map for ALE because the voxel-level and cluster-level FWE corrected maps were empty.

The only thing that seems off at this point (barring more quantitative checks) is the MKDA Chi2 analysis with FWE. I'm pretty sure I messed something up with that one.

[tyarkoni]

Very nice! The chi2 FWE map is definitely messed up; I'll try to look at that soon.

Have you tried to reproduce these results with @cmaumet's Matlab IBMA package? For at least the algorithms implemented in both packages, the results should presumably look (almost) identical, and quantitative comparison of the maps would provide a nice sanity check.

[tsalo]

Good thinking. I was able to compare the Fisher's, Stouffer's FFX, Stouffer's RFX, Weighted Stouffer's, and FFX GLM IBMAs, and even caught a couple of mistakes I made in the Fisher's and FFX GLM implementations. The RFX GLM implementation in the MATLAB toolbox seems kind of complicated (i.e., I won't be able to use it without completely setting up the toolbox), so I haven't tested that yet. The MFX GLM, Contrast Permutation, and Z Permutation IBMAs are not implemented in the MATLAB package, though there is Python code for the MFX GLM in the nidmresults-paper repo.

I converted the MFX GLM Python code to use Nipype (still need to make it an actual workflow), which is what is in the figure. A couple of things that I changed for that one are that I didn't rescale the FSL maps, nor did I take the square roots of the varcope images. I wasn't sure if either step was necessary. The scaling didn't affect the results at all and it seemed like the varcopes were already in the right units in the released datasets.

The only weird thing in the IBMA toolbox code was that the tests all appear to be one-tailed (or at least the ones I've checked). I assumed that these tests (e.g., the Stouffer's t-test on z-values) should be two-tailed and that's how I had them implemented. I can change that, but I'm not sure it makes sense to do so.

[tsalo]

I think I fixed the MKDA Chi2 code. Part of the problem is that I was using voxelwise null distributions using each value generated by the permutations rather than a general null distribution using the max value from each permutation. The other problem is that I wasn't transforming p-values to z-values correctly.

However, once I edited all of the IBMAs to use two-tailed tests, it became clear that everything in the FFX GLM was significant. Really not sure why. Another reason for the differences between the two figures is that I changed the threshold from 0.001 to 0.05 (still FDR- or FWE-corrected) to match what Dr. Maumet reported in one of her presentations on IBMAs.

[tsalo]

@cmaumet Do you think using two-sided tests for the IBMAs is okay? Also, do you have any idea why everything is significant in the FFX GLM?

[tyarkoni]

I can't speak to the first question (though I don't see why it wouldn't be okay), but assuming the z or t scores on the colorbars are correct for all analyses, it sure looks like there's a bug somewhere. The general idea that all or nearly all voxels would be significant in an FFX analysis is not intrinsically problematic (since the null is surely false everywhere in the brain, so at the limit, if you had enough data, all of these maps would be significant everywhere). But it's pretty weird that there's so little spread in the values, especially given that the max is not very large. I'll take a look at the code.

BTW, I think you can go ahead and merge your branch into master--it'll be easier to work with that way.

[tsalo]

Thanks for taking a look. @cmaumet was kind of enough to offer her expertise as well, so hopefully someone will be able to figure out what's happening.

It's possible that everything really should be significant, although there are only 21 studies, so I doubt it. My guess is that there's a bug, although the Python code (except when two-sided inference is turned on) produces equivalent results to the MATLAB version, which shows in the last two cells of this notebook. Well, there is one odd voxel with very different results, but other than that the two results are equivalent.

I'll merge into master now to make things easier.

[cmaumet]

@tsalo: thanks for pinging me! For one-tailed versus two-tailed, my understanding is that we could do both. Because FSL and SPM both use one-tailed tests, this is also the way we have implemented the meta-analyses tests in the IBMA toolbox.

Which method did you use to compute the degrees of freedom for FFX GLM? In our analyses, we had used a wrapper around FSL's FFX model, maybe that's also something you could also use via nipype?

You are probably already aware of this but the results shown in the following figure were obtained using the same dataset of 21 pain studies with p<0.05 FDR-corrected (one-tailed): https://www.frontiersin.org/10.3389/conf.fninf.2014.18.00025/event_abstract.

The RFX GLM implementation in the MATLAB toolbox seems kind of complicated (i.e., I won't be able to use it without completely setting up the toolbox), so I haven't tested that yet.

As far as I remember, for RFX GLM we reused SPM's implementation of one-sample t-tests which you could maybe also do through nipype?

The MFX GLM, Contrast Permutation, and Z Permutation IBMAs are not implemented in the MATLAB package, though there is Python code for the MFX GLM in the nidmresults-paper repo.
I converted the MFX GLM Python code to use Nipype (still need to make it an actual workflow), which is what is in the figure. A couple of things that I changed for that one are that I didn't rescale the FSL maps, nor did I take the square roots of the varcope images. I wasn't sure if either step was necessary. The scaling didn't affect the results at all and it seemed like the varcopes were already in the right units in the released datasets.

For MFX GLM, we reused FSL's flameo function with 'flame1' option. Not computing the rescaling will probably not affect the results much (this is to compensate for different solutions used by SPM and FSL to scale the data, leading to different unit scales in the contrast maps). But if you are also using flameo then the contrast variance maps should be passed as inputs rather than standard error maps which is why we squared (-sqr) the standard error maps form NIDM-Results. Does that make sense or did I miss something?

This project looks great!! Please let me know if you think I can provide further input. Thanks!

[tsalo]

@cmaumet Thank you for your help! Your responses really clarified some things for me.

For one-tailed versus two-tailed, my understanding is that we could do both.

Awesome! I’ll add an argument to the IBMA functions to switch between one- and two-tailed versions.

Which method did you use to compute the degrees of freedom for FFX GLM? In our analyses, we had used a wrapper around FSL’s FFX model, maybe that’s also something you could also use via nipype?

I translated the MATLAB code from the IBMA toolbox and adapted it to be two-tailed, so I didn't change anything in the degrees of freedom calculation. It's definitely a good idea to just use FSL, especially since we already have everything set up for the MFX GLM function. I’ll open a PR with a new version of the FFX GLM that is pretty much the same as the MFX GLM with FSL’s flameo, but using fixed effects (fe) run mode.

As far as I remember, for RFX GLM we reused SPM’s implementation of one-sample t-tests which you could maybe also do through nipype?

I’m currently using scipy.stats.ttest_1samp, but I can compare the results to SPM output.

For MFX GLM, we reused FSL’s flameo function with ‘flame1’ option. Not computing the rescaling will probably not affect the results much (this is to compensate for different solutions used by SPM and FSL to scale the data, leading to different unit scales in the contrast maps). But if you are also using flameo then the contrast variance maps should be passed as inputs rather than standard error maps which is why we squared (-sqr) the standard error maps form NIDM-Results. Does that make sense or did I miss something?

Ohhh, I didn’t realize that flameo did use a different unit than provided by the NIDM Results. That makes sense. I’ll add the squaring back into the MFX GLM/FFX GLM code.

[tsalo]

It's unclear if we'll ever take a crack at SDM, but as far as model-based CBMA methods go, per a discussion with @nicholst it's likely that one frequentist approximation of a Bayesian approach will be developed and implemented, as well as possibly one of the three above models. Closing this issue now that most methods have been implemented.