/* Thanks to Joy Smith for help with the coding */
options ls=75  ps=1000;
%let seed=123;
data one; 
   infile "pvalues.txt"; 
   input raw_p;
run;
data one;
  set one end=last;
   if last then call symput("m",_n_);
run;

%macro pi0comp;
%do i=0 %to 18;

data tmp;
   set one;
   lambda=.05*&i;
   bigger=(raw_p>=lambda)/(1-lambda);
run;

proc means data=tmp noprint; 
   id lambda;
   var bigger ;
   output out=tout&i mean=pi0hat;
run;
proc append base=tout data=tout&i;
run;
%end;

%mend;

%pi0comp;

data tout;
   set tout;  
   keep lambda pi0hat;
run;
proc tpspline ;
   model pi0hat=(lambda)/df=3;
   output out=splot;
run;

data _null_;
   set splot end=last;
   if last then call symput("pi0",P_pi0hat);
run;

proc rank data=one out=rone ties=high;
   var raw_p;
   ranks v;
run;

data rone;
   set rone;
   preqvalue=&pi0*&m*raw_p/v;
   orderp=_n_;
run;

proc sort data=rone out=orderedbyp;
   by raw_p;
run;

proc iml;
print &pi0;
use orderedbyp;
read all var{orderp} into torderp;
orderp=(torderp)`;
use rone;
read all var{preqvalue} into qvalue; 
read all var{raw_p} into pvalue;

qvalue[orderp[&m]] = min(qvalue[orderp[&m]],1);
do i=(&m-1) to 1 by -1;
qvalue[orderp[i]] = min(qvalue[orderp[i]],qvalue[orderp[i+1]],1);
end;


pqmtx=pvalue||qvalue;

/*write pvalue qvalue to new dataset;*/

create pqvalues var{pvalue qvalue}; 
append; 
close pqvalues;
quit;

data pqr;
   set pqvalues;
   if pvalue<.05;
run;
proc print data=pqr;
run;

symbol1 i=join color=black value=none;

proc  sort data=pqr;
   by pvalue;
run;

proc sort data=pqr;
   by qvalue;
run;

data pqrg3;
   set pqr;
   sigtests=_n_;
   sigXq=sigtests*qvalue;
run;

options ls=75  ps=1000;

*goptions reset=all lfactor=5 hsize=4.5 in vsize=3.5 in;
goptions reset=all; 
*goptions dev=ps lfactor=5;

symbol1 i=none color=blue value=dot;
symbol2 i=join color=black value=dot;
symbol3 i=join color=black value=dot;
symbol4 i=join color=black value=dot;
symbol5 i=join color=blue value=none;

axis1 label=(angle=90 h=2 'density');
axis2 label=(font=greek h=2 'l');

proc gplot data=splot;
   title "Plot 1";
   plot pi0hat*lambda=1 p_pi0hat*lambda=5/overlay haxis=axis2 vaxis=axis1;
run;quit;

axis1 label=(angle=90 h=1.5 '# significant tests');
axis2 label=(h=2 'q-value cut-off');
proc gplot data=pqrg3;
   title "Plot 3";
   plot sigtests*qvalue=5/haxis=axis2 vaxis=axis1;
run;quit;

axis1 label=(angle=90 h=2 'q-value');
axis2 label=(h=2 'Raw p-value');

proc gplot data=pqr;
   title "Plot 2";
   plot qvalue*pvalue=5/haxis=axis2 vaxis=axis1;
run;quit;

axis1 label=(angle=90 h=2 'E(false positives)');
axis2 label=(h=2 '# of significant tests');
proc gplot data=pqrg3;
   title "Plot 4";
   plot sigXq*sigtests=5/haxis=axis2 vaxis=axis1;
run; quit;

proc greplay igout=gseg nofs;
tc=sashelp.templt; 
template=l2r2s;
   treplay 1:gplot
           2:gplot1
           3:gplot2
           4:gplot3;
run;
quit;