% ANOVA2RM_CELL - compute F-values in cell array using repeated measure
% ANOVA.
%
% Usage:
% >> [FC FR FI dfc dfr dfi] = anova2rm_cell( data );
%
% Inputs:
% data = data consisting of PAIRED arrays to be compared. The last
% dimension of the data array is used to compute ANOVA.
% Outputs:
% FC - F-value for columns.
% FR - F-value for rows.
% FI - F-value for interaction.
% dfc - degree of freedom for columns.
% dfr - degree of freedom for rows.
% dfi - degree of freedom for interaction.
%
% Note: this function is inspired from rm_anova available at
% http://www.mathworks.se/matlabcentral/fileexchange/6874-two-way-rep
% eated-measures-anova
% It allows for fast computation of about 20 thousands ANOVA per
% second. It is different from anova2_cell which mimics the ANOVA
% function from the Matlab statistical toolbox. This function
% computes true repeated measure ANOVA.
%
% Example:
% a = { rand(1,10) rand(1,10) rand(1,10); rand(1,10) rand(1,10) rand(1,10) }
% [FC FR FI dfc dfr dfi] = anova2rm_cell(a)
% signifC = 1-fcdf(FC, dfc(1), dfc(2))
% signifR = 1-fcdf(FR, dfr(1), dfr(2))
% signifI = 1-fcdf(FI, dfi(1), dfi(2))
%
% % for comparison
% z = zeros(10,1); o = ones(10,1); t = ones(10,1)*2;
% rm_anova2( [ a{1,1}';a{1,2}';a{1,3}';a{2,1}';a{2,2}';a{2,3}' ], ...
% repmat([1:10]', [6 1]), [o;o;o;z;z;z], [z;o;t;z;o;t], {'a','b'})
%
% c = { rand(200,400,10) rand(200,400,10); ...
% rand(200,400,10) rand(200,400,10)};
% [FC FR FI dfc dfr dfi] = anova2rm_cell(c) % computes 200x400 ANOVAs
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2010
% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [fA fB fAB dfApair dfBpair dfABpair] = anova2rm_cell(data)
% compute all means and all std
% -----------------------------
a = size(data,1);
b = size(data,2);
nd = myndims( data{1} );
n = size( data{1} ,nd);
% only for paired stats
% ---------------------
if nd == 1
AB = zeros(a,b,'single');
AS = zeros(a,n,'single');
BS = zeros(b,n,'single');
sq = single(0);
for ind1 = 1:a
for ind2 = 1:b
AB(ind1,ind2) = sum(data{ind1,ind2});
AS(ind1,:) = AS(ind1,:) + data{ind1,ind2}';
BS(ind2,:) = BS(ind2,:) + data{ind1,ind2}';
sq = sq + sum(data{ind1,ind2}.^2);
end
end
dimA = 2;
dimB = 1;
elseif nd == 2
AB = zeros(size(data{1},1),a,b,'single');
AS = zeros(size(data{1},1),a,n,'single');
BS = zeros(size(data{1},1),b,n,'single');
sq = zeros(size(data{1},1),1,'single');
for ind1 = 1:a
for ind2 = 1:b
AB(:,ind1,ind2) = sum(data{ind1,ind2},nd);
AS(:,ind1,:) = AS(:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),1,n);
BS(:,ind2,:) = BS(:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),1,n);
sq = sq + sum(data{ind1,ind2}.^2,nd);
end
end
dimA = 3;
dimB = 2;
elseif nd == 3
AB = zeros(size(data{1},1),size(data{1},2),a,b,'single');
AS = zeros(size(data{1},1),size(data{1},2),a,n,'single');
BS = zeros(size(data{1},1),size(data{1},2),b,n,'single');
sq = zeros(size(data{1},1),size(data{1},2),'single');
for ind1 = 1:a
for ind2 = 1:b
AB(:,:,ind1,ind2) = sum(data{ind1,ind2},nd);
AS(:,:,ind1,:) = AS(:,:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),1,n);
BS(:,:,ind2,:) = BS(:,:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),1,n);
sq = sq + sum(data{ind1,ind2}.^2,nd);
end
end
dimA = 4;
dimB = 3;
elseif nd == 4
AB = zeros(size(data{1},1),size(data{1},2),size(data{1},3),a,b,'single');
AS = zeros(size(data{1},1),size(data{1},2),size(data{1},3),a,n,'single');
BS = zeros(size(data{1},1),size(data{1},2),size(data{1},3),b,n,'single');
sq = zeros(size(data{1},1),size(data{1},2),size(data{1},3),'single');
for ind1 = 1:a
for ind2 = 1:b
AB(:,:,:,ind1,ind2) = sum(data{ind1,ind2},nd);
AS(:,:,:,ind1,:) = AS(:,:,:,ind1,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),size(data{1},3),1,n);
BS(:,:,:,ind2,:) = BS(:,:,:,ind2,:) + reshape(data{ind1,ind2},size(data{1},1),size(data{1},2),size(data{1},3),1,n);
sq = sq + sum(data{ind1,ind2}.^2,nd);
end
end
dimA = 5;
dimB = 4;
end
A = sum(AB,dimA); % sum across columns, so result is ax1 column vector
B = sum(AB,dimB); % sum across rows, so result is 1xb row vector
S = sum(AS,dimB); % sum across columns, so result is 1xs row vector
T = sum(sum(A,dimB),dimA); % could sum either A or B or S, choice is arbitrary
% degrees of freedom
dfA = a-1;
dfB = b-1;
dfAB = (a-1)*(b-1);
dfS = n-1;
dfAS = (a-1)*(n-1);
dfBS = (b-1)*(n-1);
dfABS = (a-1)*(b-1)*(n-1);
% bracket terms (expected value)
expA = sum(A.^2,dimB)./(b*n);
expB = sum(B.^2,dimA)./(a*n);
expAB = sum(sum(AB.^2,dimA),dimB)./n;
expS = sum(S.^2,dimA)./(a*b);
expAS = sum(sum(AS.^2,dimB),dimA)./b;
expBS = sum(sum(BS.^2,dimB),dimA)./a;
expY = sq; %sum(Y.^2);
expT = T.^2 / (a*b*n);
% sums of squares
ssA = expA - expT;
ssB = expB - expT;
ssAB = expAB - expA - expB + expT;
ssS = expS - expT;
ssAS = expAS - expA - expS + expT;
ssBS = expBS - expB - expS + expT;
ssABS = expY - expAB - expAS - expBS + expA + expB + expS - expT;
ssTot = expY - expT;
% mean squares
msA = ssA / dfA;
msB = ssB / dfB;
msAB = ssAB / dfAB;
msS = ssS / dfS;
msAS = ssAS / dfAS;
msBS = ssBS / dfBS;
msABS = ssABS / dfABS;
% f statistic
fA = msA ./ msAS;
fB = msB ./ msBS;
fAB = msAB ./ msABS;
dfApair = [dfA dfAS];
dfBpair = [dfB dfBS];
dfABpair = [dfAB dfABS];
function val = myndims(a)
if ndims(a) > 2
val = ndims(a);
else
if size(a,1) == 1,
val = 2;
elseif size(a,2) == 1,
val = 1;
else
val = 2;
end
end
Datasets
Models
