Critical threshold calculation in statistical non-parameric mapping

[zof1985]

Hi Todd,

In your paper (Pataky et al. 2015) you described the implementation of the Nichols and Holmes (2002) statistical non-parametric mapping approach to 1D data.

I just want to be sure to have understood the procedure correctly, thus I think that my questions might be of interest also for other spm1d users.

Looking at the equation A.3: P(t > u) = (number of permuted t value equal or exceeding u) / (number of permutation)

1. Is the value of "u" calculated by solving this equation after having set it equal to a given alpha level?
2. Is the "t" in the left side of the equation the max t value resulting from all the permutation tests?

Best, Luca.

[0todd0000]

Hi Luca,

1. Is the value of "u" calculated by solving this equation after having set it equal to a given alpha level?

Yes. In fact, u can be calculated quite easily as the "100 x (1-alpha)th" percentile of the distribution. Imagine that you run 10,000 permutations, compute a test statistic value for each permutation, and then assemble those 10,000 test statistic values in an array called Z. If alpha=0.05 then u is simply:

import numpy as np
u = np.percentile(Z, 95)

2. Is the "t" in the left side of the equation the max t value resulting from all the permutation tests?

Yes, for 1D data that's correct; this equation represents the probability that the maximum t value will exceed u.

In spm1d it's implemented as follows:

• For each permutation compute the t trajectory and store its maximum value. (If there are 10,000 permutations there will be 10,000 maximum t values.)
• After repeating for a large number of permutations, use the resulting distribution of maximum t values to compute u with np.percentile as indicated above.

If the data are normally distributed then the u value computed in this manner should be very close to the u value computed using random field theory.

Todd

[zof1985]

Hi Todd,

Thank you for your quick answer. Now I have just one more question. About the point 1, if I'm looking for a two-tailed test. Should I divide alpha by 2 before computing the percentile? i. e. something like:

import numpy as np

alpha = 0.05
p = 1 -  (alpha/2 if test == 'two-tailed' else alpha)
u = np.percentile(Z, 100 * p)

Many thanks, Luca.

[0todd0000]

Hi Luca,

Yes, that's correct.

If you're using spm1d 0.4.1 you'll find nearly the exact same command on Line 105 of ./spm1d/stats/nonparam/_snpm.py:

>>>  alpha0 = 0.5*alpha if two_tailed else alpha

Todd

[zof1985]

Many thanks, Luca.