/* 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;
proc nlmixed data=one;
   parameters pi0=0.5 a=1.5 b=1.5;
   pi1=1-pi0;
   like=log(pi0+pi1*pdf('BETA',raw_p,a,b));
   model raw_p ~ general(like);
   ods output parameterestimates=parmz;
run;

data plot5;  
   infile "pvalues.txt";
   set one;
   array counts count1-count20;
   array midpts midpt1-midpt20;
   input raw_p;
   count1=(raw_p<=.05);
   midpt1=.025;
   do i=2 to  20;
      binedge=.05*i;
      counts{i}=(raw_p<=binedge) * (raw_p>binedge-.05);
      output;
   end;
run;

proc means data=plot5 noprint nway;
   var count1-count20; 
   output out=hts mean=ht1-ht20; 
run;
data hts;  set hts;  keep ht1-ht20;run;
proc transpose data=hts out=thts prefix=height;run;
data thts;
   set thts;
   midpt=_n_*.05-.025;
   keep height1 midpt;
run;

proc transpose data=parmz out=ests;
   id Parameter;
   var Estimate;
run;
data plot5;
   if _n_=1 then set ests;
   set thts;
   height2=pi0+(1-pi0)*pdf('BETA',midpt,a,b);
   height1=height1*20;
run;

data two; 
   set ests;
   *do i=1 to 384;
   do i=1 to &m;
      u=ranuni(&seed);
      if u<pi0 then pvalue=ranuni(&seed);
      else do; u=ranuni(&seed); 
      pvalue=betainv(u,a,b); 
      end;
   output; 
   end;
   keep pvalue;
run;

proc sort data=two;  by pvalue;run;
proc sort data=one;  by raw_p;run;
data plot6;
   merge one two;
run;
data plot7;
   if _n_=1 then set ests;
   set one;
   pdata=pi0+(1-pi0)*pdf('BETA',raw_p,a,b);
   p0data=pi0/pdata;
   p1data=1-p0data;
   if raw_p>0.1 then delete;
   keep p0data p1data raw_p pi0 a b;
run;

axis1 label=(angle=90 'density');
axis2 order=0 to 1 by .1 offset=(5) label=('p-value');
symbol1 c=black i=needle v=none w=5;
symbol2 c=black i=j v=none;
proc gplot data=plot5;
   title "Pvalue histogram, Two-component mixture"; 
   plot height1*midpt height2*midpt/overlay vaxis=axis1 haxis=axis2 ;
run; 
quit;
symbol1 i=join value=dot;

axis1 label=(angle=90 'p-value quantiles');
axis2 label=('q-value quantiles');
symbol1 i=none value=dot;
symbol2 i=joine value=none;
proc gplot data=plot6;   
   title "QQplot";
   plot pvalue*raw_p=1 raw_p*raw_p=2/overlay vaxis=axis1 haxis=axis2;
run;
quit;
symbol1 i=join value=none;
axis1 label=(angle=90 'Posterior probability');
axis2 label=('Raw p-value');

proc gplot data=plot7;
   title "Posterior probability of non-differential expression";
   plot p0data*raw_p/vaxis=axis1 haxis=axis2;
run; quit;


proc greplay igout=gseg nofs;
tc=sashelp.templt; 
*template=H3s;
template=H3;
   treplay 1:gplot4
           2:gplot5
           3:gplot6;
run;
quit;