% Compare ANOVA for linear and quadratic LSQ fits
% CompareANOVAlouismuddata.m  L Braile 3/2/05
% enter x and y data (louismud data from Davis (2002, Fig.4.27)
% Depth (ft)	Moisture (g H20/100 g solid)	Core number
X = [
0	124	1
5	78	1
10	54	1
15	35	1
20	30	1
25	21	1
30	22	1
35	18	1
0	137	2
5	84	2
10	50	2
15	32	2
20	28	2
25	24	2
30	23	2
35	20	2];
x = X(:,1); % create x vector
y = X(:,2); % create y vector
%
% for linear fit
n = length(x);
[P,S] = polyfit(x,y,1); % LSQ straight line fit
yt1 = polyval(P,x); % calculate theoretical values of y (y-hat)
ybar = mean(y);
SSR1 = sum((yt1-ybar).^2) % calculate and print SSR
MSR1 = SSR1/1 % note d.f.
SSE1 = sum((y-yt1).^2)
MSE1 = SSE1/(n-2) % note d.f.
SST1 = sum((y-ybar).^2)
R1 = sqrt(SSR1/SST1)
ftest1 = MSR1/MSE1
fvalue1 = finv(0.95,1,n-2) % f-test value from tables
p1 = fcdf(ftest1,1,n-2);
p1 = 1-p1 % Calculates the probability of
% a value as high as ftest1 resulting from a set of
% random inputs.
%
% for quadratic fit
[P,S] = polyfit(x,y,2); % LSQ quadratic line fit
yt2 = polyval(P,x); % calculate theoretical values of y (y-hat)
ybar = mean(y);
SSR2 = sum((yt2-ybar).^2) % calculate and print SSR
MSR2 = SSR2/2 % note d.f.
SSE2 = sum((y-yt2).^2)
MSE2 = SSE2/(n-3) % note d.f.
SST2 = sum((y-ybar).^2)
R2 = sqrt(SSR2/SST2)
% ftest to see significance of the quadratic fit
ftest2 = MSR2/MSE2
fvalue2 = finv(0.95,2,n-3) % f-test value from table
p2 = fcdf(ftest2,2,n-3);
p2 = 1-p2 % Calculates the probability of
% a value as high as ftest2 resulting from a set of
% random inputs.
%
% Calculate the significance of adding the quadratic term 
% (the statistical significance of the improvement of the fit) 
SSR2m1 = SSR2 - SSR1 % SS due to regression added by quadratic fit
MSR2m1 = SSR2m1/1
% ftest to see if adding the quadratic term is a significant
% improvement over the straight line fit
ftest3 = MSR2m1/MSE2
fvalue3 = finv(0.95,1,n-3) % f-test value from tables
p3 = fcdf(ftest3,1,n-3);
p3 = 1-p3 % Calculates the probability of
% a value as high as ftest3 resulting from a set of
% random inputs.
%
% plot data and fits
plot(x,y,'or','markersize',12)
hold on
set(gca,'fontsize',16,'linewidth',2)
hold on
xlabel('x-values','fontsize',16)
ylabel('y-values','fontsize',16)
title('Plot of louismud data and LSQ fits','fontsize',16)
plot(x,yt1,'-k','linewidth',2)
plot(x(1:8),yt2(1:8),'-b','linewidth',2)
hold off