Code covered by the BSD License  

Highlights from
aov2apr

aov2apr

by

 

15 Apr 2013 (Updated )

Two-way analysis of variance with a planned (a priori) factor.

aov2apr.m
function aov2apr(x,m,z)
%AOV2APR Two-way analysis of variance with a planned (a priori) factor.
% The two-way analysis of variance with a planned (a priori) factor follows
% the next statistical fundamentals (degrees of freedom are between 
% parentheses),
%
%             SSTO = SSM + SSE;  (N-1 = k-1 + N-k; k = a*b)
%
% SSTO, total sum of squares corrected; SSM, model sum of squares;SSE,
% error sum of squares (within). [N=Total data;a=Factor 1 levels;b=Factor 2
% levels].
%
%             SSM = SSA + SSB + SSAB;  [k-1 = a-1 + b-1 + (a-1)*(b-1)]  
%
% SSA, Factor A sum of squares (Model I or II); SSB, Factor B sum of 
% squares (Model I); SSAB, interaction sum of squares.
%
%             SSB = SS1 + SS2 +...+ SSC;  (b-1 = 1 + 1 +...+ 1)
%
% SS1,...,SSC, contrast sum of squares.
%                           _
%                          \
%             SSj = Sj^2 / /_ cj^2/nj; j=1,2,...,C
% 
%             Sj^2 = c1*M1 + c2*M2 +...+ cC*MC
%
% cj, contrast vector; Mj, means vector; nj, sample size vector.
%
% If there is no interaction, it is better to perform two separate one-way
% ANOVAs. Here, factor A and B with a Model I/II and I (panned), 
% respectively. An interaction effect is a change in the simple main effect
% of one variable over levels of the second.
%
% In a fixed-model (I), the treatment levels are deliberately selected and 
% will remain constant from one replication to another. Generalization of 
% such a model can be made only to the levels tested. Meanwhile, in a 
% random-model (II), the treatment levels are randomly selected and if we 
% replicated the study we would again choose the levels randomly and would
% most likely have a whole new set of levels. Results can be generalized 
% to the population levels from which the levels of the independent 
% variable were randomly selected.
%
% Planned contrasts are particular comparisions that the researcher
% is interested in examining prior to data collection or in certain
% specific contrasts, but not in the omnibus F test that examines all
% possible contrasts (Lomax, 2007). Fewer planned comparisions are usually 
% conducted (k - 1), due to their specificity. Each comparision 
% corresponds to an independent hypothesis. To be sure that this is
% satisfied, we must to have a family of orthogonal planned comparisons.
% Orthogonal contrasts are special in the sense that each contrast provides
% totally unique, nonredundant information about the means under study. If 
% the two contrasts are orthogonal, then the sum of the cross products will
% be equal to '0'. Otherwise, the contrasts are nonorthogonal.
%
% AOV2APR treats NaN values as missing values, and removes them.
%
% Syntax: function aov2apr(x,m,z)
%
% Inputs:
%      x  data nx3 matrix (Col 1 = data;Col 2 = factor 1 code [Model I/II];
%          Col 3 = factor 2 code [Model I-planned])
%      m - contrast matrix (planned)
%      z - factor 1 fixed=1;factor 1 random=2
%
% Output:
%        - Summary statistics from the analysis
%
% Taking the data of the numerical example given in the Psychology Honours:
% Statistics for Psychological Research (psy400w:Dr. Colin Tredoux) 
% handout, University of Cape Town, South Africa. Internet site
% [http://web.uct.ac.za/depts/psychology/psy400w/stats/anova/multcomp.pdf]
%
% With an theoretical modification to adjusts to our hypothetical model.
% Independent variable 1 are random effects (Model II).
%
% Data are:
%         -----------------------------------------
%                  1                     2
%         -----------------------------------------
%          1   2   3   4   5     1   2   3   4   5
%         -----------------------------------------
%         33  33   6  11  22    37  20   8   6  23
%         31  19   0   4  15    42  22  12  12  21
%         33  34   0  12  30    37  24   5  22  27
%         38  21  10  15  24
%         -----------------------------------------
%
% Contrasts matrix is:
%
%     C1: 3     3    -2    -2    -2  [Ho: mean(1,2) = mean(3,4,5)]
%     C2: 1    -1     0     0     0  [Ho: mean(1) = mean(2)]
%     C3: 0     0     1     1    -2  [Ho: mean(3,4) = mean(5)]
%     C4: 0     0     1    -1     0  [Ho: mean(3) = mean(4)]
%   
% Input data:
%   x=[33 1 1;31 1 1; 33 1 1;38 1 1;33 1 2;19 1 2;34 1 2;21 1 2;6 1 3; 
%   0 1 3;0 1 3;10 1 3;11 1 4; 4 1 4;12 1 4; 15 1 4;22 1 5;15 1 5;30 1 5;
%   24 1 5;37 2 1;42 2 1;37 2 1;20 2 2;22 2 2;24 2 2;8 2 3;12 2 3;5 2 3;
%   6 2 4;12 2 4;22 2 4;23 2 5;21 2 5;27 2 5];
%
%   m=[3 3 -2 -2 -2;1 -1 0 0 0;0 0 1 1 -2;0 0 1 -1 0];
%
% Calling on Matlab the function: 
%                aov2apr(x,m,2)
%
% Answer is:
%
% Recall that Factor A can be fixed (Model I) or random (Model II), but
% the Factor B must be necessarily the fixed and planned (Model I)
% 
% Are you sure that is it (1) or not (2)?: 1
% 
% Factor A number of levels are: 2
% Factor B number of levels are: 5
%
% -------------------------------------
% 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
% -------------------------------------
% 
% All the planned hypothesis are independent (orthogonal)
%     
% -------------------------------------------------------------
% SV               SS         df      MS          F        P
% -------------------------------------------------------------
% Model       3990.769         9
% A             23.336         1    23.336      0.896   0.3975
% B           3863.257         4
%   1         2346.686         1  2346.686     88.632   0.0000
%   2          434.571         1   434.571     16.413   0.0004
%   3          961.929         1   961.929     36.331   0.0000
%   4          120.071         1   120.071      4.535   0.0432
% AB           104.176         4    26.044      0.984   0.4344
% Error        661.917        25    26.477
% -------------------------------------------------------------
% Total       4652.686        34
% -------------------------------------------------------------
% If Factor A is Model I (fixed effects) and results significative.
% You must do a posteriori tests.
% If Factor A is Model II (random effects) and results significative.
% You must estimate the variance components.
% If interation test results not significative.
% It is recommended do two separate one-way ANOVA tests.
%
% 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.edu.mx
%
% Copyright (C) April 13, 2013. 
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A. and R. Hernandez-Walls. (2013). aov2apr: 
%    Two-way analysis of variance with a planned (a priori) factor. 
%    [WWW document]. URL http://www.mathworks.com/matlabcentral/
%    fileexchange/41305-aov2apr
%
% Reference:
% Lomax, R. G. (2007). An Introduction to Statistical Concepts. 2nd Ed.
%              Lawrence Erlbaum Associates, Inc. Pub. Mahwah:NJ
%

if  nargin < 3,
    error('aov2apr:TooFewInputs', ...
          'AOV2APR requires three input arguments.');
end

X = x;

%Remove NaN values, if any
X = X(~any(isnan(X),2),:);

disp('Recall that Factor A can be fixed (Model I) or random (Model II), but');
disp('the Factor B must be necessarily the fixed and planned (Model I)');
disp(' ');
op = input('Are you sure that is it (1) or not (2)?: ');
disp(' ');
if op == 1;
    a = max(X(:,2));
    b = max(X(:,3));
    fprintf('Factor A number of levels are: %i\n',a);
    fprintf('Factor B number of levels are: %i\n\n',b);
    
    h = b - 1; %number of planned hypothesis

    c = size(m,1);
    if c ~= h,
        error('stats:aov2apr:BadData','Error: M matrix must have %g rows.%3.2f\n',c);
    end
    
    CT = (sum(X(:,1)))^2/length(X(:,1)); %correction term
    SSTO = sum(X(:,1).^2) - CT; %total sum of squares
    dfTO = length(X(:,1)) - 1; %total degrees of freedom
    
    indice = X(:,2);
    for i = 1:a
        Xe = find(indice==i);
        d(i).X = X(Xe,1);
        d(i).nA = length(d(i).X);
        d(i).sA = ((sum(d(i).X))^2)/d(i).nA;
    end
    nA=cat(1,d.nA);sA=cat(1,d.sA);
    
    SSA = sum(sA) - CT; %sum of squares factor A
    dfA = a - 1; %degrees of freedom factor A
    MSA = SSA/dfA; %mean of squares factor A 
    
    indice = X(:,3);
    for j = 1:b
        Xe = find(indice==j);
        d(j).X = X(Xe,1);
        d(j).tB = sum(d(j).X);
        d(j).nB = length(d(j).X);
        d(j).sB = ((d(j).tB)^2)/d(j).nB;
    end
    tB = cat(1,d.tB);nB=cat(1,d.nB);sB=cat(1,d.sB);
    
    SSB = sum(sB) - CT; %sum of squares factor A
    dfB = b - 1; %degrees of freedom factor A
   
    for i = 1:a
        for j = 1:b
            Xe = find((X(:,2)==i) & (X(:,3)==j));
            d(i,j).X = X(Xe,1);
            d(i,j).nM = length(d(i,j).X);
        d(i,j).sM = ((sum(d(i,j).X))^2)/d(i,j).nM;
        end
    end
    nM=cat(1,d.nM);sM=cat(1,d.sM);
    
    SSM = sum(sM) - CT; %model sum of squares
    dfM = (a*b) - 1; %model degrees of freedom

    SSAB = SSM - SSA - SSB; %interaction sum of squeres
    dfAB = dfA*dfB; %interaction degrees of freedom
    MSAB = SSAB/dfAB; %interaction mean of square

    SSE = SSTO - SSM; %error sum of squeres
    dfE = dfTO - dfM; %error degrees of freedom
    MSE = SSE/dfE; %error mean of squeres
    
    r = c*(c-1)/2;
    C = [];
    for i = c:-1:2
        for s = i-1:-1:1
            x = [i s];
            for l = 1:r
                if s ~= i
                    d(i).CO = m(i,:);
                    d(s).CO = m(s,:);
                    d(l).O = d(i).CO * d(s).CO';
                    if d(l).O == 0
                        d(l).AL = '    It is orthogonal';
                    else
                        d(l).AL = 'It is not orthogonal';
                    end
                end
            end
            C = [C;x];
            AL = cat(1,d.AL);
        end
    end
    O = cat(1,d.O);
    
    fprintf('-------------------------------------\n');
    disp('Comparision            Orthogonality')
    fprintf('-------------------------------------\n');
    for ii = 1:r
        fprintf(' %d    %d          %s\n',C(ii,1),C(ii,2),AL(ii,:));
    end
    fprintf('-------------------------------------\n');
    disp(' ');
    
    AL = cellstr(AL);
    ix = find(strcmp(AL, 'It is not orthogonal'));
    if isempty(ix)
        disp('All the planned hypothesis are independent (orthogonal)')
    else
        disp('At least one of the planned hypothesis are not independent',...
            '(ortogonal). Check please.')
    end
    
    %Next machine it is the general procedure (balanced and unbalanced
    %designs) used for the orthogonal comparisions 
    for k = 1:c
        d(k).S = m(k,:);
        d(k).SC = d(k).S * (tB./nB);
        d(k).nB = nB(k);
        d(k).D = sum((d(k).S.*d(k).S)./d(k).nB);
        d(k).SS = (d(k).SC)^2/d(k).D;
        d(k).FT = d(k).SS./MSE;
        d(k).PT = 1 - fcdf(d(k).FT,1,dfE);
    end
    S = cat(1,d.S);SC = cat(1,d.SC);D = cat(1,d.D);SS = cat(1,d.SS);
    FT = cat(1,d.FT);PT = cat(1,d.PT);
    
    %An alternative machine can be used for the special case of balanced
    %orthogonal comparisions
    %for k = 1:c
    %    d(k).S = m(k,:);
    %    d(k).SSS = d(k).S * d(k).S';
    %    d(k).SC = d(k).S * tB;
    %    d(k).nB = nB(k);
    %    d(k).D = d(k).nB * d(k).S * d(k).S';
    %    d(k).SS = (d(k).SC)^2/d(k).D;
    %    d(k).FT = d(k).SS./MSE;
    %    d(k).PT = 1 - fcdf(d(k).FT,1,dfE);
    %end
    %S = cat(1,d.S);SC = cat(1,d.SC);D = cat(1,d.D);SS = cat(1,d.SS);
    %FT = cat(1,d.FT);PT = cat(1,d.PT); 
        
    dfST = ones(c,1);

    if z == 1
        FA = MSA/MSE;
        PrA = 1 - fcdf(FA,dfA,dfE);
    else (z == 2);
        FA = MSA/MSAB;
        PrA = 1 - fcdf(FA,dfA,dfAB);
    end
    
    FB = SS./MSE;
    FAB = MSAB/MSE;
    PrB = 1 - fcdf(FB,dfB,dfE);
    PrAB = 1 - fcdf(FAB,dfAB,dfE);

    disp('     ')
    disp('-------------------------------------------------------------')
    disp('SV               SS         df      MS          F        P')
    disp('-------------------------------------------------------------')
    fprintf('Model    %11.3f%10i\n',SSM,dfM);
    fprintf('A        %11.3f%10i%10.3f%11.3f%9.4f\n',SSA,dfA,MSA,FA,PrA);
    fprintf('B        %11.3f%10i\n',SSB,dfB);
    for ml = 1:c
        fprintf('  %d%17.3f%10d%10.3f%11.3f%9.4f\n',ml,SS(ml),dfST(ml),SS(ml),FT(ml),PT(ml))
    end
    fprintf('AB       %11.3f%10i%10.3f%11.3f%9.4f\n',SSAB,dfAB,MSAB,FAB,PrAB);
    fprintf('Error %14.3f%10i%10.3f\n',SSE,dfE,MSE);
    disp('-------------------------------------------------------------')
    fprintf('Total %14.3f%10i\n',SSTO,dfTO);
    disp('-------------------------------------------------------------')
    disp('If Factor A is Model I (fixed effects) and results significative.')
    disp('You must do a posteriori tests.')
    disp('If Factor A is Model II (random effects) and results significative.')
    disp('You must estimate the variance components.')
    disp('If interation test results not significative.')
    disp('It is recommended do two separate one-way ANOVA tests.')
    disp('     ')
else
    error('WARNING: You must to correct the input factor-data matrix.');
end

return,

Contact us