a code for stepwise regression with twelve possible fitting parameters

% litter.m


file = importdata('biomass_workbook2.csv',',');


% names of information | bio in 1980 | percenthardwood in 1980 |
% hardwood bio 1980|...etc.
bio80 = file(:,1);
perhard80 = file(:,2);
hardbio07 = file(:,3);
perhard07 = file (:,4);
litterdat = file(:,5);
anpp1 = file(:,7);
basal07 = file(:,8);
anpp6 = file(:,9);
bio07 = file(:,11);
num_tree = file(:, 12);
stemden = file(:,13);
logherb = file(:,14);
hwbio80 = file (:,15);


% poly fit finds the parmaters (pi) and the norms of the residuals (si) for
% linear and poly combinations of variables


[p1,s1] = polyfit(bio80, litterdat, 1);
[p2,s2] = polyfit(perhard80, litterdat, 1);
[p3,s3] = polyfit(hardbio07, litterdat, 1);
[p4,s4] = polyfit(perhard07, litterdat, 1);
[p5,s5] = polyfit(anpp1, litterdat, 1);
[p6,s6] = polyfit(basal07, litterdat, 1);
[p7,s7] = polyfit(anpp6, litterdat, 1);
[p8,s8] = polyfit(bio07, litterdat, 1);
[p9,s9] = polyfit(num_tree, litterdat, 1);
[p10,s10] = polyfit(stemden, litterdat, 1);
[p11,s11] = polyfit(logherb, litterdat, 1);
[p12,s12] = polyfit(hwbio80, litterdat, 1);


% polyval finds the correlation matrix
output1 = polyval(p1, bio80);
corr1 = corrcoef(litterdat,output1);


output2 = polyval(p2, perhard80);
corr2 = corrcoef(litterdat,output2);


output3 = polyval(p3, hardbio07);
corr3 = corrcoef(litterdat,output3);


output4 = polyval(p4, perhard07);
corr4 = corrcoef(litterdat,output4);


output5 = polyval(p5, anpp1);
corr5 = corrcoef(litterdat,output5);


output6 = polyval(p6, basal07);
corr6 = corrcoef(litterdat,output6);


output7 = polyval(p7, anpp6);
corr7 = corrcoef(litterdat,output7);


output8 = polyval(p8, bio07);
corr8 = corrcoef(litterdat,output8);


output9 = polyval(p9, num_tree);
corr9 = corrcoef(litterdat,output9);


output10 = polyval(p10, stemden);
corr10 = corrcoef(litterdat,output10);


output11 = polyval(p11, logherb);
corr11 = corrcoef(litterdat,output11);


output12 = polyval(p2, hwbio80);
corr12 = corrcoef(litterdat,output12);


% M is a matrix to hold thePearson's correlation coefficient
m = zeros(12,1);
m(1) = corr1(1,2);
m(2) = corr2 (1,2);
m(3) = corr3 (1,2);
m(4) = corr4 (1,2);
m(5) = corr5 (1,2);
m(6) = corr6 (1,2);
m(7) = corr7 (1,2);
m(8) = corr8 (1,2);
m(9) = corr9 (1,2);
m(10) = corr10 (1,2);
m(11) = corr11 (1,2);
m(12) = corr12 (1,2);


% yh(1) holds the fit of the predicted 


yh1 = polyval(p1,bio80);
yh2 = polyval(p2, perhard80);
yh3 = polyval(p3, hardbio07);
yh4 = polyval(p4, perhard07);
yh5 = polyval(p5, anpp1);
yh6 = polyval(p6, basal07);
yh7 = polyval(p7, anpp6);
yh8 = polyval(p8, bio07);
yh9 = polyval(p9, num_tree);
yh10 = polyval(p10, stemden);
yh11 = polyval(p11, logherb);
yh12 = polyval(p12, hwbio80);


% resid is what is not predicted


resid1 = litterdat - yh1;
resid2 = litterdat - yh2;
resid3 = litterdat - yh3;
resid4 = litterdat - yh4;
resid5 = litterdat - yh5;
resid6 = litterdat - yh6;
resid7 = litterdat -yh7;
resid8 = litterdat - yh8;
resid9 = litterdat - yh9;
resid10 = litterdat - yh10;
resid11 = litterdat - yh11;
resid12 = litterdat - yh12;


rss1 = sum((yh1- mean(litterdat)).^2);
tss1 = sum((litterdat - mean(litterdat)).^2);


% Rsq asks how much better predictedis thanmean
Rsq1 = 1-sum(resid1.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq2 = 1-sum(resid2.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq3 = 1-sum(resid3.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq4 = 1-sum(resid4.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq5 = 1-sum(resid5.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq6 = 1-sum(resid6.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq7 = 1-sum(resid7.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq8 = 1-sum(resid8.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq9 = 1-sum(resid9.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq10 = 1-sum(resid10.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq11 = 1-sum(resid11.^2)/sum((litterdat-mean(litterdat)).^2);
Rsq12 = 1-sum(resid12.^2)/sum((litterdat-mean(litterdat)).^2);


% n is a matrix of R sq
n= zeros(12,1);
n(1) = Rsq1;
n(2) = Rsq2;
n(3) = Rsq3;
n(4) = Rsq4;
n(5) = Rsq5;
n(6) = Rsq6;
n(7) = Rsq7;
n(8) = Rsq8;
n(9) = Rsq9;
n(10) = Rsq10;
n(11) = Rsq11;
n(12) = Rsq12;


% mat is a matrix of predictors
mat = zeros(length(litterdat),12);
mat(:,1:4) = file(:,1:4);
mat(:,5:7) = file(:,7:9);
mat(:,8:12) = file(:,11:15);


% this will do a stepwise fit. b isparameters, se is residuals, pval is
% values.. in model shows the process of gettingthere, statsis stats,next
% step and history are more dtailed
[b, se, pval, inmodel, stats, nextstep, history] = stepwisefit(mat, litterdat);


% model_fit is parameters
model_fit = inmodel.*b';
model_fit = model_fit';


% model 2is the predicted
model2 = zeros(length(mat),1);
model2= 4.4443 + model_fit(4).*mat(:,4) + model_fit(7)*mat(:,7) + model_fit(8).*mat(:,8);


% calculateresiduals
resid_step = litterdat-model2;


% calculate r2
Rsq_step = 1-sum(resid_step.^2)/sum((litterdat-mean(litterdat)).^2);


% calculate a fitting "line" to plot
[p_step, s_step] = polyfit(model2, litterdat,1);
a = p_step(1);
b = p_step(2);
Vcalc = a.*model2 + b;


% plot
figure(1)
plot(model2,litterdat,'k.');
hold on
plot(model2, Vcalc,'r-');
xlabel('Stepwise litterfall model')
ylabel('Litterfall data')
title ('Fitness testing for litterfall model');