% LMcs120201_4th.m
% Case Study 12.2.1, Smoking & Exercise Study
% Larsen & Marx (2006) Introduction to Mathematical Statistics 4th edition
% Page 740
% Written by Eugene.Gallagher@umb.edu 12/3/2010; revised 2/15/11
% There are at least 2 ways to analyze these data with Matlab: anova1 &
% anovan
DATA=[69 55 66 91
    52 60 81 72
    71 78 70 81
    58 58 77 67
    59 62 57 95
    65 66 79 84];
DATA=DATA;
Tdotj=sum(DATA);
Ymeandotj=mean(DATA);
[p,table,stats] = anova1(DATA);
pause
% ANOVA1 can only be used for balanced data. ANOVAN is the more
% general approach, but the data have to be restructured
multcompare(stats,'ctype','hsd','alpha',0.05)
pause
% Since a linear contrast was used below, the Scheffe procedure should be
% used for the pair-wise contrasts.
multcompare(stats,'ctype','scheffe','alpha',0.05)
pause
y=DATA(:); % convert the data into columns
group=repmat(1:4,6,1);group=group(:);
% Levene test downloaded from Matlab Central
Levenetest([y group],0.05);
% Levene's test indicates no evident problems
% Try a box-cox transformation, but 1st set up dummy variables
X=[ones(length(y),1) [ones(6,1);zeros(18,1)] ...
    [zeros(6,1);ones(6,1);zeros(12,1)] [zeros(12,1);ones(6,1);zeros(6,1)]];
PlotLogLike=1;LambdaValues=1;alpha=0.05;
[LambdaHat,LambdaInterval]=boxcoxlm(y,X,1,[-4:0.01:4]);
pause
[p,table,stats,terms] = anovan(y,group,'varnames','Smoking Level');
% Generate Figure 12.2.2 using LMex040307_4th.m as a model.
X=0:.01:4.5;
Y = fpdf(X,3,20);
plot(X,Y,'-k');
axis([0 4.5 0 0.8]);title('Figure 12.2.2','FontSize',20);
ax1=gca;
xlabel('y','FontSize',16),
ylabel('f_{F_{3,20}}(y)','FontSize',16);
ax1=gca;
set(ax1,'xtick',3.1,'FontSize',14);
hold on;
xf=3.1:.01:4.5;yf=fpdf(xf,3,20);
fill([3.1 xf 4.5]',[0 yf 0]',[.8 .8 1])
text(1,0.1,'Area=0.95','FontSize',18);
text(3.1,0.1,'Area=0.05','FontSize',18)
figure(gcf);pause
hold off;
% The following analysis uses the concept of linear contrast presented on
% page 751-758 in Larsen & Marx. The linear contrast between smokers and 
% the one non-smoking group was set a priori, so it can be tested and
% reported with an alpha level of 0.05.
format rat; LM=center(orthpoly(1:4,3,1))
LMatrix=[   -1   1/3   1/3  1/3
	    -1     0    0      1
              -3/2   -1/2   1/2   3/2
             -3/2 -11/6 -1/6 7/2
 747/310       -141/310       -909/310        303/310];
% Call's Gallagher's anova linear contrast function
anovalc(LMatrix, y, group, stats);
% To compare the slope of the 3rd linear contrast with a regression, slope
% do a linear regression.
x=[ones(length(y),1) [repmat(1,6,1);repmat(2,6,1);
   repmat(3,6,1);repmat(4,6,1)]];
[b,bint,R,RINT,STATS] = regress(y,x)
