bug in statistics_montecarlo?

[Maren Grigutsch]

Dear Robert, Eric and other fellows of the Fieldtrip team,

I'm just trying out your excellent software, in specific your extremely
useful and well designed statistics functions. Many thanks for this
great and generous gift to the neuroscience community!

I would like to point you to an inconsistency I seem to have found  in
the implementation of the statistics_montecarlo function. When I
compared the probabilities computed with statistics_analytic and
statistics_montecarlo for the depsamplesT test, the probabilities
determined via the randomization test turned out to be too low by a
factor of about 1/2, see code and Matlab output below.

-------------------------

N = 50;
M = 0.25;
alpha = 0.05;
tail = 0;     % 'both'

% construct a normally distributed difference vector
% with unit variance and mean M
d = randn(1,N);
d = d ./ sqrt(var(d,1,2));
d = d -mean(d) + M;

% compute a  paired t-test
[H,P,CI,stat] = ttest(d,0,alpha,tail)

H =

     0


P =

    0.0864


CI =

   -0.0371    0.5371


stat =

    tstat: 1.7500
       df: 49
       sd: 1.0102

% generate two sample vectors with their difference given by d
x =randn(1,N);
y = x + d;

% generate the design matrix for testing with fieldtrip
design = [
  1:N              1:N                                % subject number
  repmat(1,[1 N])  repmat(2,[1 N])     % condition number
];

cfg = [];
cfg.statistic        = 'depsamplesT';
cfg.method           = 'analytic';
cfg.correctm         = 'no';
cfg.alpha            = 0.05;
cfg.tail             = 0;
cfg.uvar             = 1;
cfg.ivar             = 2;

 [stat_analytic, cfg] = statistics_analytic(cfg,[x y],design)

stat_analytic =

       stat: -1.7500
         df: 49
    critval: [-2.0096 2.0096]
       prob: 0.0864
       mask: 0

cfg.numrandomization = 5000;
[stat_montecarlo, cfg] = statistics_montecarlo(cfg,[x y],design)

using "statfun_depsamplesT" for the single-sample statistics
constructing randomized design
total number of measurements     = 100
total number of variables        = 2
number of independent variables  = 1
number of unit variables         = 1
number of control variables      = 0
number of within-block variables = 0
repeated measurement in variable 1 over 50 levels
number of repeated measurements in each level is 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
estimated time per randomization is 0 seconds
not performing a correction for multiple comparisons

stat_montecarlo =

    prob: 0.0394
    mask: 1
    stat: -1.7500

--------------------

Maybe I'm missing something. But I think you should replace

stat.prob = min(prb_neg, prb_pos);
 by
stat.prob = 2 *min(prb_neg, prb_pos);

in line 344 of statistics_montecarlo(),  for the two-tailed test.

Maybe it would be even better to compare the absolute values of the
observed and the randomized statistics, i.e. doing something like this

for i=1:Nrand,
....
if tail == 0,
prb_pos = prb_pos + (abs(statobs)<abs(statrand));
end
...
end
if tail == 0,
stat.prob = prb_pos./Nrand;
end


What do you think?

Many thanks in advance for your comments.
Best regards,
--Maren

--
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.