Back to the main page.
Bug 1556 - bug in paired ANOVA
| Status | CLOSED FIXED |
| Reported | 2012-06-26 09:03:00 +0200 |
| Modified | 2019-08-10 12:28:38 +0200 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P3 normal |
| Assigned to: | Diego Lozano Soldevilla |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Robert Oostenveld - 2012-06-26 09:03:12 +0200
Created attachment 285 screen shot from Arno On 26 Jun 2012, at 6:15, Arnaud Delorme wrote: Hi Robert, I am trying to interface Fieldtrip statistics once again (to be able to use your cluster function). In the process, I was performing some serious testings to check the consistency between EEGLAB and Fieldtrip statistical results. One problem arose with the paired ANOVA measure of Fieldtrip I have laboriously written a pure Fieldtrip script below to show you the problem. One first problem is the degree of freedom for the denominator which should be 18 and not 8 (and other problems may be coming from these as well). I have found an interactive web page where I could enter the same data for testing and I have indicated the results as well. The results of the web page are consistent with EEGLAB statistic functions but not Fieldtrip. I hope I am not passing next to something trivial. Thanks, Arno ps: all other statistical function of Fieldtrip I have tested (unpaired Anova, paired t-test and unpaired t-test) worked ok so far. By the way, these results were obtained using fieldtrip-20120525 from last month. -- cfg.verbose = 'off'; cfg.method = 'analytic'; cfg.feedback = 'no'; cfg.alpha = 5.0000e-02; cfg.tail = 1; cfg.statistic = 'depsamplesF'; cfg.ivar = 1; cfg.uvar = 2; cfg.design = [ 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3; 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 ]; data1.powspctrm = [1 0 3 8 8 10 1 4 6 1]'; data1.dimord = 'rpt_chan_freq_time'; data1.label = {'cz'}; data1.freq = 1; data1.time = 1; data2.powspctrm = [1 8 6 1 6 8 7 4 9 0]'; data2.dimord = 'rpt_chan_freq_time'; data2.label = {'cz'}; data2.freq = 1; data2.time = 1; data3.powspctrm = [10 5 6 7 10 5 11 6 10 3]'; data3.dimord = 'rpt_chan_freq_time'; data3.label = {'cz'}; data3.freq = 1; data3.time = 1; stats = ft_freqstatistics(cfg, data1, data2, data3) stat: 6.1009e+00 dfnum: 2 dfdenom: 8 critval: 1.0033e+01 prob: 1.2617e-01 mask: 0 dimord: 'chan' label: {'cz'} freq: 1 time: 1 cfg: [1x1 struct] http://faculty.vassar.edu/lowry/corr3.html <Screen Shot 2012-06-25 at 8.55.35 PM.png>
Robert Oostenveld - 2012-07-04 16:09:27 +0200
we are working on it (including Lilla, CC) I have just copied the test code (thanks Arno!) to a test script and added some more lines. mbp> svn commit test_bug1556.m Adding test_bug1556.m Transmitting file data . Committed revision 6234. more time needs to be spent on this
Lilla Magyari - 2012-07-04 18:46:46 +0200
(In reply to comment #1) Hi, I compared the anovan function of matlab and the statfun_desamplesF in FieldTrip. To sum up (but details are below), it seems to me that the statfun_desamplesF doesn't give the same output when the experimental condition has more than two factor levels. (I will attach the script, called bug_desamplesF.) First, I compared the output of the anova function on a sample dataset of Matlab (called mileage) where I followed the steps described in the Matlab documentation (ANOVA with Random Effects). (For the computation being comparable to FT, first I had to delete some rows in mileage, because statfun_desamplesF can handle only 1 data per random factor per condition.) In this example data set, the fixed effect (called carmod) has 2-levels. The anovan gave the same F values for the fixed effect than the statfun_desamplesF. Then, I computed the same on the example data from Arno, and I got different F value with anovan and statfun_depsamplesF. When I deleted one of the conditions (there were originally 3), the F values were the same again. Lilla
Lilla Magyari - 2012-07-04 18:49:55 +0200
Created attachment 290 comparison with anovan
Irina - 2012-08-08 16:16:46 +0200
(In reply to comment #3) Hi , I checked it with my data and just want to confirm what Lilla already said: ft_depsamplesF output is different from the the anovan output if independent variable (cfg.ivar) has more than two levels.
Diego Lozano Soldevilla - 2014-01-09 17:49:15 +0100
(In reply to Irina from comment #4) I made progress with this bug. Currently there's still a discrepancy between ft_statfun_depsamplesF and what SPSS or Arnaud's site. I've not been able to discover the F dependent sample model that ft_statfun_depsamplesF is using but I found 2 common ways to compute the F stat for repeated measures ANOVA to test a main effect. Both ways gave me same results as SPSS and http://faculty.vassar.edu/lowry/corr3.html You can replicate my results as follows: %% Method 1: GLM way. % Source:visit http://www.sbirc.ed.ac.uk/cyril/glm/GLM_lectures.html#10 for an excellent tutorial on GLM and here to find formulas http://www.fil.ion.ucl.ac.uk/~wpenny/publications/rik_anova.pdf u1 = [1 0 3 8 8 10 1 4 6 1]'; u2 = [1 8 6 1 6 8 7 4 9 0]'; u3 = [10 5 6 7 10 5 11 6 10 3]'; Y = [u1; u2; u3]; % this time this not 3 gp but 10 subjects with 3 measures % create the design matrix for the different factors nb_subjects =10; nb_conditions =3; Subjects = repmat(eye(nb_subjects),nb_conditions,1); % error x = kron(eye(nb_conditions),ones(nb_subjects,1)); % effect X = [x Subjects]; % no more ones for the grand mean but a subject specific mean figure; imagesc(X); colormap('gray'); title('Repearted measure design','Fontsize',14) % Compute as usual df = nb_conditions -1; dfe = size(Y,1) - nb_subjects - df; P = X*pinv(X'*X)*X'; % our projection matrix R = eye(size(Y,1)) - P; % projection on error space SSe = diag(Y'*R*Y); % Y projected onto the error Betas = pinv(x)*Y; % compute without cst/subjects yhat = x*Betas; % yhat computed based on the treatment with subject effect - we use little x SS = norm(yhat-mean(yhat)).^2; F_values = (SS/df) ./ (SSe/dfe); p_values = 1 - fcdf(F_values, df, dfe); %% Method 2: Grandmother counting way. You can find the procedure here: https://statistics.laerd.com/statistical-guides/repeated-measures-anova-statistical-guide-2.php u1 = [1 0 3 8 8 10 1 4 6 1]'; u2 = [1 8 6 1 6 8 7 4 9 0]'; u3 = [10 5 6 7 10 5 11 6 10 3]'; Y = [u1; u2; u3]; % this time this not 3 gp but 10 subjects with 3 measures nb_subjects =10; nb_conditions =3; Y = reshape(Y,nb_subjects,nb_conditions); Ysub = mean(Y,2); Yfac = mean(Y,1); % computing sums of squeares SStot = sum(sum((Y-mean(Ysub)).^2)); SSsub = nb_conditions * sum(sum((Ysub-mean(Ysub)).^2)); SSfac = nb_subjects * sum(sum((Yfac-mean(Ysub)).^2)); SSerr = SStot - SSsub - SSfac; % mean sum of squares for factor levels df = nb_conditions - 1; dfe = size(Y(:),1) - nb_subjects - df; MSfac = SSfac/df; MSerr = SSerr/dfe; Fstat = MSfac/MSerr; pval = 1-fcdf(Fstat,df,dfe); Note that I'm using same toy data example as Arnaud getting the same results. Now the question: what is the F model that ft_statfun_depsamplesF is testing? In the meantime, I can write down a ft_statfun with one of the above methods Diego
Diego Lozano Soldevilla - 2014-01-14 11:13:33 +0100
Dear Eric, I'm continuing debuging the ft_statfun_depsamplesF function and I made some progress. The function ft_statfun_depsamplesF is providing wrong degrees of freedom for the error term and somewhat different F statistic as other software packages (SPSS and here http://faculty.vassar.edu/lowry/corr3.html). As both procedures did not provide the open source code I found 2 ways to compute the ANOVA main effect for repeated measures (see my previous comment to replicate my conclusions). Both procedures (based on GLM and the sums of squares) gave me same identical results between them and also identical to SPSS and the website. I guess that there are multiple ANOVA models to test main effects but I coudn't discover what's the ft_statfun_depsamplesF way and I'd like to ask you if you could provide us a reference on were the original implementation was based and if you can confirm my results. In case both ANOVA models are valid, I propose to implement both in different functions. best, Diego
Diego Lozano Soldevilla - 2014-01-17 14:22:37 +0100
(In reply to Diego Lozano Soldevilla from comment #6) Eric Maris pointed correctly that the model ft_statfun_depsamplesF function is testing is a MANOVA and NOT repeated measures ANOVA. Here critical point that explains the origin of the bug: it has been a documentation issue and not a wrong calculation Plans: 1) Improve ft_statfun_depsamplesF to state clearly that MANOVA is the model one's testing 2) Create a ft_statfun_depsamplesFmixedmodel to provide users with a repeated measures ANOVA testing a mixed model Fstatistic, as SPSS does. Thanks a lot for your help Eric! Diego
Diego Lozano Soldevilla - 2014-01-17 14:48:08 +0100
(In reply to Diego Lozano Soldevilla from comment #7) Eric suggested the following function names to avoid confusions: ft_statfun_depsamplesFmultivariate -> current ft_statfun_depsamplesF (MANOVA) ft_statfun_depsamplesFunivariate -> future new function for repeated measures ANOVA Before change names, I made some checking to avoid bugs to other FT function and tutorials etc: Regading wiki contents, only documentation will be affected but no issues with tutorial/FAQs/example scripts code: http://fieldtrip.fcdonders.nl/?do=search&id=depsamplesF Searching through the entire FT code, only one function can be affected. The rest is just documentation issues bash-4.1$ find -type f -name "*.m" -print0 | xargs -0 grep -l ft_statfun_depsamplesF ./statfun/ft_statfun_depsamplesF.m bash-4.1$ find -type f -name "*.m" -print0 | xargs -0 grep -l depsamplesF ./private/prepare_design.m --> ****updated needed**** ./test/test_bug1556.m ./statfun/ft_statfun_indepsamplesF.m ./statfun/ft_statfun_depsamplesF.m ./ft_statistics_montecarlo.m ./ft_statistics_analytic.m
Diego Lozano Soldevilla - 2014-05-21 10:28:19 +0200
New univariate Fstat for repeated measures added with test script Adding statfun/ft_statfun_depsamplesFunivariate.m Sending test/test_bug1556.m Transmitting file data .. Committed revision 9562.