% read in data, remove Affy controls and genes that do not vary across arrays 

addpath('./Utilities')  
load BC49DATA

x=AD7129;  x(x<=0)=1; x=log2(x); Des=Descriptions; 
           x([1430 117:174],:)=[]; Des([1430 117:174],:)=[]; 
           i=(std(x')<=0)'; x(i,:)=[]; Des(i,:)=[]; 

des=Des; X=x; [N,n]=size(X); X(X<0)=0; 

% set up validation cases 
ivalid=sort(unique([ivalidER 7 8]));



% ---------------------------------------------------------------------
% fit straight binreg MLE in matlab (R can do same, and more) 
%
X=std_rows(x); ind=select_genes_cor(X,C_ER,10); 
Des(ind,:)
[A,D,F]=svd_mw(X(ind,:)); 
[b,dev,stats]=glmfit(F',[C_ER' ones(n,1)],'binomial','probit'); 
stats






% ---------------------------------------------------------------------
% fit approximate Bayesian model, producing approx posterior means and SDs
%  of regression parameters 

[beta,li,alli]=mbinregsvd(X,C_ER,ivalid,100,5,200,1); 


% then look at top 100 genes ordered by abs(beta) ...
[b,ind]=sort(abs(beta)); ind=flipud(ind);
% beta(ind(1:10))
des(li(ind(1:25)),1:75)

% ---------------------------------------------------------------------
% fit approximate Bayesian model, producing MCMC simulation summaries
% of inferences 

[beta,li,alli]=Mbinregsvd(X,C_ER,ivalid,100,5,[2000 1000 1],0); 

% then look at top 100 genes ordered by abs(beta) ...
[b,ind]=sort(abs(beta)); ind=flipud(ind); % beta(ind(1:10))
des(li(ind(1:10)),:)