/************************************************************************/
/* SAS macro to compute sample sizes for QT studies 					*/
/* 																		*/
/* Authors:		Anne-Michelle Noone & George Luta 						*/
/* Date:		June 10, 2008											*/
/*																		*/
/* Requires: 	MVN_DIST.sas written by Genz and Bretz					*/
/* 																		*/
/* Inputs:																*/
/* 				SS	=	Variance estimate for subject					*/
/*						Default	=	12.30								*/
/*				SSD	=	Variance estimate for subject by day			*/
/*						Default	=	4.20								*/
/*				SST	=	Variance estmate for subject by time			*/
/*						Default	=	2.09								*/
/*				SSDT=	Variance estimate for subject by day by time	*/
/*						Default	=	3.13								*/
/*				SR	=	Variance estimate of the residual				*/
/*						Default	=	5.44								*/
/*				K	=	Number of replicates							*/
/*						Default	=	3									*/
/*  			TMPTS=	Number of post-dose time points							*/
/*						Default	=	5 									*/
/* 				SSLB=	Smallest sample size in power calculation		*/
/*				SSUB=	Maximum sample size in power calculation		*/
/*				TRTEFF=	List of mean differences between test drug and 	*/
/*						at each measurement								*/
/* 				ALPHA=	Significance level for power calculation 		*/
/*						Default	=	0.05								*/
/*				MVNFILE=Location of MVN_DIST.sas file (must include		*/
/*						path and file name)								*/
/*																		*/
/*	Outputs:	Prints to the screen the corresponding power under 		*/
/*				design A, B, and C for each sample size between SSLB 	*/
/*				and SSUB. 												*/
/* 				Plots power curves under each design.					*/
/*																		*/
/* Example:	Configuration 1												*/
/* %QTPower(K =	3,														*/
/*			TMPTS=	5,													*/
/*			SSLB=	30, 												*/
/*			SSUB=	100, 												*/
/*			TRTEFF=	(5 5 5 5 5), 										*/
/*			ALPHA=	0.05,												*/
/*			MVNFILE= B:\misc\QT\multnorm.sas);							*/											*/
/************************************************************************/;

%macro QTPower(	SS		=		12.30, 
				SSD		=		4.20, 
				SST		=		2.09,
				SSDT	=		3.13, 
				SR		=		5.44, 
				K		=		3,
				TMPTS	=		5,
				SSLB	=		%str(), 
				SSUB	=		%str(), 
				TRTEFF	=		%str(), 
				ALPHA	=		0.05,
				MVNFILE	=		%str());

	data one (keep=ss a1-a&&TMPTS b1-b&&TMPTS c1-c&&TMPTS d1-d&&TMPTS);
		retain ss a1-a&&TMPTS b1-b&&TMPTS c1-c&&TMPTS  d1-d&&TMPTS;

		/**** Variance Components ****/
		s2s=&SS**2; 
		s2sd=&SSD**2; 
		s2st=&SST**2; 
		s2sdt=&SSDT**2; 
		s2r=&SR**2; 

		/**** Variances ****/
		sa1=sqrt(2*(s2sd+s2sdt+s2r/&K));
		sa2=sqrt(4*(s2sdt+s2r/&K));
		sb1=sqrt(4*(s2sd+s2sdt+s2r/&K));
		sb2=sqrt(8*(s2sdt+s2r/&K));

		/**** Correlations ****/
		ra1=s2sd/(s2sd+s2sdt+s2r/&K);
			call symput('R', ra1); /* b1 has same correlation */
		ra2=0;
		rb2=0;

		/**** Assign Hypothesized Values ****/
		array DELTA(&TMPTS);
		%do i=1 %to &TMPTS;
			DELTA(&i)=%scan(&TRTEFF, &i);	
		%end;

		*** Upper Values ***;
		/* Initilaize arrays */
		array A(&TMPTS); 
		array B(&TMPTS); 
		array C(&TMPTS);
		array D(&TMPTS);
		%do ss=&SSLB %to &SSUB %by 1;
			%do i=1 %to &TMPTS;
				A(&i)=(sqrt(&ss))*(10-DELTA(&i))/sa1-probit(1-&ALPHA); /* Design A1 */
				B(&i)=(sqrt(&ss))*(10-DELTA(&i))/sa2-probit(1-&ALPHA); /* Design A2 */
				C(&i)=(sqrt(&ss))*(10-DELTA(&i))/sb1-probit(1-&ALPHA); /* Design B1 */
				D(&i)=(sqrt(&ss))*(10-DELTA(&i))/sb2-probit(1-&ALPHA); /* Design B2 */
			%end;
			ss=&ss;
			output;
		%end;
		run;

	/* File must be included each run to open IML session */
	%include "&MVNFILE";

		use one;
		read all into mat;

		xyz=j(&SSUB-&SSLB+1, 5, 0);
		N=&TMPTS;

		/* Retain previous power for use in design A2 */
		prev=0;

		/* Set up parameters for MVN_DIST */
		INFIN = J(1,N,0);
	  	MAXPTS = 2000*N*N*N; 
		ABSEPS = 0.0001; 
		RELEPS = 0;

		/* Set up covariance matrices for each design */
		/* Design A1 and B1 */
		COVAR1 = J(N, N, &R);
		COVAR1t=COVAR1;
		do i=1 to N;
			COVAR1t[i,i]=1;
		end;

		/* Design A2 and B2*/
		COVAR2 = I(N);

		do i=1 to &SSUB-&SSLB+1;
			xyz[i,1]=mat[i,1];

		 	/* Design A1 */
			LOWER = mat[i,2:N+1]; 
		  	UPPER = mat[i,2:N+1];

		  	RUN MVN_DIST(N, LOWER, UPPER, INFIN, COVAR1t, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM);
			xyz[i,2]=value;

		 	/* Design A2 */
			LOWER = mat[i,(N+2):(N*2+2)]; 
			UPPER = mat[i,(N+2):(N*2+2)];

			if (prev<0.99) then do; *If previous power is less than 0.99 then continue;
				RUN MVN_DIST(N, LOWER, UPPER, INFIN, COVAR2, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM);
				xyz[i,3]=value;
				prev=value;
			end;
			else do; *else assign power of 1;
				xyz[i,3]=1;
			end;

			/* Design B1 */ 
			LOWER = mat[i,(N*2+2):(N*3+1)]; 
			UPPER = mat[i,(N*2+2):(N*3+1)];
			INFIN = J(1,N,0);

		  	RUN MVN_DIST(N, LOWER, UPPER, INFIN, COVAR1t, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM);

			xyz[i,4]=value;

			/* Design B2 */
		  	LOWER = mat[i,(N*3+2):(N*4+1)]; 
			UPPER = mat[i,(N*3+2):(N*4+1)];
			INFIN = J(1,N,0);

		  	RUN MVN_DIST(N, LOWER, UPPER, INFIN, COVAR2, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM);

			xyz[i,5]=value;

		end;

		mattrib xyz
        colname=({SampSize A1 A2 B1 B2})
        label="Sample Size and Corresponding Power under each Design"
        format=3.2;
		print xyz;

	  	create out from xyz[colname={"Sample Size" "A1" "A2" "B1" "B2"}];
	 	append from xyz;

	QUIT;

	/* Plot Power Curves */
	data annote; 
		length function $8;

		x=&SSUB*0.89;	y=.25;	function='move'; 	xsys='2'; ysys='2'; 	output;
		x=&SSUB;		y=.25;	function='draw';	xsys='2'; ysys='2';  	output;
		x=&SSUB;		y=0.01;	function='draw'; 	xsys='2'; ysys='2';  	output;
		x=&SSUB*0.89;	y=0.01;	function='draw'; 	xsys='2'; ysys='2';  	output;
		x=&SSUB*0.89;	y=.25;	function='draw'; 	xsys='2'; ysys='2';  	output;
		x=&SSUB*0.91;	y=.2;	function='label'; 	xsys='2'; ysys='2'; 	text='A1'; 	output;
		x=&SSUB*0.91;	y=.15;	function='label'; 	xsys='2'; ysys='2';  	text='A2'; 	output;
		x=&SSUB*0.91;	y=.1;	function='label'; 	xsys='2'; ysys='2';  	text='B1'; 	output;
		x=&SSUB*0.91;	y=0.05;	function='label'; 	xsys='2'; ysys='2';  	text='B2'; 	output;
		x=&SSUB*0.99;	y=.2;	function='move'; 	xsys='2'; ysys='2';  	output;
		x=&SSUB*0.94;	y=.2;	function='draw'; 	xsys='2'; ysys='2';  	line=1; 	output;
		x=&SSUB*0.99;	y=.15;	function='move'; 	xsys='2'; ysys='2';   	output;
		x=&SSUB*0.94;	y=.15;	function='draw'; 	xsys='2'; ysys='2';   	line=2; 	output;
		x=&SSUB*0.99;	y=.1;	function='move'; 	xsys='2'; ysys='2';   	output;
		x=&SSUB*0.94;	y=.1;	function='draw'; 	xsys='2'; ysys='2';   	line=20; 	output;
		x=&SSUB*0.99;	y=0.05;	function='move'; 	xsys='2'; ysys='2';   	output;
		x=&SSUB*0.94;	y=0.05;	function='draw'; 	xsys='2'; ysys='2';   	line=40; 	output;

		axis1 minor=none label=(angle=90 "Power") order=(0 to 1 by 0.1);
		axis2 minor=none label=("Sample Size");

		symbol1 value=none color=black i=join line=1;
		symbol2 value=none color=black i=join line=2;
		symbol3 value=none color=black i=join line=20;
		symbol4 value=none color=black i=join line=40;

		goptions ftext=swiss;

		proc gplot data=out anno=annote;
		plot (a1 a2 b1 b2)*sample_size / overlay vaxis=axis1 haxis=axis2 ;
		title1 "Power curves for the three designs to define treatment differences"; 
		title2 "(True treatment differences at the &TMPTS post-baseline time points equal: &TRTEFF)"; 
		run;quit;title1;title2;

%mend QTPower;