0todd0000 / spm1dmatlab

One-Dimensional Statistical Parametric Mapping in Matlab.
GNU General Public License v3.0
28 stars 13 forks source link

Non-parametric ANOVA2rm #45

Closed NicolasFlores14 closed 7 years ago

NicolasFlores14 commented 7 years ago

Dear Todd,

I have a question about the non-parametric procedure for ANOVA2rm.

I have performed a non-parametric ANOVA2rm, which has shown 2 clusters exceeding the significant threshold for a given factor (factor A for example). However, these 2 clusters have exactly the same p value. Then, to check whether this kind of output was systematic, I performed another non-parametric ANOVA2rm for another variable, which showed again 2 clusters exceeding the significant threshold. As in the previous case, the p value was identical for both clusters.

Also, even if my experimental design didn't allow it, I performed a non-parametric paired t-test on the factor A. With this analysis, 2 clusters exceeded the significant threshold, and with 2 different p values.

Do you know if the problem comes from my data (or my inputs in the function) or from the non-parametric ANOVA2rm?

Thanks a lot for your help, Best regards.

Nicolas Flores

0todd0000 commented 7 years ago

Hi Nicolas,

I don't think it's necessarily a problem, I think it's just a limitation of spm1d's non-parametric permutation tests. For these tests p values mean the following: p percent of permutations produced at least one cluster at least as large as the given cluster. That interpretation implies that different clusters can have the same p value, even when those clusters look quite different geometrically.

More specifically, there are two relevant factors:

  1. Minimum p value. If there are N permutations then the minimum possible p value is 1/N. Even if clusters are vastly different sizes they will have the same p value if

  2. Cluster metric. By default a simple "maximum extent" metric is used, which is the maximum number of supra-threshold nodes in a single cluster for each permutation. It is possible to use two other metrics: maximum cluster height and maximum cluster integral. These can be specified using the "cluster_metric" keyword argument as follows:

snpm       = spm1d.stats.nonparam.ttest(y, 0);
snpmi      = snpm.inference(alpha, 'iterations', 1000, 'cluster_metric',  'MaxClusterExtent');
snpmi      = snpm.inference(alpha, 'iterations', 1000, 'cluster_metric',  'MaxClusterHeight');
snpmi      = snpm.inference(alpha, 'iterations', 1000, 'cluster_metric',  'MaxClusterIntegral');

However, changing the number of permutations and/or the cluster metric does not guarantee that the p values will be different for different clusters. If just one permutation yields multiple large suprathreshold clusters, or if multiple permutations happen to yield clusters with identical metrics, then clusters can have identical p values even if they have quite different sizes.

Todd

NicolasFlores14 commented 7 years ago

Hello Todd,

Thanks for your explanations, that is clear now.

Best regards.

Nicolas

NicolasFlores14 commented 7 years ago

Hello Todd,

I reopen this issue for further questions about the permutation procedure. I would like to use it but I also would like to be able to well explain it.

I read your paper "Zero- vs. one-dimensional, parametric vs. non-parametric, and confidence interval vs. hypothesis testing procedures in one-dimensional biomechanical trajectory analysis" (2015) and some questions remain.

In the figure 4, you explained principles of the permutation. You mentioned that only one permutation exceeded the critical value t* (the original trajectory). You formed a second permutation PDF, and here, there were 2 bin categories: one between 0 and 1 (about 19 in Frequency), and one between 7 and 8 (about 1 in Frequency). What did represent the x-label? The number of suprathreshold clusters (question 1)?

If I understand, the value of the frequency was the value of the integral of the suprathreshold cluster, is it right (question 2)? And this value of frequency could be different according to the cluster metric (maximum cluster extent, or maximum cluster height, or maximum cluster integral) isn't it (question 3)?

If it is, what did represent the bin category between 7 and 8 (question 4)? Indeed, you mentioned that only the original trajectory exceeded the suprathreshold cluster, so if my summarise is correct, I don't understand what was the bin category between 7 and 8.

In the Appendix A, you illustrated permutations with the labeling (AAAAA-BBBBB) and associated values (1.14 1.21 1.25 1.43 1.57 - 1.37 1.52 1.61 1.74 1.54). You explained that each permutation resulted in a critical value t*. I wonder if the permutation pertains to the labelling and its associated value, such as BAAAA-ABBBB (1.37 1.21 1.25 1.43 1.57 - 1.14 1.52 1.61 1.74 1.54) (question 5)?

To finish, this is more personal for my study, how can I determine the number of permutation according to my experimental design (question 6)? There were 19 participants and 4 conditions, which combined 2 factors (A and B) with two levels each (0 and 1). 5 trials were performed for each condition.

Thanks in advance for your useful help, Best regards.

Nicolas

0todd0000 commented 7 years ago

Hi Nicolas,

(Question 1) In Fig.4 the "Primary PDF" and "Secondary PDF" are both histograms, where the x axis represents the 1D geometrical feature of interest upon which statistical inference is based. In the Primary PDF the feature is the maximum t value (across the whole 1D continuum). In the Secondary PDF the feature is the maximum cluster integral, where "clusters" are threshold-surviving "upcrossings". In other words, thresholding an nD continuum generally leaves a number of nD upcrossings, and the Secondary PDF describes probabilities associated with those upcrossing's geometry. In this case upcrossing integral was used, but other features are possible including extent and height. Incidentally, spm1d's parametric inference uses extent, for which relatively easy analytical probabilities exist.

(Question 2) Yes

(Question 3) Yes

(Question 4) In this case the upcrossing happened to have an integral value between 7 and 8. If you look at Fig.4d (left side) you'll see the upcrossing, and also that it has a height of about one and an extent of about ten; if it were a rectangle its integral would be about 1 x 10 = 10, but since it is not a rectangle its integral is slightly smaller than 10. Thus, rather than "bin category" think of the bins as forming a discrete (histogram) approximation to a continuous probability density. In this case there are very few discrete permutations, implying that the discrete approximation is probably not very similar to the true, continuous probability density. Increasing the number of permutations generally causes the histogram to converge to the continuous density. Also, when the data are normally distributed the discrete histogram approximation will converge to the analytical densities (i.e. Student's t distribution, F distribution, random field theory distributions, etc.)

(Question 5) Each permutation produces a t value, not a critical t value. The critical t value can only be computed after t values are computed for a large number of permutations. In other words, each permutation contributes a single t value to the ultimate t distribution. The critical t value is the 100 x (1 - alpha) th percentile of that distribution. Your interpretation of permutations is correct: each permutation involves assigning the observations to different groups then computing the t statistic.

(Question 6) Conventionally 10,000 or 100,000 permutations are used to provide a reasonable approximation to the underlying continuous probability density. If the study has a small number of observations, and the total possible number of permutations is less than 10,000 then all permutations should be used (like in Fig.4). I'd suggest selecting 10,000 as a default, then repeating the analysis a number of times with different random number generator states, which can be set using MATLAB's rng function:

>>   rng(0)
>>   rand(1,3)    % random vector #1
>>   rand(1,3)    % different from random vector #1
>>   rng(0)
>>   rand(1,3)    % same as random vector #1

Setting the rng value before starting permutation will make the permutation result precisely repeatable. If you repeat permutation for a few different rng values you should notice that the critical t or F value changes slightly, but also that those changes probably do not affect the results qualitatively. If the results don't change qualitatively you can be reasonably confident that 10,000 iterations is sufficient.

Todd

NicolasFlores14 commented 7 years ago

Hello Todd,

Thanks for all your answers.

I performed permutations on my data (the 3 components of the GRF) with 10,000 iterations. However, as I notified you in the first post, the p value is always 1/10000=0.0001 no matter what the GRF component (3 in my example), the number of clusters (a total of 8 clusters on the thrice components with p=0.0001 for each one in my example) and the cluster size (very different sizes in my example).

Maybe I did not understand yet how works the permuation... but I think it is strange to always obtain the same p value even if the cluster is very very small.

In comparison on one very small cluster on a GRF component, when I have performed the parametric procedure (because the non-parametric toolbox was not available), this cluster had a p value p=0.048.

Does exist cases where the p value of the cluster is not equal to the minimum p=1/N (with N the number of permutations)?

Thanks again for your help, Best regards.

Nicolas

0todd0000 commented 7 years ago

Hi Nicolas,

Yes, there should be cases where p values are not 1/N. Please see the example below, which is a modification of: ./spm1d/examples/nonparam/1d/ex1d_ttest2.m In this example you should see two clusters and two p values: one is 0.013 and one is 0.026, and both are larger than 1/N (which is 1/924 = 0.0011).

clear;  clc

% Load data:
dataset = spm1d.data.uv1d.t2.SimulatedTwoLocalMax();
[yA,yB] = deal(dataset.YA, dataset.YB);

% Conduct non-parametric test:
rng(0)
alpha      = 0.05;
two_tailed = false;
iterations = -1;
snpm       = spm1d.stats.nonparam.ttest2(yA, yB);
snpmi      = snpm.inference(alpha, 'two_tailed', two_tailed, 'iterations', iterations);

% Plot:
close all
snpmi.plot();
snpmi.plot_threshold_label();
snpmi.plot_p_values();

Please confirm that you see the same results. If you see the same results then I think the p value computations are OK.

To help explain the results in your dataset please send:

  1. A screenshot of the result showing a low p value for a small cluster.
  2. The output from "disp( snpmi )" (i.e. similar to the output below).
SnPM{t} inference (1D)
              z: [1×101 double]
    nPermUnique: 924
    nPermActual: 924
          alpha: 0.0500
          zstar: 3.5782
       h0reject: 1
              p: [0.0130 0.0260]

Todd

NicolasFlores14 commented 7 years ago

Hello Todd,

I confirm that I get the same results.

The small cluster below: plot

and the output from "disp": disp

Thanks for your help, Best regards.

Nicolas

0todd0000 commented 7 years ago

Hi Nicolas,

In this case I suspect the reason for the low p value is that the original permutation (i.e. the original, unpermuted data) is the only permutation to have produced a cluster. The p value will only be greater than 1/N if at least two permutations produce a supra-threshold cluster. Thus a different dataset with a larger cluster can actually have a higher p value if more than one permutation produces a similar or larger cluster.

It is worth noting that, unlike parametric results, permutation-based results are dataset specific, and not necessarily generalizable across different datasets. This is because different datasets can have different distributions (e.g. normal, skewed left, skewed right, bimodal, etc.). Parametric results are perfectly generalizable across all datasets (and p values decrease as cluster sizes increase) because the parametric approach assumes that all datasets have the identical (normal) distribution.

Todd