%Box, Hunter, Hunter Coagulation Times Example
%Here 24 animals (Box, Hunter, Hunter; Statistics for Experimenters,
%p. 166) are randomly allocated to 4 different diets,
%but the numbers allocated to different diets are not the same.
%The coagulation time for blood is measured for each animal.
%Are the diet-based differences significant?
clear all
close all
  disp('ANOM Coagulation Times')
  lw =2.5; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 13;


 times = [62, 60, 63, 59, ...
          63, 67, 71, 64, 65, 66, ...
          68, 66, 71, 67, 68, 68, ...
          56, 62, 60, 61, 63, 64, 63, 59];
 diets = {'dietA', 'dietA','dietA','dietA', ...
'dietB','dietB','dietB','dietB','dietB','dietB',...
'dietC','dietC','dietC','dietC','dietC','dietC',...
'dietD','dietD','dietD','dietD','dietD','dietD','dietD','dietD'};
 [p,table,stats] = anova1(times, diets,'off')
%  gnames: {4x1 cell}
%  n: [4 6 6 8]
%  source: 'anova1'
%  means: [61 66 68 61]
%  df: 20
%  s: 2.3664
ni = stats.n; NN = sum(ni);
k=length(ni);
MSE = stats.s^2;
mui=stats.means;
mu=mean(mui);

si=[];
for i=1:k
    si = [si  sqrt(  MSE * (1/ni(i) *  ((k-1)/k)^2 + ...
            (sum(1./ni) - 1/ni(i))*1/(k^2) ) )  ];
end
alpha=0.05; tquant = tinv(1 - alpha/(2*k), NN-k);

for i = 1:k
mui(i)
[mu -  tquant * si(i), mu + tquant * si(i)] 
end

%PLOTTING 
figure;
xx = (1-0.5):0.01:(k+0.5); xy = -0.5:0.01:0.5;
onex = ones(size(xx));  onexy=ones(size(xy));
% plot of grand mean
plot(xx,  onex*mu, 'linewidth', lw)
hold on
for (i = 1:k) %ANOM bands
    plot(xy+i, (mu  -  si(i)*tquant)*onexy,'r-','linewidth', lw)
    plot(xy+i, (mu  + si(i)*tquant)*onexy,'r-','linewidth', lw)
end
%treatment means as green circles
plot(1:k, mui,'o','markersize',msize, 'MarkerEdgeColor','k',...
                'MarkerFaceColor','g') 
 axis([0.5 4.5 60 69])
%print -depsc 'C:\BESTAT\ANOVA\ANOVAeps\anomct.eps'

% http://statbook.gatech.edu/
% Published under a Creative Commons Attribution-Non-Commercial-Share Alike
% 3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/

