/* Thanks to Joy Smith, North Carolina State University, 
   for help with the coding */
%LET seed=123;

DATA one; 
   INFILE "pvalues.txt" END=last; 
   INPUT raw_p;
   IF last THEN CALL SYMPUT("m",_n_);
RUN;

%MACRO picomp(dataset,criticalt);
%DO i=0 %TO 18;

DATA tmp;
   SET &dataset;
   lambda=.05*&i;
   *lambda=.8+.01*&i;  *for another choice of tuning parameter values;
   pgtlambda=(raw_p>lambda)/(1-lambda);
   plesst=(raw_p<&criticalt);
RUN;
PROC SORT DATA=TMP;  BY isample; RUN;
PROC MEANS DATA=tmp NOPRINT; 
   BY isample;
   ID lambda;
   VAR pgtlambda plesst;
   OUTPUT OUT=tout&i MEAN=pi0hat junk1 SUM=junk2 sum_plesst;
RUN;
PROC APPEND BASE=tout DATA=tout&i;
RUN;
%END;

PROC SORT DATA=tout;  BY isample; RUN;

ODS LISTING CLOSE;
PROC TPSPLINE;
   BY isample;
   MODEL pi0hat=(lambda)/DF=3;
   OUTPUT OUT=splot PRED;
RUN;
DATA tout;  SET tout;  IF lambda=0.9;RUN;
DATA pi0s;
   SET splot(RENAME=(p_pi0hat=pi0smooth)); 
   IF lambda=0.9;
   KEEP isample pi0smooth;
RUN;

DATA both; /* merge bootstrap samples and pi0hats */
   MERGE tout pi0s;
   BY isample;
   fdrhat = pi0smooth*&m*&criticalt/sum_plesst;
RUN;
%MEND;

%MACRO bsample(odata,boot,b);
DATA &boot;
   DO isample=1 TO &b;
   DO i= 1 TO &m;
      n=FLOOR(&m*RANUNI(&seed)+1);
      SET &odata POINT=n;
      OUTPUT;
   END;
   END;
   STOP;
RUN;
%MEND bsample;

%bsample(one,b1,500);
%picomp(b1,0.01);
ODS LISTING;
PROC MEANS DATA=both MEAN STD; 
   VAR fdrhat pi0smooth;
RUN;
PROC UNIVARIATE DATA=both NOPRINT;
   VAR fdrhat pi0smooth;
   OUTPUT OUT=summarize MEDIAN=fdrmed pi0med PCTLPTs=2.5 97.5 PCTLPRE=qfdr qpi0; 
run;
PROC PRINT DATA=summarize;
    TITLE "summary";
RUN;