Code covered by the BSD License

# aov2apr

### Antonio Trujillo-Ortiz (view profile)

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)
%
%
% 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
%           Universidad Autonoma de Baja California
%           Apdo. Postal 453
%           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',...
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,```