function [efs,F,cdfs,p,eps,dfs,b,y2,sig,eta2]=repanova(d,D,fn,gg,alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% N-way repeated measures ANOVAs
% (all factors must be repeated measures)
%
% Input:
%
% d = data A matrix with rows = replications (eg subjects) and
% columns = conditions
% D = factors A vector with as many entries as factors, each entry being
% the number of levels for that factor
%
% Data matrix d must have as many columns (conditions) as
% the product of the elements of the factor matrix D
%
% First factor rotates slowest; last factor fastest
%
% Eg, in a D=[2 3] design: factor A with 2 levels; factor B with 3:
% data matrix d must be organised:
%
% A1B1 A1B2 A1B3 A2B1 A2B2 A2B3
% rep1
% rep2
% ...
%
% fn = factor names (string array) [optional]
% gg = Greenhouse-Geiser correction flag [default=1]
% alpha = nominal alpha value [default=.05]
%
% Output:
%
% efs = effect, eg [1 2] = interaction between factor 1 and factor 2
% F = F value
% cdfs = corrected df's (using Greenhouse-Geisser)
% p = p-value
% eps = epsilon
% dfs = original dfs
% b = betas
% eta2 = partial eta-squared (measure of effect size)
%
% Rik.Henson@mrc-cbu.cam.ac.uk, 17/3/03
%
% Added eta^2 R.Henson, 17/8/08
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<5
alpha=0.05;
end
if nargin<4
gg=1; % G-G correction
end
if nargin<3 % No naming of factors provided (added 26/3/03)
for f=1:length(D)
fn{f}=sprintf('%d',f);
end
end
Nf = length(D); % Number of factors
Nd = prod(D); % Number of conditions
Ne = 2^Nf - 1; % Number of effects
Nr = size(d,1); % Number of replications (eg subjects)
sig=cell(2,1);
if size(d,2) ~= Nd
error(sprintf('data has %d conditions; design only %d',size(d,2),Nd))
end
sc = cell(Nf,2); % create main effect/interaction component contrasts
for f = 1 : Nf
sc{f,1} = ones(D(f),1);
sc{f,2} = detrend(eye(D(f)),0);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
sy = cell(Nf,2); % create main effect/interaction component contrasts
for f = 1 : Nf
sy{f,1} = ones(D(f),1)/D(f);
sy{f,2} = eye(D(f));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tss = sum(detrend(d(:),0).^2); % total sum of squares (only needed for eta2 below)
for e = 1 : Ne
cw = num2binvec(e,Nf)+1;
c = sc{1,cw(Nf)}; % create full contrasts
for f = 2 : Nf
c = kron(c,sc{f,cw(Nf-f+1)});
end
y = d * c; % project data to contrast sub-space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
cy = sy{1,cw(Nf)}; % create full contrasts
for f = 2 : Nf
cy = kron(cy,sy{f,cw(Nf-f+1)});
end
y2{e} = mean(d * cy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nc = size(y,2);
df1 = rank(c);
df2 = df1*(Nr-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% long, SPM way (GG doesn't work):
%
% y = y(:);
% X = kron(eye(nc),ones(Nr,1));
% b{e} = pinv(X)*y;
% Y = X*b{e};
% R = eye(nc*Nr)- X*pinv(X);
% r = y - Y;
%% V = r*r';
% V = y*y';
% eps(e) = trace(R*V)^2 / (df1 * trace((R*V)*(R*V)));
%
%% ss = Y'*y;
%% mse = (y'*y - ss)/df2;
%% mss = ss/df1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computationally simpler way
b{e} = mean(y); % betas
ss = sum(y*b{e}');
mse = (sum(diag(y'*y)) - ss)/df2;
mss = ss/df1;
if gg
V = cov(y); % sample covariance
eps(e) = trace(V)^2 / (df1*trace(V'*V));% Greenhouse-Geisser
else
eps(e) = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta2(e,1) = ss/tss; % Approx measure of effect size (total eta^2)
eta2(e,2) = ss/(ss+mse*df2); % Another approx measure of effect size (partial eta^2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
efs{e} = Nf+1-find(cw==2); % codes which effect
F(e) = mss/mse;
dfs(e,:) = [df1 df2];
cdfs(e,:) = eps(e)*dfs(e,:);
p(e) = 1-spm_Fcdf(F(e),cdfs(e,:));
if p(e) < alpha; sig{1}=[sig{1} e]; end
if p(e) < alpha/Ne; sig{2}=[sig{2} e]; end
en=fn{efs{e}(1)}; % Naming of factors modified (added 21/6/06)
for f = 2:length(efs{e})
en = [fn{efs{e}(f)} en];
end
disp(sprintf('Effect %02d: %-18s F(%3.2f,%3.2f)=%4.3f,\tp=%4.3f (eta2=%4.3f)',...
e,en,cdfs(e,1),cdfs(e,2),F(e),p(e),eta2(e,1)))
% disp(sprintf('Effect %02d: %-18s F(%3.2f,%3.2f)=%4.3f,\tp=%4.3f',...
% e,en,dfs(e,1),dfs(e,2),F(e),1-spm_Fcdf(F(e),dfs(e,:))))
end
disp(sprintf('\n'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
%for e = 1 : Ne
% en=fn{efs{e}(1)};
% for f = 2:length(efs{e})
% en = [en fn{efs{e}(f)}];
% end
% disp(sprintf('Effect %02d: %-18s \n',e,en))
% disp(y2{e})
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-function to code all main effects/interactions
function b = num2binvec(d,p)
if nargin<2
p = 0; % p is left-padding with zeros option
end
d=abs(round(d));
if(d==0)
b = 0;
else
b=[];
while d>0
b=[rem(d,2) b];
d=floor(d/2);
end
end
b=[zeros(1,p-length(b)) b];