%LMcs130203_4th.m
% Larsen & Marx (2006) Introduction to Mathematical Statistics 4th edition
% Case Study 13.2.3 Page 778-779, The Transylvannia5 effect
% An example of linear contrasts for Randomized Block data.
% Written by Eugene.Gallagher@umb.edu 2/15/11
% Calls Trujillo-Ortiz et al. adTukeyAOV2.m from Matlab file central
% Tukey's test for additivity.
fprintf('\nAnalysis of Case Study 13.2.3: The Transylvannia Effect\n')
DATA=[6.4 5.0 5.8
    7.1 13.0 9.2
    6.5 14.0 7.9
    8.6 12 7.7
    8.1 6 11
    10.4 9 12.9
    11.5 13.0 13.5
    13.8 16.0 13.1
    15.4 25.0 15.8
    15.7 13.0 13.3
    11.7 14.0 12.8
    15.8 20 14.5];
[R,C]=size(DATA);
boxplot(DATA,'sym','r*','labels',...
{'Before Full Moon','During Full Moon','After Full Moon'});
figure(gcf);pause

plot(DATA');
ax1=gca;
set(ax1,'xtick',[1 2 3]);
set(ax1,'XtickLabel',...
   {'Before Full Moon','During Full Moon','After Full Moon'},'FontSize',9);
legend('Au','Se','Oc','Nv','De','Ja','Fe','Mr','Ap','My','Jn','Jl');
xlabel('Moon Phase');ylabel('Hospital Admission Rates');
figure(gcf);pause
pause % Needed so can see the box plots
sumDATA=sum(DATA);
meanDATA=mean(DATA);
% Since the design is balanced, either anova2 or anovan can be used.
[p,table,stats]=anova2(DATA,1)
% This will produce the ANOVA table as a figure. The results are
% printed out in table. stats could be sued for multcompare
pause
% The data can also be analyzed using anovan, producing identical results.
% anovan allows labeling of the ANOVA table.
y=DATA(:);
% convert the data into columns; drop the NaN elements
g1=repmat(['Au';'Se';'Oc';'Nv';'De';'Ja';'Fe';'Mr';'Ap';'My';'Jn';'Jl'],3,1);
g2=[repmat('BFM',12,1);repmat('DFM',12,1);repmat('AFM',12,1)];
group=[repmat(1,12,1) repmat(2,12,1) repmat(3,12,1)];
% find and delete any NaN elements, if any
i=~isnan(y);y=y(i);g1=g1(i,:);g2=g2(i,:);group=group(i);
% Use Trujillo-Ortiz's Levenestest
levenetest([y group],0.05);
[p,table,stats] = anovan(y,{g1 g2},'model','linear',...
    'varnames',{'Months';'Lunar Cycles'})
% This is Example 13.2.1, comparing treatments
[c,m,h] = multcompare(stats,'ctype','tukey-kramer','display','on',...
    'dimension',2);
fprintf('Pairwise Difference \tLower 95\tEstimate\tUpper 95\n');
fprintf('\t\t%1.0f - %1.0f\t\t\t   %4.1f\t\t%4.1f\t\t%4.1f\n',c')
figure(h);
title(' '),xlabel(' ');xlabel('Hospital Admission Rates');
title('Case Study 13.2.3'),
figure(gcf)
pause;

% Since the treatments were not replicated within blocks, a test
% for additivity should be run. If replicates were available a formal
% block by interaction test could have been run. 
r1=repmat([1:R]',C,1);
r2=[repmat(1,R,1);repmat(2,R,1);repmat(3,R,1)];
X=[y r1 r2];
adTukeyAOV2(X,2,0.05)
% Note that there is evidence (p=0.046) to reject the additivity assumption

fprintf('\nCheck the additivity assumption with just 2 groups:\n')
% Reanalyze the data pooling 2 non-full moon periods.
D=mean(DATA(:,[1 3])')';D=[D DATA(:,2)];
plot(D');
ax1=gca;
set(ax1,'xtick',[1 2]);
set(ax1,'XtickLabel',...
   {'Not Full Moon','Full Moon'},'FontSize',9);
legend('Au','Se','Oc','Nv','De','Ja','Fe','Mr','Ap','My','Jn','Jl');
xlabel('Moon Phase');ylabel('Hospital Admission Rates');
figure(gcf);pause
pause % Needed so can see the box plots
[r,c]=size(D);
r1=repmat([1:r]',c,1);
r2=[repmat(1,r,1);repmat(2,r,1)];
X=[D(:) r1 r2];
adTukeyAOV2(X,2,0.05);
% p=0.0367; so still a strong interaction evident.

[p2,table2,stats2] = anovan(D(:),[r1 r2],'model','linear',...
    'varnames',{'Months';'Lunar Cycles'})


% Not covered in Larsen & Marx, but now it is possible to test formally
% for the interaction term.
Y=DATA(:);
G1=repmat(['Au';'Se';'Oc';'Nv';'De';'Ja';'Fe';'Mr';'Ap';'My';'Jn';'Jl'],3,1);
% set two groups: Not full moon and During Full moon
G2=[repmat('NFM',12,1);repmat('DFM',12,1);repmat('NFM',12,1)];
Group=[repmat(1,12,1) repmat(2,12,1) repmat(1,12,1)];
% find and delete any NaN elements, if any
i=~isnan(Y);Y=Y(i);G1=G1(i,:);G2=G2(i,:);Group=Group(i);
% Use Trujillo-Ortiz's Levenestest
levenetest([Y Group],0.05);
[p,table,stats] = anovan(Y,{G1 G2},'model',2,...
    'varnames',{'Months';'Lunar Cycles'})
% There should be no formal analysis of main effects of the main effects if
% I was taught in my graduate statistics class that if there is a 
% significant interaction, show the interactions in an effects plot,discuss
% them and end the analysis.

% If there were no interactions, this would be a valid post hoc analysis:
% The following analysis uses the concept of linear contrast presented on
% page 751-758 in Larsen & Marx. The linear contrast between the full moon 
% period and the other two phases was set a priori, so it can be tested and
% reported with an alpha level of 0.05.
LMatrix=[-1/2 1 -1/2];
planned=0;
anovalc(LMatrix, y, group,stats,planned)