CircularANOVA - One or two-way ANOVA on circular data.
Circular ANOVA can be performed using one of three methods with different
assumptions and limitations:
1) Watson & Williams (1956)
The data are drawn from von Mises distributions with equal concent-
ration parameters, the significance levels rely on X² approximations,
the extension to multi-way analysis is difficult to interpret because
interaction factors can be negative
2) Harrison & Kanji (1986)
The data are drawn from von Mises distributions with equal concent-
ration parameters, the significance levels rely on X² approximations,
the test can be extended to multi-way analysis (positive interaction
factors), the statistics are affected by both location and dispersion
3) Anderson & Wu (1995)
No distributional assumptions, the significance levels are determined
using a bootstrap approach, the test can be extended to multi-way
analysis
USAGE
[p,F] = CircularANOVA(angles,groups,method)
[p,F] = CircularANOVA(angles,[factor1 factor2],method)
angles angles in radians
groups group indices
factor1 levels for factor 1
factor2 levels for factor 2
method method to compute ANOVA:
'ww' Watson & Williams (1956) (default)
'l2' Harrison & Kanji (1988), L2 distance
'lr' Anderson & Wu (1995), likelyhood ratio
SEE
See also CircularVariance, CircularConfidenceIntervals, Concentration,
ConcentrationTest.
REFERENCE
Anderson & Wu (1992) Decomposition of variation in directional data,
IIQP Report RR-92-10.0001 function [p,F] = CircularANOVA(angles,factors,method) 0002 0003 %CircularANOVA - One or two-way ANOVA on circular data. 0004 % 0005 % Circular ANOVA can be performed using one of three methods with different 0006 % assumptions and limitations: 0007 % 0008 % 1) Watson & Williams (1956) 0009 % The data are drawn from von Mises distributions with equal concent- 0010 % ration parameters, the significance levels rely on X² approximations, 0011 % the extension to multi-way analysis is difficult to interpret because 0012 % interaction factors can be negative 0013 % 0014 % 2) Harrison & Kanji (1986) 0015 % The data are drawn from von Mises distributions with equal concent- 0016 % ration parameters, the significance levels rely on X² approximations, 0017 % the test can be extended to multi-way analysis (positive interaction 0018 % factors), the statistics are affected by both location and dispersion 0019 % 0020 % 3) Anderson & Wu (1995) 0021 % No distributional assumptions, the significance levels are determined 0022 % using a bootstrap approach, the test can be extended to multi-way 0023 % analysis 0024 % 0025 % USAGE 0026 % 0027 % [p,F] = CircularANOVA(angles,groups,method) 0028 % [p,F] = CircularANOVA(angles,[factor1 factor2],method) 0029 % 0030 % angles angles in radians 0031 % groups group indices 0032 % factor1 levels for factor 1 0033 % factor2 levels for factor 2 0034 % method method to compute ANOVA: 0035 % 'ww' Watson & Williams (1956) (default) 0036 % 'l2' Harrison & Kanji (1988), L2 distance 0037 % 'lr' Anderson & Wu (1995), likelyhood ratio 0038 % 0039 % SEE 0040 % 0041 % See also CircularVariance, CircularConfidenceIntervals, Concentration, 0042 % ConcentrationTest. 0043 % 0044 % REFERENCE 0045 % 0046 % Anderson & Wu (1992) Decomposition of variation in directional data, 0047 % IIQP Report RR-92-10. 0048 0049 % Copyright (C) 2009-2011 by Michaël Zugaro 0050 % 0051 % This program is free software; you can redistribute it and/or modify 0052 % it under the terms of the GNU General Public License as published by 0053 % the Free Software Foundation; either version 3 of the License, or 0054 % (at your option) any later version. 0055 0056 warning('This function is not fully implemented yet (not all methods are available).'); 0057 0058 p = []; 0059 F = []; 0060 if isempty(angles), return; end 0061 0062 nIterations = 500; 0063 0064 % Check parameters 0065 if nargin < 2, 0066 error('Incorrect number of parameters (type ''help <a href="matlab:help CircularANOVA">CircularANOVA</a>'' for details).'); 0067 end 0068 if ~isdvector(angles), 0069 error('Incorrect angles (type ''help <a href="matlab:help CircularANOVA">CircularANOVA</a>'' for details).'); 0070 end 0071 isradians(angles); 0072 0073 if nargin == 2, 0074 method = 'ww'; 0075 else 0076 method = lower(method); 0077 if ~isastring(method,'ww','l2','lr'), 0078 error(['Unknown ''' method ''' method (type ''help <a href="matlab:help CircularANOVA">CircularANOVA</a>'' for details).']); 0079 end 0080 end 0081 0082 A = exp(j*angles); 0083 R = abs(mean(A)); 0084 N = length(A); 0085 mu = angle(mean(A)); 0086 0087 if isvector(factors), 0088 % One-way ANOVA 0089 groups = factors(:); 0090 groupIDs = unique(groups); 0091 if any(diff(groupIDs)~=1), 0092 error('Incorrect groups numbers (should be 1..N).'); 0093 end 0094 switch(method), 0095 case 'ww', 0096 % Watson & Williams (1956) 0097 q = length(groupIDs); 0098 for i = 1:length(groupIDs), 0099 ok = groups == groupIDs(i); 0100 Ri(i) = abs(mean(A(ok))); 0101 Ni(i) = sum(ok); 0102 ki(i) = Concentration(angles(ok)); 0103 end 0104 SSw = N - sum(Ni.*Ri); % ignore 2k factor 0105 SSb = sum(Ni.*Ri) - N*R; % ignore 2k factor 0106 F = (N-q)/(q-1)*SSb/SSw; 0107 k = sum(ki.*Ni)/N; 0108 if 2 < k && k < 10, 0109 % Improve X² approximation for moderate concentration (Stephens, 1969) 0110 F = F * (1+3/(8*k)); 0111 end 0112 p = 1-fcdf(F,q-1,N-q); 0113 otherwise, 0114 error('One-way ANOVA: only the ''ww'' method is implemented.'); 0115 end 0116 else 0117 % Two-way ANOVA 0118 factor1 = factors(:,1); 0119 factor2 = factors(:,2); 0120 levels1 = unique(factor1); 0121 if any(diff(levels1)~=1), 0122 error('Incorrect levels for factor 1 (should be 1..N).'); 0123 end 0124 levels2 = unique(factor2); 0125 if any(diff(levels2)~=1), 0126 error('Incorrect levels for factor 2 (should be 1..N).'); 0127 end 0128 switch(method), 0129 case 'lr', 0130 % Anderson & Wu (1995), location-only 0131 % Check requirements: 2x2 balanced design 0132 if length(levels1)~=2, 0133 error('Two-way ANOVA: factors can only have 2 levels (2x2 design).'); 0134 end 0135 if length(levels2)~=2, 0136 error('Two-way ANOVA: factors can only have 2 levels (2x2 design).'); 0137 end 0138 m = Accumulate([factor1 factor2]); 0139 if any(m~=m(1)), 0140 error('Two-way ANOVA: all cells must have the same count (balanced design).'); 0141 end 0142 m = m(1); 0143 % F values for observations 0144 [Fa,Fb,Fab] = lo(A,m,factor1,factor2); 0145 % F values for surrogate data 0146 for i = 1:nIterations, 0147 p = randperm(N); 0148 [sFa(i),sFb(i),sFab(i)] = lo(A(p),m,factor1,factor2); 0149 end 0150 p(1) = 1-npcdf(sFa,Fa); 0151 p(2) = 1-npcdf(sFb,Fb); 0152 p(3) = 1-npcdf(sFab,Fab); 0153 F = [Fa;Fb;Fab]; 0154 otherwise, 0155 error('Two-way ANOVA: only the ''lr'' method is implemented.'); 0156 end 0157 end 0158 0159 0160 %----------------------------------------------------------------------------------------------------------------- 0161 % Compute F values for 2-way ANOVA (location only method) 0162 %----------------------------------------------------------------------------------------------------------------- 0163 0164 function [Fa,Fb,Fab] = lo(A,m,factor1,factor2) 0165 0166 % Mean angle 0167 theta = angle(mean(A)); 0168 levels1 = unique(factor1); 0169 levels2 = unique(factor2); 0170 0171 % Mean angles for each level of factor 1 0172 for i = 1:length(levels1), 0173 theta_i(i) = angle(mean(A(factor1==i))); 0174 end 0175 0176 % Mean angles for each level of factor 2 0177 for i = 1:length(levels2), 0178 theta_j(i) = angle(mean(A(factor2==i))); 0179 end 0180 0181 % Mean angles for each level of factor 1 and factor 2 0182 for i1 = 1:length(levels1), 0183 for i2 = 1:length(levels2), 0184 theta_ij((i1-1)*length(levels2)+i2) = angle(mean(A(factor1==i1&factor2==i2))); 0185 end 0186 end 0187 t_ij = repmat(theta_ij,1,m); 0188 t_ij = t_ij'; 0189 t_ij = t_ij(:); 0190 0191 % Interaction angles 0192 theta_11 = mean(A(factor1==1&factor2==1)); 0193 theta_22 = mean(A(factor1==2&factor2==2)); 0194 theta_12 = mean(A(factor1==1&factor2==2)); 0195 theta_21 = mean(A(factor1==2&factor2==1)); 0196 phi_l(1) = angle(mean([theta_11 theta_22])); 0197 phi_l(2) = angle(mean([theta_12 theta_21])); 0198 0199 % Sum of squares 0200 SSa = 4*m*sum(1-cos(theta_i-theta)); 0201 SSb = 4*m*sum(1-cos(theta_j-theta)); 0202 SSab = 4*m*sum(1-cos(phi_l-theta)); 0203 SSr = 2*sum(1-cos(angle(A)-t_ij)); 0204 0205 % F values 0206 Fa = SSa/SSr; 0207 Fb = SSb/SSr; 0208 Fab = SSab/SSr;