/*  The PHPOW Macro  */

%macro phpow(t = , tau = , alpha = .05, n = ., power = ., var = ,   delta = , s0 = );

/*This is a SAS macro for calculating power and sample size. It is based on Hsieh and Lavori (2000), Controlled Clinical Trials, 21: 552-560*/

proc iml;
file print;
title 'Power and Sample Size Results';

 /*The inner integrand is calculated*/

   start integrand(time) global(yv,delta,var);
      pi = 4*atan(1);
      part1 = 1/sqrt(2*pi*var)*exp(-yv**2/(2*var));
	 part2 = 1 - (&s0)**(delta**yv); 	
       p = part1*part2;  
	 return(p);
   finish;

 /*The inner integral is calculated*/

   start marginal(v) global(yv,t, tau );
   tt = t;
   ttau = tau;
   upper = tt + ttau;
      interval = ttau||upper;
	yv = v;
	call quad(pm,"integrand",interval); 
	return(pm);
   finish;

 /*The outer integral and probability of death are calculated*/

start outer(tt,ttau,  svar, ddelta) global( t, tau, var,                                                                                         delta );
      t = tt;
      tau = ttau;
	var = svar;
	delta = ddelta;
      interval= .M ||.P;
      call quad(per,"MARGINAL",interval); 
	prob_d = per/t;
      return(prob_d);
   finish;
   prob_d =  outer(&t, &tau,  &var, &delta);
   
 /* Power is calculated if it is missing */

   %if &power = . %then %do;
   deaths = prob_d*&n;
   zpower = sqrt(&var)*abs(log(&delta))*sqrt(deaths)- 
   probit(1 - &alpha/2);
   power_calculated = round(probnorm(zpower), .1);
   sample_size = &n;
   %end;

 /* Sample size is calculated if it is missing */

   %else %do;
   zpower = probit(&power);
   deaths = 
   (zpower+probit(1-&alpha/2))**2/(abs(log(&delta)))**2/&var;
   sample_size_calculated = round(deaths/prob_d, 1);
   power = &power;
   %end;

 /*Output results*/

     Alpha = α
     hazard_ratio = δ
     accrual_time = &t;
     followup_time = τ
     covariate_variance = &var;
     baseline_survival = "&s0";
     put 'Alpha   ='alpha;
     put;
     put 'Hazard Ratio   ='hazard_ratio;
     put;
     put 'Accrual Time   ='accrual_time;
     put;
     put 'Followup Time   ='followup_time;
     put;
     put 'Covariate Variance   ='covariate_variance;
     put;
     put 'Baseline Survival =' baseline_survival;
     put;  
   %if &power=. %then %do;
     put 'Sample Size   =' sample_size;
     put;
     put 'Power(Calculated)   ='power_calculated;
     put;
    %end;
   %else %do;
    put 'Power   ='power;
    put;
    put 'Sample Size (Calculated)   =' sample_size_calculated;
   %end;
   run;
 %mend;