Home > atmlab > circular > circ_hktest.m

circ_hktest

PURPOSE ^

SYNOPSIS ^

function [pval table] = circ_hktest(alpha, idp, idq, inter, fn)

DESCRIPTION ^

 [pval, stats] = circ_hktest(alpha, idp, idq, inter, fn)
   Parametric two-way ANOVA for circular data with interations.

   Input:
     alpha   angles in radians
     idp     indicates the level of factor 1 (1:p)
     idq     indicates the level of factor 2 (1:q)
     inter   0 or 1 - whether to include effect of interaction or not
     fn      cell array containing strings with the names of the factors
               

   Output:
     pval    vector of pvalues testing column, row and interaction effects
     table   cell array containg the anova table

   The test assumes underlying von-Mises distributrions.
   All groups are assumed to have a common concentration parameter k,
   between 0 and 2.

 PHB 7/19/2009 with code by Tal Krasovsky, Mc Gill University

 References:
   Harrison, D. and Kanji, G. K. (1988). The development of analysis of variance for
   circular data. Journal of applied statistics, 15(2), 197-223.

 Circular Statistics Toolbox for Matlab

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

circ_hktest.m

SOURCE CODE ^

0001 function [pval table] = circ_hktest(alpha, idp, idq, inter, fn)
0002 
0003 %
0004 % [pval, stats] = circ_hktest(alpha, idp, idq, inter, fn)
0005 %   Parametric two-way ANOVA for circular data with interations.
0006 %
0007 %   Input:
0008 %     alpha   angles in radians
0009 %     idp     indicates the level of factor 1 (1:p)
0010 %     idq     indicates the level of factor 2 (1:q)
0011 %     inter   0 or 1 - whether to include effect of interaction or not
0012 %     fn      cell array containing strings with the names of the factors
0013 %
0014 %
0015 %   Output:
0016 %     pval    vector of pvalues testing column, row and interaction effects
0017 %     table   cell array containg the anova table
0018 %
0019 %   The test assumes underlying von-Mises distributrions.
0020 %   All groups are assumed to have a common concentration parameter k,
0021 %   between 0 and 2.
0022 %
0023 % PHB 7/19/2009 with code by Tal Krasovsky, Mc Gill University
0024 %
0025 % References:
0026 %   Harrison, D. and Kanji, G. K. (1988). The development of analysis of variance for
0027 %   circular data. Journal of applied statistics, 15(2), 197-223.
0028 %
0029 % Circular Statistics Toolbox for Matlab
0030 
0031 % process inputs
0032 alpha = alpha(:); idp = idp(:); idq = idq(:);
0033 
0034 if nargin < 4
0035   inter = true;
0036 end
0037 
0038 if nargin < 5
0039   fn = {'A','B'};
0040 end
0041   
0042 
0043 % number of groups for every factor
0044 pu = unique(idp);
0045 p = length(pu);
0046 qu = unique(idq);
0047 q = length(qu);
0048 
0049 % number of samples
0050 n = length(alpha);
0051 
0052 % compute important sums for the test statistics
0053 cn = zeros(p,q); cr = cn;
0054 pm = zeros(p,1); pr = pm; pn = pm;
0055 qm = zeros(q,1); qr = qm; qn = qm;
0056 for pp = 1:p
0057     p_id = idp == pu(pp); % indices of factor1 = pp
0058     for qq = 1:q
0059         q_id = idq == qu(qq); % indices of factor2 = qq
0060         idx = p_id & q_id;
0061         cn(pp,qq) = sum(idx);     % number of items in cell
0062         cr(pp,qq) = cn(pp,qq) * circ_r(alpha(idx)); % R of cell
0063     end
0064     % R and mean angle for factor 1
0065     pr(pp) = sum(p_id) * circ_r(alpha(p_id));
0066     pm(pp) = circ_mean(alpha(p_id));
0067     pn(pp) = sum(p_id);
0068 end
0069 
0070 % R and mean angle for factor 2
0071 for qq = 1:q
0072     q_id = idq == qu(qq);
0073     qr(qq) = sum(q_id) * circ_r(alpha(q_id));
0074     qm(qq) = circ_mean(alpha(q_id));
0075     qn(qq) = sum(q_id);
0076 end
0077 
0078 % R and mean angle for whole sample (total)
0079 tr = n * circ_r(alpha);
0080 
0081 % estimate kappa
0082 kk = circ_kappa(tr/n);
0083 
0084 % different formulas for different width of the distribution
0085 if kk > 2
0086   % large kappa
0087   
0088   % effect of factor 1
0089   eff_1 = sum(pr.^2 ./ sum(cn,2)) - tr.^2/n;
0090   df_1 = p-1;
0091   ms_1 = eff_1 / df_1;
0092 
0093   % effect of factor 2
0094   eff_2 = sum(qr.^2 ./ sum(cn,1)') - tr.^2/n;
0095   df_2 = q-1;
0096   ms_2 = eff_2 / df_2;
0097 
0098   % total effect
0099   eff_t = n - tr.^2/n;
0100   df_t = n-1;
0101 
0102   m = mean(cn(:));
0103   
0104   if inter
0105 
0106     % correction factor for improved F statistic
0107     beta = 1/(1-1/(5*kk)-1/(10*(kk^2)));    
0108     
0109     % residual effects
0110     eff_r = n - sum(sum(cr.^2./cn));
0111     df_r = p*q*(m-1);
0112     ms_r = eff_r / df_r;
0113     
0114     % interaction effects
0115     eff_i = sum(sum(cr.^2./cn)) - sum(qr.^2./qn) ...
0116                   - sum(pr.^2./pn) + tr.^2/n;
0117     df_i = (p-1)*(q-1);
0118     ms_i = eff_i/df_i;
0119     
0120     % interaction test statistic
0121     FI = ms_i / ms_r;
0122     pI = 1-fcdf(FI,df_i,df_r);
0123     
0124   else
0125     
0126     % residual effect
0127     eff_r = n - sum(qr.^2./qn)- sum(pr.^2./pn) + tr.^2/n;
0128     df_r = (p-1)*(q-1);
0129     ms_r = eff_r / df_r;   
0130     
0131     % interaction effects
0132     eff_i = [];
0133     df_i = [];
0134     ms_i =[];
0135     
0136     % interaction test statistic
0137     FI = [];
0138     pI = NaN;
0139     beta = 1;
0140   end
0141   
0142   % compute all test statistics as
0143   %  F = beta * MS(A) / MS(R);
0144 
0145   F1 = beta * ms_1 / ms_r;
0146   p1 = 1 - fcdf(F1,df_1,df_r);
0147 
0148   F2 = beta * ms_2 / ms_r;
0149   p2 = 1 - fcdf(F2,df_2,df_r);
0150   
0151 else
0152   % small kappa
0153   
0154   % correction factor
0155   rr = besseli(1,kk) / besseli(0,kk);
0156   f = 2/(1-rr^2);
0157   
0158   chi1 = f * (sum(pr.^2./pn)- tr.^2/n);
0159   df_1 = 2*(p-1);
0160   p1 = 1 - chi2cdf(chi1, df_1);
0161 
0162   chi2 = f * (sum(qr.^2./qn)- tr.^2/n);
0163   df_2 = 2*(q-1);
0164   p2 = 1 - chi2cdf(chi2, df_2);
0165   
0166   chiI = f * (sum(sum(cr.^2 ./ cn)) - sum(pr.^2./pn) ...
0167             - sum(qr.^2./qn)+ tr.^2/n); 
0168   df_i = (p-1) * (q-1);
0169   pI = 1 - chi2pdf(chiI, df_i);
0170   
0171 end
0172 
0173 na = nargout;
0174 if na < 2
0175   printTable;
0176 end
0177 
0178 prepareOutput;
0179 
0180 
0181 
0182 
0183   function printTable
0184     
0185     if kk>2
0186     
0187       fprintf('\nANALYSIS OF VARIANCE TABLE (HIGH KAPPA MODE)\n\n');
0188 
0189       fprintf('%s\t\t\t\t%s\t%s\t\t%s\t\t%s\t\t\t%s\n', ' ' ,'d.f.', 'SS', 'MS', 'F', 'P-Value');
0190       fprintf('--------------------------------------------------------------------\n');
0191       fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{1}, df_1 , eff_1, ms_1, F1, p1);
0192       fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{2}, df_2 , eff_2, ms_2, F2, p2);
0193       if (inter)
0194           fprintf('%s\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', 'Interaction', df_i , eff_i, ms_i, FI, pI);
0195       end
0196       fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', df_r, eff_r, ms_r);
0197       fprintf('--------------------------------------------------------------------\n');
0198       fprintf('%s\t\t%u\t\t%.2f', 'Total   ',df_t,eff_t);
0199       fprintf('\n\n')
0200     else
0201       fprintf('\nANALYSIS OF VARIANCE TABLE (LOW KAPPA MODE)\n\n');
0202 
0203       fprintf('%s\t\t\t\t%s\t%s\t\t\t%s\n', ' ' ,'d.f.', 'CHI2', 'P-Value');
0204       fprintf('--------------------------------------------------------------------\n');
0205       fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{1}, df_1 , chi1, p1);
0206       fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{2}, df_2 , chi2, p2);
0207       if (inter)
0208           fprintf('%s\t\t%u\t\t%.2f\t\t\t%.4f\n', 'Interaction', df_i , chiI, pI);
0209       end
0210       fprintf('--------------------------------------------------------------------\n');
0211       fprintf('\n\n')
0212       
0213     end
0214     
0215 
0216   end
0217 
0218   function prepareOutput
0219     
0220     pval = [p1 p2 pI];    
0221     
0222     if na > 1
0223       if kk>2
0224         table = {'Source','d.f.','SS','MS','F','P-Value'; ...
0225                   fn{1}, df_1 , eff_1, ms_1, F1, p1; ...
0226                   fn{2}, df_2 , eff_2, ms_2, F2, p2; ...
0227                   'Interaction', df_i , eff_i, ms_i, FI, pI; ...
0228                   'Residual', df_r, eff_r, ms_r, [], []; ...
0229                   'Total',df_t,eff_t,[],[],[]};
0230       else
0231         table = {'Source','d.f.','CHI2','P-Value'; ...
0232           fn{1}, df_1 , chi1, p1;
0233           fn{2}, df_2 , chi2, p2;
0234           'Interaction', df_i , chiI, pI};
0235       end
0236     end
0237     
0238   end
0239 
0240 
0241 end
0242 
0243 
0244 
0245 
0246 
0247 
0248 
0249 
0250 
0251

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005