function [pval table] = circ_hktest(alpha, idp, idq, inter, fn)
%
% [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
% process inputs
alpha = alpha(:); idp = idp(:); idq = idq(:);
if nargin < 4
inter = true;
end
if nargin < 5
fn = {'A','B'};
end
% number of groups for every factor
pu = unique(idp);
p = length(pu);
qu = unique(idq);
q = length(qu);
% number of samples
n = length(alpha);
% compute important sums for the test statistics
cn = zeros(p,q); cr = cn;
pm = zeros(p,1); pr = pm; pn = pm;
qm = zeros(q,1); qr = qm; qn = qm;
for pp = 1:p
p_id = idp == pu(pp); % indices of factor1 = pp
for qq = 1:q
q_id = idq == qu(qq); % indices of factor2 = qq
idx = p_id & q_id;
cn(pp,qq) = sum(idx); % number of items in cell
cr(pp,qq) = cn(pp,qq) * circ_r(alpha(idx)); % R of cell
end
% R and mean angle for factor 1
pr(pp) = sum(p_id) * circ_r(alpha(p_id));
pm(pp) = circ_mean(alpha(p_id));
pn(pp) = sum(p_id);
end
% R and mean angle for factor 2
for qq = 1:q
q_id = idq == qu(qq);
qr(qq) = sum(q_id) * circ_r(alpha(q_id));
qm(qq) = circ_mean(alpha(q_id));
qn(qq) = sum(q_id);
end
% R and mean angle for whole sample (total)
tr = n * circ_r(alpha);
% estimate kappa
kk = circ_kappa(tr/n);
% different formulas for different width of the distribution
if kk > 2
% large kappa
% effect of factor 1
eff_1 = sum(pr.^2 ./ sum(cn,2)) - tr.^2/n;
df_1 = p-1;
ms_1 = eff_1 / df_1;
% effect of factor 2
eff_2 = sum(qr.^2 ./ sum(cn,1)') - tr.^2/n;
df_2 = q-1;
ms_2 = eff_2 / df_2;
% total effect
eff_t = n - tr.^2/n;
df_t = n-1;
m = mean(cn(:));
if inter
% correction factor for improved F statistic
beta = 1/(1-1/(5*kk)-1/(10*(kk^2)));
% residual effects
eff_r = n - sum(sum(cr.^2./cn));
df_r = p*q*(m-1);
ms_r = eff_r / df_r;
% interaction effects
eff_i = sum(sum(cr.^2./cn)) - sum(qr.^2./qn) ...
- sum(pr.^2./pn) + tr.^2/n;
df_i = (p-1)*(q-1);
ms_i = eff_i/df_i;
% interaction test statistic
FI = ms_i / ms_r;
pI = 1-fcdf(FI,df_i,df_r);
else
% residual effect
eff_r = n - sum(qr.^2./qn)- sum(pr.^2./pn) + tr.^2/n;
df_r = (p-1)*(q-1);
ms_r = eff_r / df_r;
% interaction effects
eff_i = [];
df_i = [];
ms_i =[];
% interaction test statistic
FI = [];
pI = NaN;
beta = 1;
end
% compute all test statistics as
% F = beta * MS(A) / MS(R);
F1 = beta * ms_1 / ms_r;
p1 = 1 - fcdf(F1,df_1,df_r);
F2 = beta * ms_2 / ms_r;
p2 = 1 - fcdf(F2,df_2,df_r);
else
% small kappa
% correction factor
rr = besseli(1,kk) / besseli(0,kk);
f = 2/(1-rr^2);
chi1 = f * (sum(pr.^2./pn)- tr.^2/n);
df_1 = 2*(p-1);
p1 = 1 - chi2cdf(chi1, df_1);
chi2 = f * (sum(qr.^2./qn)- tr.^2/n);
df_2 = 2*(q-1);
p2 = 1 - chi2cdf(chi2, df_2);
chiI = f * (sum(sum(cr.^2 ./ cn)) - sum(pr.^2./pn) ...
- sum(qr.^2./qn)+ tr.^2/n);
df_i = (p-1) * (q-1);
pI = 1 - chi2pdf(chiI, df_i);
end
na = nargout;
if na < 2
printTable;
end
prepareOutput;
function printTable
if kk>2
fprintf('\nANALYSIS OF VARIANCE TABLE (HIGH KAPPA MODE)\n\n');
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');
fprintf('--------------------------------------------------------------------\n');
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);
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);
if (inter)
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);
end
fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', df_r, eff_r, ms_r);
fprintf('--------------------------------------------------------------------\n');
fprintf('%s\t\t%u\t\t%.2f', 'Total ',df_t,eff_t);
fprintf('\n\n')
else
fprintf('\nANALYSIS OF VARIANCE TABLE (LOW KAPPA MODE)\n\n');
fprintf('%s\t\t\t\t%s\t%s\t\t\t%s\n', ' ' ,'d.f.', 'CHI2', 'P-Value');
fprintf('--------------------------------------------------------------------\n');
fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{1}, df_1 , chi1, p1);
fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{2}, df_2 , chi2, p2);
if (inter)
fprintf('%s\t\t%u\t\t%.2f\t\t\t%.4f\n', 'Interaction', df_i , chiI, pI);
end
fprintf('--------------------------------------------------------------------\n');
fprintf('\n\n')
end
end
function prepareOutput
pval = [p1 p2 pI];
if na > 1
if kk>2
table = {'Source','d.f.','SS','MS','F','P-Value'; ...
fn{1}, df_1 , eff_1, ms_1, F1, p1; ...
fn{2}, df_2 , eff_2, ms_2, F2, p2; ...
'Interaction', df_i , eff_i, ms_i, FI, pI; ...
'Residual', df_r, eff_r, ms_r, [], []; ...
'Total',df_t,eff_t,[],[],[]};
else
table = {'Source','d.f.','CHI2','P-Value'; ...
fn{1}, df_1 , chi1, p1;
fn{2}, df_2 , chi2, p2;
'Interaction', df_i , chiI, pI};
end
end
end
end
Datasets
Models
