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