Improve Brain_Data.randomise

[ejolly]

Currently Brain_Data.randomise has a naive implementation that's a bit slow and requires manually thresholding results afterwards. FSL's randomise on the other hand automatically computes TFCE maps or voxel-wise FWER maps. Ignoring TFCE for now, it would be nice to have a comparable FWER computation.

In my reading of Nichols & Holmes (2001), randomise computes FWER while it does permutation, by keeping track of the maximum t-stat over the whole brain for of each permutation iteration and then thresholding all voxels that exceed the Nth t-stat of that distribution according to the formula N = floor(alpha * num_perm) + 1. All voxels that are below that threshold don't survive correction at FWER alpha.

For example for 1000 permutations at a FWER alpha of .05 we would use the (.05 * 1000) + 1 = 51st largest t-stat from the distribution of maximal permuted t-stats as the critical cutoff for thresholding all voxels.

I have an implementation on the randomise branch of my fork that is reasonably fast through parallelization (~ randomise is 1.25-2x faster because C). However, it currently gives more conservative results relative to calling randomise directly through nipype.

All the action is in the new onesamp_max_stat_permutation function under nltools.stats here and I've documented each step reasonably well if someone could help out to see what might be going wrong.

You can call the function directly, but more conveniently you'd want to call it using the new Brain_Data.randomise_fwer() method noted here

[ljchang]

Have you tried using Numba to see if that helps speed up the permutations?

I'll take a look at the maximum t-stat function in the next few weeks while I prep for the fMRI class.

[ejolly]

Haven't played with numba because this implementation was fast enough, but maybe after we figure out discrepant results with FSL we can add it

[markallenthornton]

What's the current status of these efforts? The links to branches in Eshin's original post all return 404, so I assume they've moved to a different branch but I'm not sure where. We'd love to use a cluster-threshold or TFCE-based maximal statistic permutation method - happy to help with implementation, if I can!

[ljchang]

I had started prototyping a TFCE implementation awhile back based on @markallenthornton's matlab version, but never finished. I would love to have both randomise and TFCE implemented if anyone has time to work on this. My branch is so old, it's probably not worth looking at, but here is a gist with the key parts. No recollection of where i was at, and it may make more sense to start from scratch.

We currently are only supporting FDR correction at the moment and after looking at several implementations of cluster correction, I decided it wasn't worth the effort given the payoff.

@ejolly : had some amazing speedups with numba on a different toolbox we have been developing, so perhaps he might have new insights if it would be worth it here. There are also some other parallelization and GPU options we could explore.