clear all
close all

 disp('Linear Regression: Body Fat Example One Variable')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;

load  'fat.dat'
broz = fat(:,2);
abdomen  = fat(:,12);  
n = length(broz);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 2; %number of parameters (beta0, beta1)
% beta0 estimated by b0; beta1 estimated by b1
%sums of squares
% 'Abdomen measurement is "x" and the response brox is "y".
% Sums of Squares
SXX = sum( (abdomen - mean(abdomen)).^2 )
SYY = sum(  (broz - mean(broz)).^2 )
SXY = sum( (abdomen - mean(abdomen)).* (broz - mean(broz)) )

% estimators of coefficients beta1 and beta0
b1 = SXY/SXX
b0 = mean(broz) - b1 * mean(abdomen)
% predictions
brozhat = b0 + b1 * abdomen;
%residuals
res = broz - brozhat;
% ANOVA Identity
SST = sum(  (broz - mean(broz)).^2 ) %which is SYY
SSR = sum(  (brozhat - mean(broz)).^2 )
SSE = sum(  (broz - brozhat).^2 ) %which is sum(res.^2)
% forming F and test of adequacy of linear regression
MSR = SSR/(p - 1)
MSE = SSE/(n - p)  %should be sigma2hat
F = MSR/MSE
pvalue = 1-fcdf(F, p-1, n-p) %H_0: regression has beta1=0, no need for 
% linear fit
%  Other measures of goodness of fit
R2 = SSR/SST
R2adj = 1 - (n-1)/(n-p)* SSE/SST
s = sqrt(MSE)
%  Standard error of coefficient estimators
sb1 = s/sqrt(SXX)
sb0 = s * sqrt(1/n + (mean(broz))^2/SXX )

% are the coefficients equal to 0?
t1 = b1/sb1
pb1 = 2 * (1-tcdf(abs(t1),n-p) )
t0 = b0/sb0
pb1 = 2 * (1-tcdf(abs(t0),n-p) )

% predicting y for the new observation x, CI and PI
newx = 230; %abdomen = 230
brozpred = b0 + b1 * newx
sym = s * sqrt(1/n + (mean(abdomen) - newx)^2/SXX ) %s for y mean
syp = s * sqrt(1 + 1/n + (mean(abdomen) - newx)^2/SXX ) %s for y prediction
%intervals CI and PI
alpha = 0.05;
%mean response interval
lbym = brozpred - tinv(1-alpha/2, n-p) * sym;
rbym = brozpred + tinv(1-alpha/2, n-p) * sym;
% prediction interval
lbyp = brozpred - tinv(1-alpha/2, n-p) * syp;
rbyp = brozpred + tinv(1-alpha/2, n-p) * syp;
%print the intervals
[lbym rbym]
[lbyp rbyp]


figure(1)
% scatterplot of (x,y)'s and regression
plot(abdomen, broz, 'ro')
hold on
plot(abdomen, brozhat, 'k-')
hold off

figure(2)
%residuals against x
stem(abdomen, res) 

figure(3)
%residuals against order
stem((1:n), res) 

%%
% built in matlab's black-box m-file


vecones = ones(size(broz));
X = [vecones abdomen];
[bx,bintx,rx,rintx,statsx] = regress(broz, X);

bx
bintx
rx
rintx
statsx

figure(4)
errorbar((1:n),rx(1:n), rintx(1:n,1), rintx(1:n,2))
hold on
plot((1:n),rx(1:n),'ro','markersize',5,'MarkerFaceColor','g') 
plot((1:n),1.96 * sqrt(statsx(4))*ones(1,n), 'r:')
plot((1:n),-1.96 * sqrt(statsx(4))*ones(1,n), 'r:')

axis tight