/*The SURVPOW  macro code appears on page 103 in the
book, "SAS(R) Survival Analysis Techniques for Medical
Research, Second Edition" by Alan B. Cantor, Ph.D.  */


%macro survpow(s1=  , s2=    , nsub=365, actime= ,futime=  ,
rate=      ,p=.5, loss1=0, loss2=0, w=1, siglevel=.05) ;

/* Find number of points in data set for group 1 and convert to       vectors */

data _null_; 
	set &s1;
	i=_n_;
	call symput('counta',left(i));
run;
data y; 
	set &s1; 
	retain sa1-sa&counta ta1-ta&counta ;
	array surv{*} sa1-sa&counta;
	array ttime{*} ta1-ta&counta;
	t=t*⊄
	all=1;
	i=_n_;
	surv{i}=s; 
	ttime{i}=t; 
	output;
	proc sort; 
	by all;
data y; 
	set y; 
	by all; 
	if last.all;
	keep  all ta1-ta&counta sa1-sa&counta;
/*  Find number of points in data set for group 2 and convert to       vector */
data _null_; 
	set &s2;
	i=_n_;
	call symput('countb', left(i));
run;
data yy; 
	set &s2; 
	retain sb1-sb&countb tb1-tb&countb;
	array surv{*} sb1-sb&countb;
	array ttime{*} tb1-tb&countb;
	t=t*⊄
	all=1;
	i=_n_;
	surv{i}=s; ttime{i}=t; 
	output;
proc sort; 
	by all;
data yy; 
	set yy; 
	by all; 
	if last.all;
	keep  all tb1-tb&countb sb1-sb&countb;
/*  Find hazards at each partition point  */
data z; 
	all=1;
	do t=0 to (&actime+&futime)*&nsub - 1 ; 
		output; 
		end;
proc sort; 
	by all;
data merged; 
	merge z y yy; 
	by all;
	if trim("&counta") = "1" then lam1=-log(sa1)/ta1;
	%do i=1 %to &counta -1 ;
        %let j = %eval(&i+1);
        if ta&i le t lt ta&j then lam1 =
        (sa&i-sa&j)/((sa&j-sa&i)*(t-ta&i)+sa&i*(ta&j-ta&i));
	%end;
	if trim("&counta") = "2" and ta2 = . then do;
        lambda = -log((sa1-sa2)/(1-sa2))/ta1;
        lam1 = lambda*(1-sa2)*exp(-lambda*t)/(sa2+(1-sa2)*exp(-lambda*t));
        end;
	if trim("&countb") = "1" then lam2=-log(sb1)/tb1;
	%do i=1 %to &countb -1 ;
		%let j = %eval(&i+1);
			if tb&i le t lt tb&j then lam2 =
		(sb&i-sb&j)/((sb&j-sb&i)*(t-tb&i)+sb&i*(tb&j-tb&i));
		%end;
		if trim("&countb") = "2" and tb2 = . then do;
        lambda = -log((sb1-sb2)/(1-sb2))/tb1;
        lam2 = lambda*(1-sb2)*exp(-lambda*t)/(sb2+(1-sb2)*exp(-lambda*t));
        end;

/*  Calculate ratio of hazards and number at risk at each partition point and 
		accumulate needed sums */

data next; 
	set merged;  
	by all;
	retain n1 n2 n;
	if _n_ = 1 then do;
        n1=&rate*&p*&actime;
        n2=&rate*(1-&p)*&actime;
        n=n1+n2;
	end;
	tau=&futime*⊄
	psi1=&loss1/⊄
	psi2=&loss2/⊄
	phi=n1/n2;
	theta=lam1/lam2;
	d1=lam1*n1;
	d2=lam2*n2;
	d=d1+d2;
	c1=psi1*n1;
	c2=psi2*n2;
	if t > tau then do;
        c1=c1+n1/(&actime*&nsub+tau-_n_+1);
        c2=c2+n2/(&actime*&nsub+tau-_n_+1);
        end;
	n1=n1-d1-c1;
	n2=n2-d2-c2;
	sum1+(d*&w*(phi*theta/(1+phi*theta) - phi/(1+phi)));
	sum2+d*&w**2*phi/(1+phi)**2;
	n=n1+n2;

/*  Calculate e and power */

	if last.all then do;
        e=sum1/sqrt(sum2);
        z=-probit(&siglevel/2);
        power = 1 - probnorm(z-e) +  probnorm(-z-e);
        ac_time=symget('actime');
        fu_time=symget('futime');
        ac_rate=symget('rate');
        n=ac_rate*ac_time;
        alpha = symget('siglevel');
        prop=symget('p');
        los_rat1=symget('loss1');
        los_rat2=symget('loss2');
        weights = symget('w');
        output;
        end;
	label ac_time='Accrual Time';
	label power='Power';
	label fu_time = 'Followup Time';
	label ac_rate = 'Accrual Rate';
	label n = 'N';
	label prop = 'Prop in Grp 1';
	label los_rat1 = 'Loss Rate 1';
	label los_rat2 = 'Loss Rate 2';
	label weights = 'Weights';

/*  Print results */

proc print l noobs;
	var ac_time fu_time ac_rate n alpha prop los_rat1 los_rat2 	weights power;

/*
data results; set results next; 
*/
	run;
%mend;