% LMcs110301_4th.m
% Hypothesis testing with general linear model
% Larsen & Marx Case Study 11.3.1, Page 685-688 in 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by Eugene.Gallagher@umb.edu; written 11/23/10; revised 11/23/10
DATA=[3900 256.9;3350 211.6;3220 238.1;3220 211.8;2790 194.1;2780 124.5
      2770 187.3; 2290 110.5;2160 233.1;1890 150.3;1810 124.7;1800 41.2
      1770 182.1;1700 118.1;1680 31.9;1510 114.3;1500 144.9;1410 59.7
      1270 126.9;1200 43.9;1090 136.3];
[r,c]=size(DATA);
 X=[ones(r,1) DATA(:,1)];Y=DATA(:,2);
B=X\Y;
fprintf('\nThe Y intercept = %8.3f and slope = %5.3f\n',B(1),B(2));
Yest=X*B;
Resid=Y-Yest;
% sort by the independent variable so the lines are smooth
[Xs,k]=sort(X(:,2));
% sort the estimated values using the index vector k:
Ys=Yest(k);
plot(X(:,2),Y,'o',Xs,Ys,'m',...
                        'LineWidth',2,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',8) 
xlabel('Cigarette consumption per adult per year','FontSize',20),
ylabel('CHD Mortality per 100,000 (ages 35-64)','FontSize',20)
title('Case Study 11.3.1','FontSize',22)
figure(gcf)
pause
plot(Yest,Resid,'o',...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',8) 
xlabel('Expected value','FontSize',20),ylabel('Residual','FontSize',20);
title('Case Study 11.3.1','FontSize',22)
grid;figure(gcf)
% Now use Matlab's regress
alpha=0.05;
[b,bint,r,rint,stats] = regress(Y,X,alpha);
% The help regress notation produces the following components of the
% stats variable:
% [B,BINT,R,RINT,STATS] = REGRESS(Y,X) returns a vector STATS containing, in
%    the following order, the R-square statistic, the F statistic and p value
%    for the full model, and an estimate of the error variance.
fprintf('R^2=%4.1f%%; F for overall regression =%4.1f with p-value=%7.2g;\n',...
    stats(1)*100,stats(2),stats(3));
df=length(Y)-2;
% Now, test whether the slope is different from 0. The F statistic
% produced by regress with 1 and df is the square of the t statistic
% in the t test for slope =1.
P=2*(1-tcdf(sqrt(stats(2)),df));
fprintf('t test for B_1=0 is %5.2f with % 3.0f df and p-value=%7.2g\n',...
    sqrt(stats(2)),df,P)
% Generate a 95% CI for s^2, also called Root Mean Square Error or RMSE:
% See Larsen & Marx page 691 for equation. Change alpha above for CI's
% other than 95% CI
L95CIRMSE=df*stats(4)/chi2inv(1-alpha/2,df);
U95CIRMSE=df*stats(4)/chi2inv(alpha/2,df);
fprintf('RMSE= %7.2f with %3.1f%% CI= [%7.2f %7.2f]\n',stats(4),(1-alpha)*100,L95CIRMSE,...
    U95CIRMSE)
