function [adTukeyAOV2] = adTukeyAOV2(X,r,alfa)
%Tukey's test of additivity for a two-way classification Analysis of Variance.
%(The two-way ANOVA design with replication can be with equal or unequal sample sizes.)
%
% Syntax: function [adTukeyAOV2] = adTukeyAOV2(X,r,a)
%
% Inputs:
% X - data matrix.
% r - two-way ANOVA design (with replication,r = 1; without replication,r = 2).
% alfa - significance level (default = 0.05).
% Outputs (table with all the completed results of the analysis)
% P - Probability that null Ho: is true (additivity).
%
% Example: From the example for a two-way (2x3) classification Analysis of Variance with
% replication (r = 1) and equal sample sizes on Box 11.1 and 13.4 of Sokal and
% Rohlf (1981, p.333-337,416-417), to test the Tukey's additivity of data with
% a significance level = 0.05.
% Factor 1
% ---------------------------------
% 1 2
% ---------------------------------
% 7.16 8.26 6.14 6.14
% 1 6.78 14.00 3.86 10.00
% 13.60 16.10 10.40 11.60
% F 8.93 9.66 5.49 5.80
% a ---------------------------------
% c 5.20 13.20 4.47 4.95
% t 2 5.20 8.39 9.90 6.49
% o 7.18 10.40 5.75 5.44
% r 6.37 7.18 11.80 9.90
% ---------------------------------
% 2 11.11 10.50 9.63 14.50
% 3 9.74 14.60 6.38 10.20
% 18.80 11.10 13.40 17.70
% 9.74 11.80 14.50 12.30
% ---------------------------------
%
% Data matrix must be:
% X=[7.16 1 1;6.78 1 1;13.6 1 1;8.93 1 1;8.26 1 1;14 1 1;16.1 1 1;9.66 1 1;
% 5.2 1 2;5.2 1 2;7.18 1 2;6.37 1 2;13.2 1 2;8.39 1 2;10.4 1 2;7.18 1 2;
% 11.11 1 3;9.74 1 3;18.8 1 3;9.74 1 3;10.5 1 3;14.6 1 3;11.1 1 3;11.8 1 3;
% 6.14 2 1;3.86 2 1;10.4 2 1;5.49 2 1;6.14 2 1;10 2 1;11.6 2 1;5.8 2 1;
% 4.47 2 2;9.9 2 2;5.75 2 2;11.8 2 2;4.95 2 2;6.49 2 2;5.44 2 2;9.9 2 2;
% 9.63 2 3;6.38 2 3;13.4 2 3;14.5 2 3;14.5 2 3;10.2 2 3;17.7 2 3;12.3 2 3];
%
% Calling on Matlab the function:
% adTukeyAOV2(X,1)
%
% Answer is:
%
% It is a procedure with replication.
% The number of levels of factor 1 are: 2
% The number of levels of factor 2 are: 3
%
% -------------------------------------------------------------
% SOV SS df MS F P
% -------------------------------------------------------------
% AB 23.926 2
% Nonaddit. 4.240 1 4.240 0.22 0.7234
% Remainder 19.686 1 19.686
% -------------------------------------------------------------
% The hypothesis of additivity is tenable.
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls and R. Castro-Valdez
% Facultad de Ciencias Marinas
% Universidad Autonoma de Baja California
% Apdo. Postal 453
% Ensenada, Baja California
% Mexico.
% atrujo@uabc.mx
%
% January 10, 2003.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls and R. Castro-Valdez. (2003).
% adTukeyAOV2: Tukey's test of additivity for a two-way classification
% Analysis of Variance. A MATLAB file.
% [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
% loadFile.do?objectId=2931&objectType=file
%
% References:
%
% Sokal, R. R. and Rohlf, F. J. (1981), Biometry. 2nd. ed.
% New-York:W. H. Freedman and Company. p. 333-337,414-417.
% Snedecor, G. W. and Cochran, W. G. (1973), Statistical Methods. 6th. ed.
% Iowa:The Iowa State Univ. Press. p. 331-334.
%
if nargin < 3,
alfa = 0.05;
end
if nargin < 2,
error('Requires at least two input arguments.');
end
a=max(X(:,2)); %number of levels factor 1.
b=max(X(:,3)); %number of levels factor 2.
indice=X(:,2);
for i=1:a
Xe=find(indice==i);
eval(['A' num2str(i) '=X(Xe,1);']);
end
indice=X(:,3);
for j=1:b
Xe=find(indice==j);
eval(['B' num2str(j) '=X(Xe,1);']);
end
C=(sum(X(:,1)))^2/length(X(:,1)); %correction term.
SSTO=sum(X(:,1).^2)-C; %total sum of squares.
dfTO=length(X(:,1))-1; %total degrees of freedom.
%procedure related to the factor 1.
A=[];
for i=1:a
eval(['x =((sum(A' num2str(i) ').^2)/length(A' num2str(i) '));']);
A=[A,x];
end
SSA=sum(A)-C; %sum of squares for the factor 1.
dfA=a-1; %degrees of freedom for the factor 1.
%procedure related to the factor 2.
B=[];
for j=1:b
eval(['x =((sum(B' num2str(j) ').^2)/length(B' num2str(j) '));']);
B=[B,x];
end
SSB=sum(B)-C; %sum of squares for the factor 2.
dfB=b-1; %degrees of freedom for the factor 2.
if r==1;
disp('It is a procedure with replication.');
disp(' ')
fprintf('The number of levels of factor 1 are:%2i\n', a);
fprintf('The number of levels of factor 2 are:%2i\n', b);
%procedure related to the treatments (T).
for i=1:a
for j=1:b
Xe=find((X(:,2)==i) & (X(:,3)==j));
eval(['T' num2str(i) num2str(j) '=X(Xe,1);']);
end
end
T=[];
for i=1:a
for j=1:b
eval(['x =((sum(T' num2str(i) num2str(j) ').^2)/length(T' num2str(i) num2str(j) '));']);
T=[T,x];
end
end
SST=sum(T)-C; %sum of squares of the treatments.
dfT=(a*b)-1; %degrees of freedom of the treatments.
%procedure related to the interaction (AB).
SSAB=SST-SSA-SSB; %sum of squares of the interaction.
dfAB=dfA*dfB; %degrees of freedom of the interaction.
MSAB=SSAB/dfAB; %sum of squares of the interaction.
%procedure related to the nonadditivity (NA) for a procedure with replication.
CA=(sum(X(:,1)))/(a*b);
AA=[];
for i=1:a
eval(['x =((sum(A' num2str(i) ')/b)-CA);']);
AA=[AA,x];
end
BA=[];
for j=1:b
eval(['x =((sum(B' num2str(j) ')/a)-CA);']);
BA=[BA,x];
end
TA=[];
for i=1:a
for j=1:b
eval(['x =sum(T' num2str(i) num2str(j) ');']);
TA=[TA,x];
end
end
TAr=reshape(TA,b,a);
for i=1:a
eval(['TAr' num2str(i) '=TAr(:,i);']);
end
QA=[];
for i=1:a
eval(['x =(TAr' num2str(i) ')''*BA'';']);
QA=[QA,x];
end
Q=QA*AA';
K=sum(AA.^2)*sum(BA.^2);
no1=[]; no2=[]; %procedure for estimation of sample size for equal or unequal
for i=1:a %sample sizes.
for j=1:b
eval(['x1 =length(T' num2str(i) num2str(j) ');']);
eval(['x2 =length(T' num2str(i) num2str(j) ')^2;']);
no1=[no1,x1]; no2=[no2,x2];
end
end
no=((1/((a*b)-1))*(sum(no1)-(sum(no2)/sum(no1))));
SSNA=Q^2/(K*no); %sum of squares of nonadditivity.
dfNA=1; %degrees of freedom of nonadditivity.
MSNA=SSNA/dfNA; %mean square of nonadditivity.
SSRem=SSAB-SSNA; %sum of squares of the remainder.
dfRem=dfAB-dfNA; %degrees of freedom of the remainder.
MSRem=SSRem/dfRem; %mean square of the remainder.
FNA=MSNA/MSRem; %F statistic for the nonadditivity.
PrNA= 1 - fcdf(FNA,dfNA,dfRem); %probability associated to the F statistic for the nonadditivity.
disp(' ')
fprintf('-------------------------------------------------------------\n');
disp(' SOV SS df MS F P')
fprintf('-------------------------------------------------------------\n');
fprintf('AB %11.4f%8i\n\n',SSAB,dfAB);
fprintf('Nonaddit. %11.4f%8i%12.3f%11.2f%9.4f\n\n',SSNA,dfNA,MSNA,FNA,PrNA);
fprintf('Remainder %11.4f%8i%12.3f\n\n',SSRem,dfRem,MSRem);
fprintf('-------------------------------------------------------------\n');
disp(' ')
if PrNA >= alfa
disp('The hypothesis of additivity is tenable.');
else
disp('The hypothesis of additivity is untenable.');
end
else r==2;
disp('It is a procedure without replication.');
disp(' ')
fprintf('The number of levels of factor 1 are:%2i\n', a);
fprintf('The number of levels of factor 2 are:%2i\n', b);
%procedure related to the nonadditivity (NA) for a procedure without replication.
SSR=SSTO-SSA-SSB; %sum of squares of residual.
dfR=dfTO-dfA-dfB; %degrees of freedom of residual.
c2A=[];
for i=1:a
eval(['x =(mean(A' num2str(i) ' - mean(X(:,1))).^2);'])
c2A=[c2A,x];
end
CA=sum(c2A);
c2B=[];
for j=1:b
eval(['x =(mean(B' num2str(j) ' - mean(X(:,1))).^2);'])
c2B=[c2B,x];
end
CB=sum(c2B);
cA=[];
for i=1:a
eval(['m =mean(A' num2str(i) ' - mean(X(:,1)));'])
cA=[cA,m];
end
cB=[];
for j=1:b
eval(['m =mean(B' num2str(j) ' - mean(X(:,1)));'])
cB=[cB,m];
end
DA=[];
for j=1:b
eval(['x =cA*B' num2str(j) ';'])
DA=[DA,x];
end
D=cB*DA';
SSNA=D^2/(sum(c2A)*sum(c2B)); %sum of squares of nonadditivity.
dfNA=1; %degrees of freedom of nonadditivity.
MSNA=SSNA/dfNA; %mean square of nonadditivity.
SSRem=SSR-SSNA; %sum of squares of the remainder.
dfRem=dfR-dfNA; %degrees of freedom of the remainder.
MSRem=SSRem/dfRem; %mean square of the remainder.
FNA=MSNA/MSRem; %F statistic for the nonadditivity.
PrNA= 1 - fcdf(FNA,dfNA,dfRem); %probability associated to the F statistic for the nonadditivity.
disp(' ')
fprintf('-------------------------------------------------------------\n');
disp(' SOV SS df MS F P')
fprintf('-------------------------------------------------------------\n');
fprintf('Residual %11.4f%8i\n\n',SSR,dfR);
fprintf('Nonaddit. %11.4f%8i%12.3f%11.2f%9.4f\n\n',SSNA,dfNA,MSNA,FNA,PrNA);
fprintf('Remainder %11.4f%8i%12.3f\n\n',SSRem,dfRem,MSRem);
fprintf('-------------------------------------------------------------\n');
disp(' ')
if PrNA >= alfa
disp('The hypothesis of additivity is tenable.');
else
disp('The hypothesis of additivity is untenable.');
end
end