function [beta,v,prob,haty,e,varind]=mregauto(y,x,psig)
% fits a multiple linear regression model to response data y 
%  using predictor matrix x 
%
%  Uses backward selection to drop out predictors at 100.psig% significance level
%
%  ouputs are the LSE beta, the SDs of elements of beta in V, 
%  approximate probability level for each coeff, fitted values haty and residuals e
%  varind is the index vector of which variables in x are included
%
%  y is nx1 column vector, x is pxn matrix 
%  plots betas and approximate 95% intervals


[p,n]=size(x); y=reshape(y,n,1); ind=1:p; xw=x; 

u=1; 
for j=1:p

% model fit ...
X=xw; H=[ones(n,1),X']; V=inv(H'*H); beta=V*H'*y; 		
haty=H*beta; e=y-haty; s=sum(e'*e)/(n-length(beta)); v=sqrt(s*diag(V));              	
df=n-length(beta); prob=2*(1-pt(abs(beta./v),df));             
[a,i]=sort(prob(2:(p+1))); 
u=max(a); 
if (u>psig) 
   ind(i(p))=[]; p=p-1; xw=x(ind,:); 
else
   break
end;
end;

varind=ind; qtc=qt(0.975,df); 

xw=x(ind,:); 
X=xw; H=[ones(n,1),X']; V=inv(H'*H); beta=V*H'*y; 		
haty=H*beta; e=y-haty; s=sum(e'*e)/(n-length(beta)); v=sqrt(s*diag(V));              	
df=n-length(beta); prob=2*(1-pt(abs(beta./v),df));             


figure(4) 
plot([0 1+p],[0 0]); hold on
for j=1:p
   plot([j,j],[beta(j+1)+qtc*v(j+1),beta(j+1)-qtc*v(j+1)],'k')
   text(j,beta(j+1),'+','HorizontalAlignment','Center')
end;
xlabel('Predictor gene index'); ylabel('Coeff')
title('Regression coefficients') 
hold off;