/*NOTE THIS CODE WILL CALCULATE THE ESTIMATED SURVIVAL TIME FOR THE SAMPLE RANDOMIZED CONTROLED TRIAL DATASET (QAS_DATA). THE DATASET CONTAINS 6 VARIABLES:
1) ID - UNIQUE STUDY IDS
2) TRTMT - TREATMENT ASSIGNMENT; TREATMENT A (TRTMT=0) AND TREATMENT B (TRTMT=1)
3) DEATHDAY - TIME (YEARS) FROM START OF FOLLOW-UP UNTIL EITHER DEATH OR CENSORING
4) DEATH - INDICATOR FOR DEATH AS THE OUTCOME (DEATH=1) OR CENSORING (DEATH-=0)
5) T_COMBO - TIME (YEARS) FROM START OF FOLLOW-UP UNTIL EITHER THE COMBINED OUTCOME OR CENSORING
6) COMBO - INDICATOR FOR THE COMBINED OUTCOME OF ALL-CAUSE HOSPITALIZATION OR DEATH (COMBO=1) OR CENSORING (COMBO=0)*/

/********************************************************PROCEDURE 1*********************************************************************************/
/*Using PHREG, calculate the survival estimates for the treatment groups for 'death' as the outcome of interest*/ 
	

proc phreg noprint data = sasuser.QAS_Data;

	   model deathday*death(0) = trtmt / 

       ties = efron;

       output out = model0 survival = s1; /*output the survival probability estimates (S1) into a new dataset (model0)*/

       run;

/*Sort the new dataset (model0) by treatment group (trtmt) and time until death or censoring (deathday)*/
proc sort data=model0;

	   	by trtmt deathday;

	   run;

/*You will notice in the dataset model0 that when the survival probabilities were outputted they only begin according to the first event. 
	   As all people were event free at the start of the study (probability of survival = 1), we need to add the starting survival probability
	   at time equal to zero into our dataset model0 to allow full integration of the survival curve*/ 

data first; /*make new dataset 'first' containing a survival probability (s1)=1 at time=0*/

			deathday=0;

			s1=1;

		run;

/*Create new dataset (trtmt_A) of our survival probabilities by combining dataset 'modle0' and 'first'. We now integrate the 
		survival curve to calculate the mean event free survival time until either death or censoring. 
		We first do this for the subjects who were receiving Treatment A. */

  data trtmt_A; 

			set first model0(where=(trtmt=0));

			t=dif(deathday);

			s=lag(s1);

			integral=t*s;

			if integral=. then delete;

		run;

		proc  means noprint data=trtmt_A sum;

			var integral;

			output out=sum1 sum(integral)=integral;

		run;



proc means data=sum1 sum;/*get the summation of the area. This is the integral value for the area under the all-cause mortality curve*/

var integral;

title 'Intergral of Treatment A Survival Curve for All-Cause Mortality';

run;

/*We now integrate the survival curve to calculate the mean event free survival time until either death or censoring for subjects 
	receiving Treatment B*/


  data Trtmt_B; 

			set first model0(where=(trtmt=1));

			t=dif(deathday);

			s=lag(s1);

			integral=t*s;

			if integral=. then delete;

		run;

		proc  means noprint data=Trtmt_B sum;

			var integral;

			output out=sum2 sum(integral)=integral;

		run;



proc means data=sum2 sum;/*get the summation of the area. This is the integral value for the area under the all-cause mortality curve*/

var integral;

title 'Intergral of Treatment B Survival Curve for Mortality';

run;

/********************************************************PROCEDURE 2*********************************************************************************/

/*********************************************************************************************/
/*Repeat the above procedures to calculate the survival estimates for the treatment groups for the combined endpoint of 'all-cause hospitalization or death'
	as the outcome of interest and output the survival probability estimates (S1) into a new dataset (model1)*/
/***********************************************************************************************/

 proc phreg noprint data = sasuser.QAS_Data;

	   model t_combo*combo(0) = trtmt / 

       ties = efron;

       output out = model1 survival = s1;

       run;

	   proc sort data=model1;

	   	by trtmt t_combo;

	   run;



data first;

			t_combo=0;

			s1=1;

		run;



  data Trtmt_A; 

			set first model1(where=(trtmt=0));

			t=dif(t_combo);

			s=lag(s1);

			integral=t*s;

			if integral=. then delete;

		run;

		proc  means noprint data=Trtmt_A sum;

			var integral;

			output out=sum2 sum(integral)=integral;

		run;



proc means data=sum2 sum;

var integral;

title 'Intergral of Treatment A Survival Curve for Combined Outcome';

run;

/************************RUN FOR Treatment B for the Combined Endpoint**************************************/



  data Trtmt_B; 

			set first model1(where=(trtmt=1));

			t=dif(t_combo);

			s=lag(s1);

			integral=t*s;

			if integral=. then delete;

		run;

		proc  means noprint data=Trtmt_B sum;

			var integral;

			output out=sum3 sum(integral)=integral;

		run;



proc means data=sum3 sum;

var integral;

title 'Intergral of Treatment B Survival Curve for Combined Outcome';

run;

/********************************************************PROCEDURE 3*********************************************************************************/

/*To calculate the time spent in health state H1 for the treatment groups it is simply the 
	integrated time for the combined outcome for treatment A and B, respectively.*/
/*To calculate the time spent in health state H2 simply subtract the integrated time for the combined 
	endpoint from the integrated time for the death outcome*/ 
/*As observed in this sample dataset, despite a significant reduction in the risk of the combined endpoint for treatment B 
	compared to treatment A (HR 0.81, 95%CI 0.74-0.88, )there is what many would consider a minimal benefit of treatment B 
	over treatment A with respect to improvement in quality adjusted survival (0.06 years)QAS: 
		Treatment A = [(3.35 years * utility of 0.81) + ((4.38-3.35)*0.57)] = 3.30 QAS Years
		Treatment B = [(3.59 years * utility of 0.81) + ((4.39-3.59)*utility of 0.57))] = 3.36 QAS Years
		A difference of 0.06 QAS for Treatment B over Treatment A