****************************************************************;
*  Preposterior Analysis for Bernoulli noninferiority trial     ;
*                                                               ;
*  Instructions                                                 ;
*     1) Highlight and run the lines from %MACRO to %MEND       ;
*        This needs to be done only once per SAS session        ; 
*     2) The prior is beta(a,b)                                 ;
*        k = TypeI/TypeII cost                                  ;
*        c = cost per observation / Type II cost                ; 
*        p0 = upper boundary of H0 (inferiority)                ;
*        p1 = lower boundary of H2 (superiority)                ;
*     3) For further details see                                ;
*        Biostatistics, A Bayesian Introduction                 ;
*        G.G. Woodworth                                         ;
*        Section 13.4.1                                         ; 
****************************************************************;

/* This macro call produces Figure 13.9 and table 13.3 */
%RatePrePost(a=14,b=1,k=12,c=.0003,p0=.92,p1=.96, n=100 TO 300 BY 1);

/* Step 1. Highlight from %Macro to %MEND RatePrePost; and click the runner  */

/* Step 2.  Fill in your own values here                                     */
%RatePrePost(a= ,b= ,k=  ,c= ,p0=  ,p1=  , n=  TO   BY  );

/* Step 3. Highlight the command you created in step 1 and click the runner  */
/* Repeat steps 2 and 3 as many times as you like                            */

%MACRO RatePrePost(a,b,k,c,p0,p1,n);
%LET ni=1;
PROC DATASETS LIBRARY=WORK;
	DELETE OUTMAT;
RUN;
QUIT;

PROC IML;
	a=&a; b=&b;
	_lbeta=lgamma(a+b)-lgamma(a)-lgamma(b);
	PriM0 = cdf('beta',&p0,a,b);
	PriM1=1-cdf('beta',&p1,a,b);
	PEL=(&k*PriM0||PriM1)[><]; NIPEL=PEL;
	pbar=(&p0+&p1)/2; sig=sqrt(pbar*(1-pbar));
	d0=(&p1-&p0);
	eta = d0/sig; nom=16/eta**2;
* Preposterior;
	DO n = &n;
		x = (0:n)`;
		fx =exp(_lbeta+lgamma(n+1)-lgamma(x+1)-lgamma(n-x+1)
		    +lgamma(a+x)+lgamma(b+n-x)-lgamma(n+a+b));
	* Joint Distribution of Dhat, Sed, and Delta ;
		PRM0=0; PRM1=0; MaxTail=0;
		PostM0=cdf('beta',&p0,a+x,b+n-x);
		PostM1=1-cdf('beta',&p1,a+x,b+n-x);
		Rej =(&k*PostM0 < PostM1);
		_len=nrow(Rej); 
		CV=_LEN - min(loc((Rej//{1})));
		CR=MIN({101}||(100*(1-CV/n))); 
		MAXTAIL = (REJ#PostM0)[<>];
		PRM0=(Rej#PostM0#fx)[+]; 
		PRM1=(Rej#PostM1#fx)[+];
		BTI=PRM0/PriM0; 
		BAP=PRM1/PriM1;
		ThisPEL=&k*PRM0+(PriM1-PRM1);
		DPEL=PEL-ThisPEL; 
		PEL=ThisPEL;

		OUTMAT = OUTMAT//(N||PriM0||PriM1||PRM0||PRM1||BTI||BAP||MAXTAIL||PEL||DPEL||CV||CR);
	END;
	COLNAMES={"N" "PM0" "PM1" "PRM0" "PRM1" "BTI" "BAP" "MAXTAIL" "PEL" "DPEL" "CritValue" "CritRate"};

	CREATE OUTMAT FROM OUTMAT [COLNAME=COLNAMES];
		APPEND FROM OUTMAT;
	QUIT;
RUN; 

DATA OUTMAT;
	SET OUTMAT;
	COST=PEL+&c*N;
RUN;
SYMBOL1 VALUE=DOT COLOR=BLACK I=J;
SYMBOL2 VALUE=DOT COLOR=RED I=J;

PROC MEANS NOPRINT MIN MAX DATA=OUTMAT;
	BY CritValue;
	VAR N;
	OUTPUT OUT=CritValues
		MIN=Min_N
		MAX=Max_N;
RUN;
OPTIONS NOCENTER;
PROC PRINT NOOBS DATA=CritValues;
	TITLE2 "Critical Values (-1 = Never Reject)";
	VAR CritValue Min_N Max_N;
RUN;
	

PROC GPLOT DATA=OUTMAT;
	TITLE2 "Minimum expected cost by sample size";
	PLOT (cost)*N / OVERLAY;
PROC GPLOT DATA=OUTMAT;
	TITLE2 "Probability of H0 (inferiority) at the CV";
	PLOT (PRM0)*N / OVERLAY;
PROC GPLOT DATA=OUTMAT;
	TITLE2 "Probability of H2 (superiority) at the CV";
	PLOT (PRM1)*N / OVERLAY;
PROC GPLOT DATA=OUTMAT;
	TITLE2 "Minimum Success Rate to Reject H0 by sample size";
	PLOT (CritRate)*N / OVERLAY;
RUN; 
QUIT;
%MEND RatePrePost;