********************************************************************************
* Purpose:  Compare the efficiency between GLM model and Rank-ANCOVA model
*			For mixed distribution of one Normal population and one single point distribution of 0;
* 
* Macro variables: N_LOCF=# of subjects whose endpoint score is used in the analysis;
* 	 N1_BOCF=# of subjects in the active group whose Baseline score is used in the analysis;
*    N2_BOCF=# of subjects in the placeboe group whose Baseline score is used in the analysis;
*	 DELTA=Effect size;
*    STD=Standard deviation of change from baseline;
*    N_Simu= # of simulations to be performed;
*
* Last modificaton: 14-Jan-2010
*
********************************************************************************;


%macro mBOCF(N_LOCF=,N1_BOCF=,N2_BOCF=,N_site=,DELTA=,STD=,N_simu=);

proc printto print='j.temp' new;run;

 %do j=1 %to &N_Simu;
  ****** Generate random data set *****;
	Data LOCF(keep=trt painchg);
	 do i=1 to &N_LOCF;
	  painchg=&std*normal(0);
	  if mod(i,2)=0 then do;
	   trt='ACTIVETX';
	   painchg=painchg+Δ
	   output;
	  end;
	  else do;
	   trt='PLACEBO';
	   output;
	  end;
     end;
	 run;

	 data BOCF(keep=trt painchg);
	  painchg=0;
	  do i=1 to &N1_BOCF;
	   trt='ACTIVETX';
	   output;
	  end;
	  do i=1 to &N2_BOCF;
	   trt='PLACEBO';
	   output;
	  end;
     run;

	 data all;
      set LOCF BOCF; 
	  N= &N_LOCF+&N1_BOCF+&N2_BOCF;
	  do i=1 to N;
	  baseline=abs(6.6+1.9*normal(0));
	  if baseline >10 then baseline=10;
	   SITE=mod(int(N*ranuni(0)),&N_site); *** Assuming N_site sites ***;
	   if ranuni(0)<0.6 then rescue_med='NO';	*** Assuming 60% did not take rescue med. ***;
	   else rescue_med='YES';
	   if ranuni(0)<0.3 then D_drug='YES'; *** Assuming 30% used D-drug ***;
	   else D_drug='NO';
	  end;
     run;

	 ****************** Perform RANK ANCOVA *************;
	 *** Step1 Compute the ranks of the response variable (Y) and covariate (X);

 	Proc sort data=all;by site Rescue_med D_drug;run;

 	proc rank data=all  nplus1 ties=mean out=ranks(keep=trt baseline site painchg D_drug Rescue_med); 
  		by site Rescue_med D_drug;
  		var baseline painchg;
 	run;

 *** Step2 Calculate the residuals from the linear regression of ranks of Y on ranks of X;

 	proc reg data=ranks noprint;
  		by site Rescue_med D_drug;;
  		model painchg=baseline;
  		output out=residual r=resid;
 	run;

	ods output CMH=CMH&j;
 *** Step3 CMH mean score statistics is used to compare the mean values between groups;

 	proc freq data=residual;
 	 	tables site*rescue_med*D_drug*trt*resid/noprint cmh2;
 	run;

	ods output DIFFS=DIFFS&j;
 ***************** Perform ANCOVA ************************;
 	PROC MIXED DATA= all noclprint nobound noinfo noitprint noprofile;
		CLASS    trt rescue_med D_drug site;
		MODEL    painchg = trt baseline rescue_med D_drug site;
		lsmeans trt/cl pdiff;
 	run; 

	data cmh&j(keep=prob);set cmh&j;if statistic=1;run;
	data diffs&j(keep=probt estimate);set diffs&j;run;
	data cmh_glm&j;merge cmh&j diffs&j;run;
  %end;

	proc printto;run;
  ******* Concatenats all output datasets ******;

	data outdt(drop=prob probt);	format trial_CMH trial_glm $8.;
	  set %do j=1 %to &n_simu; cmh_glm&j %end;;
	  *if estimate<0 then delete;
	  rename estimate=delta prob=p_CMH probt=p_glm;
	  if prob<0.05 and estimate > 0 then trial_CMH='Positive';
	  else trial_CMH='Negative'; 
	  if probt<0.05 and estimate > 0 then trial_glm='Positive';
	  else trial_glm='Negative';
	run;

	title "Mean p-values from Rank-ANCOVA vs. Parametric GLM from &N_Simu Simulations";
	title2 "Data were from 372 numbers randomly generated from N(&delta,2.2^2) and 44 Zero's";
	title3 "Number of Sites=&N_site";
	proc means data=outdt;run;
	title 'Probability of Positive Trial using  Rank-ANCOVA vs. Parametric GLM from &N_Simu Simulations';
	title2 "Data were from 372 numbers randomly generated from N(&delta,2.2^2) and 44 Zero's";
	title3 "Number of Sites=&N_site";
	proc freq data=outdt;tables trial_CMH trial_glm; run;
%mend;

%mBOCF(N_LOCF=372,N1_BOCF=24,N2_BOCF=20,N_site=10,DELTA=0.8,STD=2.2,N_Simu=50)