bug in statistics_montecarlo?
Maren Grigutsch
grigu at CBS.MPG.DE
Fri Jul 6 15:27:14 CEST 2007
Eric Maris wrote:
>Dear Maren,
>
>
>Most likely, there will be a bug in statistics_montecarlo somewhere, but I
>do not think that your numerical demonstration is evidence for such a bug.
>
>The following points are relevant in this respect:
>
>
> 1. For a permutation test to be correct (in the sense that it controls the
> false alarm rate), it is NOT required that it produces the same p-value as
> the corresponding parametric test (i.e., using the same test statistic and
> calculating the p-value under auxiliary assumption of normality).
>
>
>
Dear Eric,
Many thanks for your fast response! You may be right that the data I
sent you cannot be regarded as an unequivocal evidence for a bug in
statistics_montecarlo.
What I had in mind and wanted to demonstrate to you was the following:
The dependent samples t-test gives accurate P-values if the sampling
distribution of the difference of means is roughly normal, as was the
case in my example. The permutation distribution of the dependent
samples t-statistic approximates the sampling distribution of t-values
under the null hypothesis. For my example with normally distributed
differences and 5000 permutations, the permutation distribution should
closely resemble the theoretical t distribution underlying the analytic
t-test. Thus, I had expected te get roughly the same p-values from both
the analytic and the permutation t-test. However, when I repeated my
simulations many times, the p-values obtained from
statistics_montecarlo, using a two-tailed test, were always smaller by a
factor of about 1/2 as compared to the p-values obtained from
statistics_analytic or Matlab's ttest function. In contrast, when I
applied a one-sided test the P-values were roughly the same between the
different functions.
>
>2. To show that there is a bug in the statistics_montecarlo, you must show
>that its permutation p-values do not have the properties they should have
>according to the theory behind permutation tests (see, Maris and Oostenveld,
>in the most recent issue of the Journal of Neuroscience Methods). For
>instance, you could do the following:
>
>1. Take a test statistic, e.g., the dependent samples t-statistic.
>
>2. Generate N pairs of observations from a bivariate distribution with
>identical marginal distributions (this is the null hypothesis of a
>permutation test for a within-units study).
>
>3. Calculate the permutation p-value under the permutation distribution,
>which is obtained by randomly permuting the elements of the N pairs.
>
>4. Repeat steps 2-3 a large number of times and calculate the proportion
>of times that it is less than 0.05.
>
>5. If the proportion in 4. is significantly different from 0.05 (the
>nomical alpha-level) you have shown that there is bug in
>statistics_montecarlo.
>
>
>
>
Well, this sounds very convincing. To follow your suggestions I ran the
following script (using only 500 permutations, for the sake of time):
---------
Nrep = 500;
Nperm = 500;
N = 50;
design = [
1:N 1:N % subject number
repmat(1,[1 N]) repmat(2,[1 N]) % condition number
];
cfg = [];
cfg.statistic = 'depsamplesT';
cfg.method = 'montecarlo';
cfg.correctm = 'no';
cfg.alpha = 0.05;
cfg.tail = 0;
cfg.uvar = 1; % "subject" is unit of observation (unit
factor)
cfg.ivar = 2; % "condition" is the independent variable
(repeated measures)
cfg.numrandomization = Nperm;
cfg.feedback = 'no';
prob = zeros(1,Nrep);
T = zeros(1,Nrep);
wb = waitbar(0,'Please wait...');
for k=1:Nrep,
x = randn(1,N);
y = randn(1,N);
[stat_montecarlo, cfg] = statistics_montecarlo(cfg,[x y],design);
prob(k) = stat_montecarlo.prob;
T(k) = stat_montecarlo.stat;
waitbar(k/Nrep,wb);
end
close(wb);
% should be approximately 0.05
sum(prob<0.05) / Nrep
-------
And here is the result:
ans =
0.1000
If I got it correctly, this means that the observed FA rate (0.1) is
twice the expected value of 0.05 (the nominal alpha level). So it still
seems to me that there is something wrong with statistics_montecarlo.
I'm just about reading the paper of yours mentioned above. In section
3.1.1 you describe cluster-based nonparametric statistical tests, and
you give the following statement: "for a two-sided test, we select test
statistics whose absolute value is larger than some threshold".
In my opinion, this is exactly what should be done in
statistics_montecarlo for a two-sided test: comparing the *absolute*
t-values.
>A simulation study of the type described here was performed by Maris,
>Schoffelen, & Fries (2007, also in the Journal of Neuroscience Methods) to
>show that a permutation test for coherence differences controls the false
>alarm rate whereas the corresponding parametric test does not.
>
>
Thank you for pointing me to that paper. I'm sure it will help me a lot
to improve my understanding of nonparametric statistical tests.
I'm looking forward to your response. Hope not to bother you too much.
But I think the issue has some impact and is worth to be thoroughly
discussed.
Best regards,
--Maren
>
>
>
>
>Kind regards,
>
>
>
>
>
>Eric Maris
>
>
>
>dr. Eric Maris
>
>NICI/Biological Psychology and
>
>F.C. Donders Center for Cognitive NeuroImaging
>
>University of Nijmegen
>
>P.O. Box 9104
>
>6500 HE Nijmegen
>
>The Netherlands
>
>T:+31 24 3612651 (NICI)
>
>T:+31 24 3610754 (FCDC)
>
>F:+31 24 3616066 (NICI)
>
>E: maris at nici.ru.nl
>
>MSc Cognitive Neuroscience :www.ru.nl/master/cns/
>
>
--
Max Planck Institute for Human Cognitive and Brain Sciences
Dipl.-Phys. Maren Grigutsch mailto:grigu at cbs.mpg.de
Stephanstr.1a http://www.cbs.mpg.de
04103 Leipzig, Germany phone/fax: +49 341 9940-136/-260
----------------------------------
The aim of this list is to facilitate the discussion between users of the FieldTrip toolbox, to share experiences and to discuss new ideas for MEG and EEG analysis. See also http://listserv.surfnet.nl/archives/fieldtrip.html and http://www.ru.nl/fcdonders/fieldtrip.
More information about the fieldtrip
mailing list