function [AOVp] = AOVp(X,alpha)
%AOVP Planned (a priori) Analysis of Variance Test.
%
% Syntax: function [AOVp] = AOVp(X,alpha)
%
% Inputs:
% X - data matrix (Size of matrix must be n-by-2; data=column 1, sample=column 2).
% alpha - significance level (default = 0.05).
% Outputs:
% - Checking of orthogonality.
% (If it is not, the file WARNING you [Warning: One or more of the contrasts are
% nonorthogonal. You must to check the planned comparisions or the vector integers.
% Some of its sum of squares could be unreal. The addition don't corresponds to
% the treatments sum of squares.])
% - Comparision means table.
% - Complete Analysis of Variance table.
%
% Example: From the example of Sokal and Rohlf (1981, pp.234-238), to compare the
% homogeneity among the 5(=k) means in a planned (a priori) analysis of
% variance with a significance level = 0.001.
%
% T r e a t m e n t
% ---------------------------------------
% 1 2 3 4 5
% 0 2%G 2%F 1%G+1%F 2%S
% ---------------------------------------
% 75 57 58 58 62
% 67 58 61 59 66
% 70 60 56 58 65
% 75 59 58 61 63
% 65 62 57 57 64
% 71 60 56 56 62
% 67 60 61 58 65
% 67 57 60 57 65
% 76 59 57 57 62
% 68 61 58 59 67
% ---------------------------------------
%
% Here, before the experiment has been carried out and the results
% obtained, we were interested to test whether the control was different
% from the rest of the treatments representing addition of sugar
% on the length of pea sections. There are two main groups, one the
% control and the other the sugars group. Next we test whether the
% mixture of sugars is different from the pure sugars. Finally we
% test any difference from the pure sugars.
% On planned comparisions the restriction is that they must to be
% orthogonal: each comparision tests have an independent relationship
% among the means. For it is necessary to extract the coefficients of
% linear comparisions. For convenience, they are generally expressed as
% integers thought they may be fractional. This is done to indicate
% the contrast. Thus, each treatment mean is weighted by its coefficient.
%
% The number of comparisions to test are k-1.
%
% T r e a t m e n t
% ---------------------------------------
% 1 2 3 4 5
% ---------------------------------------
% 1 4 -1 -1 -1 -1
% Comparision 2 0 1 1 -3 1
% 3 0 1 1 0 -2
% 4 0 1 -1 0 0
% ---------------------------------------
%
% Two comparisions are orthogonal if only if the sum of the sample
% sizes times the products of their coefficients equals zero.
%
% Data matrix must be:
% X=[75 1;67 1;70 1;75 1;65 1;71 1;67 1;67 1;76 1;68 1;
% 57 2;58 2;60 2;59 2;62 2;60 2;60 2;57 2;59 2;61 2;
% 58 3;61 3;56 3;58 3;57 3;56 3;61 3;60 3;57 3;58 3;
% 58 4;59 4;58 4;61 4;57 4;56 4;58 4;57 4;57 4;59 4;
% 62 5;66 5;65 5;63 5;64 5;62 5;65 5;65 5;62 5;67 5];
%
% Calling on Matlab the function:
% AOVp(X,0.001)
%
% -Immediately it display:
% The number of treatments are: 5
% -and:
% The number of comparisions are at least: 4
% -Then it asks:
% Please, confirm the number of comparisions: 4
% -Also it asks (be careful to give the values among brackets):
% Give me the coefficient vectors to check orthogonality:
% Give me the comparision 1:[4 -1 -1 -1 -1]
% Give me the comparision 2:[0 1 1 -3 1]
% Give me the comparision 3:[0 1 1 0 -2]
% Give me the comparision 4:[0 1 -1 0 0]
%
% -------------------------------------
% Comparision Orthogonality
% -------------------------------------
% 4 3 It is orthogonal
% 4 2 It is orthogonal
% 4 1 It is orthogonal
% 3 2 It is orthogonal
% 3 1 It is orthogonal
% 2 1 It is orthogonal
% -------------------------------------
%
% -Finally it asks (be careful to give the elements among brackets and
% with the transpose symbol.):
% Give me the comparision pair 1:[T1']
% Give me the comparision pair 2:[T2' T3' T4' T5']
%
% Give me the comparision pair 1:[T2' T3' T5']
% Give me the comparision pair 2:[T4']
%
% Give me the comparision pair 1:[T2' T3']
% Give me the comparision pair 2:[T5']
%
% Give me the comparision pair 1:[T2']
% Give me the comparision pair 2:[T3']
%
% Answer is:
% ---------------------------------------
% M e a n s
% Comparision 1 2
% ---------------------------------------
% 1 70.100 59.900
% 2 60.533 58.000
% 3 58.750 64.100
% 4 59.300 58.200
% ---------------------------------------
%
% Planned Analysis of Variance Table.
% -------------------------------------------------------------------------------
% SOV SS df MS F P Decision
% -------------------------------------------------------------------------------
% Treatments 1077.320 4
% 1 832.320 1 832.320 152.564 0.0000 S
% 2 48.133 1 48.133 8.823 0.0048 NS
% 3 190.817 1 190.817 34.977 0.0000 S
% 4 6.050 1 6.050 1.109 0.2979 NS
% Error 245.500 45 5.456
% Total 1322.820 49
% -------------------------------------------------------------------------------
% With a given significance level of: 0.0010
% The results are significant (S) and/or not significant (NS).
%
% Created by A. Trujillo-Ortiz and R. Hernandez-Walls
% Facultad de Ciencias Marinas
% Universidad Autonoma de Baja California
% Apdo. Postal 453
% Ensenada, Baja California
% Mexico.
% atrujo@uabc.mx
%
% Copyright (C) June 5, 2003.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). AOVp: Planned (a priori) analysis of variance
% test. A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
% loadFile.do?objectId=3578
%
% References:
%
% Sokal, R. R. and Rohlf, F. J. (1981), Biometry. 2nd. ed.
% New-York:W. H. Freeman & Co. pp. 234-238.
%
if nargin < 2,
alpha = 0.05; %(default)
end;
if (alpha <= 0 | alpha >= 1)
fprintf('Warning: significance level must be between 0 and 1\n');
return;
end;
k=max(X(:,2));
fprintf('The number of treatments are:%2i\n\n', k);
n=[];s2=[];A=[];
indice=X(:,2);
for i=1:k
Xe=find(indice==i);
eval(['T' num2str(i) '=X(Xe,1);']);
eval(['n' num2str(i) '=length(T' num2str(i) ') ;'])
eval(['s2' num2str(i) '=(std(T' num2str(i) ').^2) ;'])
eval(['x =(sum(T' num2str(i) ').^2)/(n' num2str(i) ');']);
eval(['tn= n' num2str(i) ';'])
eval(['ts2= s2' num2str(i) ';'])
n=[n;tn];s2=[s2;ts2];A=[A,x];
end;
C=(sum(X(:,1)))^2/length(X(:,1)); %correction term
SST=sum(A)-C; %treatments sum of squares
dfT=k-1; %degrees of freedom of tratments
dfE=sum(n)-k; %degrees of freedom of error term
summ=0;
for i=1:k
eval(['summ =summ + (n' num2str(i) '-1)*s2' num2str(i) ';']);
end;
SSE=summ; %error sum of squares
MSE = SSE/dfE; %error mean square
co = k-1;
fprintf('The number of comparisions are at least:%2i\n', co);
c = input('Please, confirm the number of comparisions: ');
disp(' ');
disp('Give me the coefficient vectors to check orthogonality:');
for i = 1:c
eval(['c' num2str(i) '=input(''Give me the comparision ' num2str(i) ':'');']);
end;
disp(' ');
C = [];
for i = 1:c
eval(['x= c' num2str(i) ';']);
C=[C;x];
end;
O = [];AL = [];CO = [];
for i = c:-1:2
for s = i-1:-1:1
if s ~= i
Co = [i s];
CO = [CO;Co];
eval(['o' num2str(i) '= c' num2str(i) '*c' num2str(s) ''';']);
eval(['x= o' num2str(i) ';']);
O = [O;x];
if x == 0;
aL = ' It is orthogonal';
else
aL = 'It is not orthogonal';
end;
AL = [AL;aL];
end;
end;
end;
[r o] = size(CO);
fprintf('-------------------------------------\n');
disp('Comparision Orthogonality')
fprintf('-------------------------------------\n');
for ii = 1:r
fprintf(' %d %d %s\n',CO(ii,1),CO(ii,2),AL(ii,:));
end;
fprintf('-------------------------------------\n');
disp(' ');
if any(aL ~= ' It is orthogonal')
fprintf('Warning: One or more of the contrasts are nonorthogonal.');
disp(' ');
fprintf('You must to check the planned comparisions or the vector integers.');
disp(' ');
fprintf('Some of its sum of squares could be unreal: the addition don''t');
disp(' ');
fprintf('corresponds to the treatments sum of squares.');
disp(' ');
disp(' ');
end;
Ms = [];Fs = [];SSts = [];v = [];MSts = [];P = [];Si = [];
[nC,mC] = size(C);
for kC = 1:nC
t = 2;
for pr = 1:t
eval(['comparision' num2str(pr) '=input(''Give me the comparision pair ' num2str(pr) ':'');']);
end;
disp(' ');
Nt = [];
for pr = 1:t
eval(['n' num2str(pr) '=length(comparision' num2str(pr) ');']);
eval(['x= n' num2str(pr) ';']);
Nt = [Nt;x];
end;
M = [];
for pr = 1:t
eval(['m' num2str(pr) '=mean(comparision' num2str(pr) ');']);
eval(['x= m' num2str(pr) ';']);
M = [M,x];
end;
Ms=[Ms;M];
S = [];
for pr = 1:t
eval(['s' num2str(pr) '=sum(comparision' num2str(pr) ');']);
eval(['x= s' num2str(pr) ';']);
S = [S,x];
end;
CT = (sum(S))^2/sum(Nt);
D = [];
for pr = 1:t
eval(['d =((sum(comparision' num2str(pr) ').^2)/length(comparision' num2str(pr) '));']);
D = [D,d];
end;
SSt = sum(D)-CT; %sum of squares of the comparision
SSts = [SSts;SSt];
dft = t-1; %degrees of freedom of the comparision
MSt = SSt/dft; %mean of square of the comparision
MSts = [MSts;MSt];
F = MSt/MSE; %F-statistic for the comparision
Fs = [Fs;F];
v1 = dft;
v = [v;v1];
v2 = dfE;
Pr = 1 - fcdf(F,v1,v2); %probability associated to the F-statistic
P = [P;Pr];
if Pr >= alpha;
si ='NS';
else
si ='S ';
end;
Si = [Si;si];
end;
SSTO = SST + SSE; %total sum of squares
dfTO = dfT + v2; %total degrees of freedom
fprintf('---------------------------------------\n');
disp(' M e a n s ')
disp(' Comparision 1 2')
fprintf('---------------------------------------\n');
for c = 1:k-1
fprintf(' %d %.3f %.3f\n',c,Ms(c,1),Ms(c,2));
end;
fprintf('---------------------------------------\n');
disp(' ')
disp('Planned Analysis of Variance Table.')
fprintf('-------------------------------------------------------------------------------\n');
disp('SOV SS df MS F P Decision')
fprintf('-------------------------------------------------------------------------------\n');
fprintf('Treatments%14.3f%9i\n',SST,dfT);
for c = 1:k-1
fprintf(' %d %14.3f%9i%12.3f%12.3f%11.4f %s\n',(c),SSts(c),v(c),MSts(c),Fs(c),P(c),Si(c,:))
end;
fprintf('Error %18.3f%9i%12.3f\n',SSE,v2,MSE);
fprintf('Total %18.3f%9i\n',SSTO,dfTO);
fprintf('-------------------------------------------------------------------------------\n');
fprintf('With a given significance level of: %.4f\n', alpha);
disp('The results are significant (S) and/or not significant (NS).');
return;