
%clear
randn('state',0);		% Reset random seed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create synthetic data for 4 conditions (repeated measures)

N = 12;				% 12 subjects
P = 4;				% 4 conditions
J = N*P;			% Total number of datapoints

sigma = 1.5;			% Noise

b_c = [1 1 2 4]';		% Condition means (betas, "true" effects)
b_s = round(randn(12,1)*2);	% Subject means

y = kron(b_c,ones(N,1)) + kron(ones(P,1),b_s) + sigma*randn(J,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ONE-WAY 1x4 BETWEEN-SUBJECTS ANOVA, M=1 factor, P=4 conditions
disp('------------------------------------------------');
disp('1-WAY BETWEEN-SUBJECTS ANOVA');
disp(' ');

K = [4];
P = prod(K);

X = [kron(eye(K),ones(N,1)) ones(J,1)];

B = pinv(X)*y;			% Beta_hats (estimated effects)

figure
imagesc(X)
colormap([.5 .5 .5; 1 1 1])
%print -dtiff X_1x4_between.tif

% figure
% plot(y)
% hold on
% plot(X*B,'r')


% Two models approach
% -------------------
disp('Two models approach');
X_1 = X; B_1 = B;			% Full model
X_0 = ones(J,1); B_0 = pinv(X_0)*y;	% Reduced model (grand mean only)

r   = y - X_1*B_1;			% Residuals for full model
R_1 = r'*r;				% Residuals sum of squares for full

r   = y - X_0*B_0;			% Residuals for reduced model
R_0 = r'*r;				% Residuals sum of squares for reduced

df  = [rank(X_1)-rank(X_0) J-rank(X_1)];	
F   = ((R_0 - R_1)/df(1)) / (R_1/df(2));
p   = 1-spm_Fcdf(F,df);

sprintf('F(%d,%d)=%3.2f, p=%4.3f',K-1,J-K,F,p)
disp(' ');

% Some equivalences in Classical Sum of Squares (SS) approach
% -----------------------------------------------------------
disp('Sums of Squares approach');
SSm = y'*y;				% SS total
SSe = R_1;				% SS error

B_0 = mean(y);				% Grand mean (intercept)
SSi = J*(B_0^2);			% SS due to intercept

SSt = SSm - SSe - SSi;			% SS due to treatments
% y = y-mean(y);			
% B = pinv(X)*y;
% SSt = N*(B(1:K)'*B(1:K));		% (SS due to treatments)

df = [K-1 J-K];
MSt = SSt / df(1);			% Mean square treatments
MSe = SSe / df(2);			% Mean square error

F   = MSt/MSe;
p   = 1-spm_Fcdf(F,df);

sprintf('F(%d,%d)=%3.2f, p=%4.3f',K-1,J-K,F,p)

disp(' ');

% Classical approach using reduced F-tests
% ----------------------------------------

disp('Classical approach using F tests');
r = y - X*B;			% Residuals
RSS = r'*r;			% Residual sum of squares for "full" model

X_r = ones(J,1);		% Reduced model

B_r = pinv(X_r)*y;
r_r = y - X_r*B_r;

RSS_r = r_r'*r_r;		% Residual sum of squares for "reduced" model

df = [K-1 J-K];
F  = ((RSS_r - RSS)/df(1)) / (RSS/df(2));
p  = 1-spm_Fcdf(F,df);

sprintf('F(%d,%d)=%3.2f, p=%4.3f',df(1),df(2),F,p)
disp(' ');

% Equivalent generalisation (a la SPM, using residual forming matrices R)
% -----------------------------------------------------------------------
disp('SPM approach');

c = [eye(K) zeros(K,1)]';	% Effects of interest contrast
c = [diff(eye(4)) zeros(3,1)]';	% Equivalent contrast

c_0 = eye(size(X,2)) - c*pinv(c);
X_0 = X*c_0;

R   = eye(size(X,1)) - X*pinv(X);
R_0 = eye(size(X,1)) - X_0*pinv(X_0);

M = R_0 - R;

df = [rank(X)-rank(X_0) size(X,1)-rank(X)];
F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
p  = 1-spm_Fcdf(F,df);

sprintf('F(%d,%d)=%3.2f, p=%4.3f',df(1),df(2),F,p)
disp(' ');


disp(' ');
disp(' ');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ONE-WAY 1x4 WITHIN-SUBJECTS ANOVA, M=1 factor, K=4 conditions
disp('------------------------------------------------');
disp('1-WAY WITHIN-SUBJECTS ANOVA');


X = [kron(eye(K),ones(N,1)) kron(ones(K,1),eye(N))];

B = pinv(X)*y;			% Beta_hats (estimated effects)

figure
imagesc(X)
colormap([.5 .5 .5; 1 1 1])
%print -dtiff X_1x4_within.tif

figure
plot(y)
hold on
plot(X*B,'r')

c = [eye(K) zeros(K,N)]';		% Effects of interest contrast

c_0 = eye(size(X,2)) - c*pinv(c);
X_0 = X*c_0;

R   = eye(size(X,1)) - X*pinv(X);
R_0 = eye(size(X,1)) - X_0*pinv(X_0);

M = R_0 - R;

df = [rank(X)-rank(X_0) size(X,1)-rank(X)];
F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
p  = 1-spm_Fcdf(F,df);

sprintf('F(%d,%d)=%3.2f, p=%4.3f (assuming sphericity)',df(1),df(2),F,p)

% G-G - get trace(RV) working....

Y = reshape(y,N,P);
V = cov(Y);

eps = trace(V)^2 / (df(1)*trace(V'*V));	
%eps = trace(R*V)^2 / trace((R*V)*(R*V));

df = [rank(X)-rank(X_0) size(X,1)-rank(X)]*eps;
F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
p  = 1-spm_Fcdf(F,df);

sprintf('F(%3.2f,%3.2f)=%3.2f, p=%4.3f (episilon=%3.2f)',df(1),df(2),F,p,eps)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TWO-WAY 2x2 WITHIN-SUBJECTS ANOVA, M=2 factors, K_1=K_2=2 => 4 conditions
% POOLED ERROR TERM! (BUT NONSPHERICAL?)

K = [2 2];			% Levels of each factor
M = length(K);			% Number of factors
P = prod(K);			% Number of conditions

for m = 1:M
	com{m} = ones(1,K(m));
	dif{m} = diff(eye(K(m)));
end

% ANOVA effect contrasts spelled out below
% (For generalisation to K_1 x K_2 x ... K_M ANOVAs, see repanova.m)

ci = [	kron(dif{1},com{2});		% Main effect of factor 1
 	kron(com{2},dif{1});		% Main effect of factor 2
	kron(dif{1},dif{2});		% Interaction of factor 1 and factor 2
	kron(com{1},com{2})		% Constant
     ];

X = [kron(eye(P),ones(N,1)) kron(ones(P,1),eye(N))];

figure
imagesc(X)
colormap([.5 .5 .5; 1 1 1])
%print -dtiff X_2x2_within.tif

B = pinv(X)*y;			% Beta_hats (estimated effects)

% nonsphericity
Y = reshape(y,N,P);
V = cov(Y);
df = P-1;
eps = trace(V)^2 / (df(1)*trace(V'*V));	

for d = 1:(P-1)			% Constant not interesting (first 3 effects)

	c = [ci(:,d); zeros(N,1)];

	c_0 = eye(size(X,2)) - c*pinv(c);
	X_0 = X*c_0;

	R   = eye(size(X,1)) - X*pinv(X);
	R_0 = eye(size(X,1)) - X_0*pinv(X_0);

	M = R_0 - R;

	df = [rank(X)-rank(X_0) size(X,1)-rank(X)];

%	F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
%	p  = 1-spm_Fcdf(F,df);
%	sprintf('%s: F(%d,%d)=%3.2f, p=%4.3f',mat2str(ci(d,:)),df(1),df(2),F,p)

	df = df * eps;
	F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
	p  = 1-spm_Fcdf(F,df);

	sprintf('F(%3.2f,%3.2f)=%3.2f, p=%4.3f (episilon=%3.2f)',df(1),df(2),F,p,eps)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TWO-WAY 2x2 WITHIN-SUBJECTS ANOVA, M=2 factors, K_1=K_2=2 => 4 conditions
% EQUIVALENT ROTATED DESIGN MATRIX

cX = [ci zeros(P,N); zeros(N,P) eye(N)]';

% Rotation of design matrix into standard ANOVA contrasts

X = [kron(eye(P),ones(N,1)) kron(ones(P,1),eye(N))];
X = X*cX;

figure
imagesc(X)
colormap([0 0 0; .5 .5 .5; 1 1 1])
%print -dtiff X_2x2_within.tif

B = pinv(X)*y;			% Beta_hats (estimated effects)

ci = eye(P);			% Though only columns 1-3 in X are of interest

for d = 1:P-1

	c = [ci(:,d); zeros(N,1)];

	c_0 = eye(size(X,2)) - c*pinv(c);
	X_0 = X*c_0;

	R   = eye(size(X,1)) - X*pinv(X);
	R_0 = eye(size(X,1)) - X_0*pinv(X_0);

	M = R_0 - R;

	df = [rank(X)-rank(X_0) size(X,1)-rank(X)];
	F  = ((B'*X'*M*X*B)/df(1)) / ((y'*R*y)/df(2));
	p  = 1-spm_Fcdf(F,df);

	sprintf('%s: F(%d,%d)=%3.2f, p=%4.3f',mat2str(ci(d,:)),df(1),df(2),F,p)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TWO-WAY 2x2 WITHIN-SUBJECTS ANOVA, M=2 factors, K_1=K_2=2 => 4 conditions
% FULL DESIGN MATRIX

K = [2 2];			% Levels of each factor
M = length(K);			% Number of factors
P = prod(K);			% Number of conditions

for m = 1:M
	com{m} = ones(1,K(m));
	dif{m} = diff(eye(K(m)));
end

% ANOVA effect contrasts spelled out below
% (For generalisation to K_1 x K_2 x ... K_M ANOVAs, see repanova.m)

ci = [	kron(dif{1},com{2});		% Main effect of factor 1
 	kron(com{2},dif{1});		% Main effect of factor 2
	kron(dif{1},dif{2});		% Interaction of factor 1 and factor 2
	kron(com{1},com{2})		% Constant
     ];

X = [kron(ci,ones(N,1)) kron(ci,eye(N))];

figure
imagesc(X)
colormap([0 0 0; .5 .5 .5; 1 1 1])
%print -dtiff X_2x2_within.tif

B = pinv(X)*y;			% Beta_hats (estimated effects)

ci = eye(P);			% Though only columns 1-3 in X are of interest

for d = 1:(P-1)			% Constant not interesting (first 3 effects)

%	c   = [ci(:,d) zeros(P,N); zeros(J,1) kron(ci(:,d),eye(N))];
	c   = [ci(d,:) zeros(1,N*P); zeros(N,P) kron(ci(d,:),eye(N))]';
	c_0 = eye(size(X,2)) - c*pinv(c);
	X_0 = X*c_0;

%	c   = [zeros(P,N); kron(ci(:,d),eye(N))];
	c   = [zeros(N,P)  kron(ci(d,:),eye(N))]';
	c_1 = eye(size(X,2)) - c*pinv(c);
	X_1 = X*c_1;

	R_0 = eye(size(X,1)) - X_0*pinv(X_0);
	R_1 = eye(size(X,1)) - X_1*pinv(X_1);

	M = R_0 - R_1;

	df = [rank(X_1)-rank(X_0) size(X_1,1)-rank(X_1)];

	B = pinv(X_1)*y;

	F  = ((B'*X_1'*M*X_1*B)/df(1)) / ((y'*R_1*y)/df(2));
	p  = 1-spm_Fcdf(F,df);

	sprintf('F(%d,%d)=%3.2f, p=%4.3f',df(1),df(2),F,p)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TWO-WAY 2x2 WITHIN-SUBJECTS ANOVA, M=2 factors, K_1=K_2=2 => 4 conditions
% CONTRAST SPECIFIC ERROR TERMS, REPRODUCED IN TWO-STAGE MODEL

K = [2 2];			% Levels of each factor
M = length(K);			% Number of factors
P = prod(K);			% Number of conditions

for m = 1:M
	com{m} = ones(1,K(m));
	dif{m} = diff(eye(K(m)));
end

ci = [	kron(dif{1},com{2});		% Main effect of factor 1
 	kron(com{2},dif{1});		% Main effect of factor 2
	kron(dif{1},dif{2});		% Interaction of factor 1 and factor 2
	kron(com{1},com{2})		% Constant
     ];

X_1 = [kron(ci,ones(N,1)) kron(ci,eye(N))];

for d = 1:3

	c_1 = kron(ci(:,d),eye(12));

	y_2 = pinv(c_1)*y;

%	X_2 = pinv(c_1)*X_1		% Equivalent reduction of full model

	X_2 = ones(N,1);
	B_2 = pinv(X_2)*y_2;

	c_2 = [1];
	c_0 = eye(size(X_2,2)) - c_2*pinv(c_2);
	X_0 = X_2*c_0;

	R_2  = eye(size(X_2,1)) - X_2*pinv(X_2);
	R_0 = eye(size(X_2,1)) - X_0*pinv(X_0);

	M = R_0 - R_2;

	df = [rank(X_2)-rank(X_0) size(X_2,1)-rank(X_2)];
	F  = ((B_2'*X_2'*M*X_2*B_2)/df(1)) / ((y_2'*R_2*y_2)/df(2));
	p  = 1-spm_Fcdf(F,df);

	sprintf('%s: F(%d,%d)=%3.2f, p=%4.3f',mat2str(ci(d,:)),df(1),df(2),F,p)
end

% 2nd-level model for 3x3 interaction

X = kron(eye(4),ones(12,1));
figure
imagesc(X)
colormap([.5 .5 .5; 1 1 1])


