/*
 * DRAFT
 * Title CSSP Method - Analysis Macro
 * Programmer: Tracey Marsh
 *
 * Description
 * 		Performs a sequential test on an analytic data set - returns result (signal or not) and summary report.
 *      Interim files are stored and required for subsequent looks.
 *
 * Usage
 *		%CSSP_test(look_summary,this_look,info_prop, alpha, bnd_type="Poc", interimLib=work, detail=N);
 * 
 * Arguments
 *		look_summary = sas data set - analytic data summary of data new since previous look 
 *            contains variables as described:
 *					look, numeric - positive integers of the form 1-n: defines which look the data is from
 *					strata, numeric - positive integers of the form 1-n: defines analysis level subpopulations
 *					exp_doi, numeric - positive: total drug of interest exposure since previous look
 *				    exp_cmp, numeric - positive, >= exp_doi: total exposure (includes drug of interest) since previous look
 *					evt_doi, numeric - non-negative: number of adverse events with drug of interest exposure
 *                  evt_cmp, numeric - non-negative, >= evt_doi: number of adverse events (w or w/o exposure to drug of interest)
 *		this_look = an integer describing the look number - used for labelling and identifying segment of "new" data within look_summary table
 *		info_prop = proportion of information cumulatively inclusive of data new to this look. 
 *					Should be in (0,1], increasing with each look, equal to 1 at final look.
 *					Used in determination of boundary for given look.
 *		alpha = total amount of type I error cumulative over all looks (i.e. when info_prop = 1).
 *		bnd_type = type of error spending function.
 *				   Supported types: "Poc", "OBF", "Uni" ("Uni" assumed as default).
 * 		interimLib = SAS library for storing interim data sets and accessing interim data sets from previous looks.
 *		detail = Y to produce a more detailed summary report (includes strata level information).
 *
 * Details
 *		The distribution of adverse events for the drug of interest exposure group is simulated under the null 
 *		hypothesis (no risk difference between exposure to drug of interest and exposure to the comparator) as a binomial  
 *		event of trials equal to the total number of adverse events observed, with probability equal to the proportional exposure to
 *   	the drug of interest. The simulation is performed for each strata independently, and the total events aggregated. Events are also 
 *		simulated at each look independently, and a cumulative total maintained. (Hence the need to store and access interim tables.) A 
 *		cumulative p-value is inferred from the rank of the number of adverse events oserved with exposure to the drug of interest within 
 *		the simulated total events, conditional on having not previously signalled. If the p-value for observed events is less than the 
 *		cumulative error spending function, a signal is indicated.
 *
 * Output
 * 		Summary 
 *
 * References
 *	Li, L. A conditional sequential sampling procedure for drug safety surveillance. Statistics in Medicine 2009; 28:3124-3138. 
 *
 * Example
 *	%CSSP_test(look_summary=sasout.example, this_look=4, info_prop=4/8, alpha=0.05, bnd_type="Poc", interimLib=sasout, detail=Y);
 *  see accompanying files "CSSP Method - Example.sas" and "example_data.csv"
 */


%macro CSSP_test(look_summary=, this_look=, info_prop=, alpha=0.05, bnd_type="Poc", interimLib=work, detail=N);

%let nboot = 10000;

***DEFINE BOUNDARY***;
data _null_; 
	select (UPCASE(&bnd_type.));
		when ("OBF") alpha_spend = 2*( 1- cdf('NORMAL', Quantile('NORMAL',1-&alpha./2)/(&info_prop.)**(0.5)));
		when ("POC") alpha_spend = &alpha.*log(1 + (exp(1)-1)*(&info_prop.));
		when ("UNI") alpha_spend = &alpha.*(&info_prop.);
		otherwise alpha_spend = &alpha.*(&info_prop.);
	end; 
	boundary_rank = &nboot.*(1-alpha_spend);
	call symput("cutrank", boundary_rank); 
run;

%put look = &this_look. information proportion = &info_prop. cutrank = &cutrank.;


***CREATE WIDE DATA SET WITH STRATA STATS***;

*determine number of strata in dataset;
proc sort data = &look_summary. (keep=strata look where = (look=&this_look.))
          out = strata_unq nodupkey;
  by strata;
run;

%let dsnid = %sysfunc(open(strata_unq));
%let nstrata = %sysfunc(attrn(&dsnid.,NOBS));
%let rc = %sysfunc(close(&dsnid.));
%put nstrata  &nstrata.;
	
data strata_stats (keep = look strata e_sum p e_d);
	set &look_summary. (where = (look=&this_look.));
	length esum p1 8.;

	e_sum = evt_doi + evt_cmp;
	p = exp_doi/(exp_doi + exp_cmp);
	e_d = evt_doi;
run;	

proc transpose data=strata_stats out=strata_stats_t;
	by look;
run;

data tmp1; set strata_stats_t; if _NAME_='e_sum'; rename col1-col&nstrata. = e_sum1-e_sum&nstrata.; run;
data tmp2; set strata_stats_t; if _NAME_='p'; rename col1-col&nstrata. = p1-p&nstrata.; run;
data tmp3; set strata_stats_t; if _NAME_='e_d'; rename col1-col&nstrata. = e_d1-e_d&nstrata.; run;

data strata_stats (drop = _NAME_);
	merge tmp1 tmp2 tmp3;
	by look;
run;


***SIMULATE DISTRIBUTION OF NEW EVENTS UNDER NULL***;
data dist (keep = look iter events_new) ;
	set strata_stats;
	look = &this_look.;
	array a_esum {&nstrata.} e_sum1-e_sum&nstrata.;
	array a_p {&nstrata.} p1-p&nstrata.;
	array a_ed {&nstrata.} e_d1-e_d&nstrata.;

	*cumulative counts - summed over strata;
	events_new = 0;

	**iter1  holds data counts;
	iter = 1; 
	do strata = 1 to &nstrata.;
		if a_esum{strata} >. & a_p{strata} >. & a_ed{strata} >. then 
			events_new = events_new + a_ed{strata};
	end;
	output;
	**iter2 on holds simulation of expected counts**;
	do iter = 2 to &nboot.;
		events_new = 0;
		do strata = 1 to &nstrata.;
			if a_esum{strata} >. & a_p{strata} >. & a_ed{strata} >. then do;
				events = a_esum{strata};
				prob = a_p{strata};
				if prob = 1 then events_new = events_new + events;
				else if prob > 0 & events >0 then events_new = events_new + ranbin(0, events, prob);
			end;
		end;
		output;
	end;	
run;

***CUMULATIVE SIMULATED DISTRIBUTION***;
%if &this_look. = 1 %then %do;
	data dist_lim; set dist; rename events_new = evt_doi_cum; run;
%end;
%else %do;
	data dist_lim (keep = iter look evt_doi_cum); 
		merge &interimLib..distsim (in = a) dist;
		by iter;
		if a;
		evt_doi_cum = evt_doi_cum + events_new; *evt_doi_cum is the cumulative number of events up until time t;
	run;
%end;

***APPLY STOPPING ALGORITHM**;
proc rank data = dist_lim out = dist_w_rank ties = low;
	var evt_doi_cum;
	ranks eventrank;
run;

data &interimLib..distsim (keep = iter look evt_doi_cum) rej(keep = iter look signal) condp (keep = iter signal look condp eventrank cutrank evt_doi_cum);
	set dist_w_rank;
	retain signal 0;
	cutrank = &cutrank;

	attrib _all_ label = '';
	
	*first iteration within a sim is the observed data counts;
	if iter=1 then do;
		if eventrank > cutrank then do;
			signal = 1;
			put "SIGNAL - SIGNAL - SIGNAL";
			output rej;	
		end;
		else signal = 0;
		condp = 1-(eventrank-1)/&nboot.;
			*the estimated probability of observing "more or equal extreme" outcomes;
			*this is conditional because only events less than boundary stay in set - see next if clause;
		output condp;
	end;
	*once it signals, nothing will be output - this prevents continued running;
	if signal = 0 & eventrank <= cutrank then do;
		output &interimLib..distsim;
	end;
run;

***ALPHA ACTUALLY SPENT***;
%let dsnid = %sysfunc(open(&interimLib..distsim));
%let dist_obs = %sysfunc(attrn(&dsnid.,NOBS));
%let rc = %sysfunc(close(&dsnid.));

data condp;
	set condp;
	if (signal) then escp = .;
	else escp = 1-&dist_obs./&nboot.;
run;

***STORE INTERIM TABLES***;	
%if &this_look.=1 %then %do;
	data &interimLib..condp; 
		retain look evt_doi_cum eventrank cutrank condp escp signal;
		set condp (keep = look evt_doi_cum eventrank cutrank condp escp signal); run;
	data &interimLib..analysis_stats; set strata_stats; run;
	data &interimLib..analysis_raw; set &look_summary. (where = (look=&this_look.)); run;
%end; 
%else %do;
	proc append
		base = &interimLib..condp
		data = condp (keep = look evt_doi_cum eventrank cutrank condp escp signal);
	quit;
	proc append
		base = &interimLib..analysis_stats
		data = strata_stats;
	quit;
	proc append
		base = &interimLib..analysis_raw
		data = &look_summary. (where = (look=&this_look.));
	quit;
%end;
proc sort data =  &interimLib..condp
          out = &interimLib..condp;
  by descending look;
run;

***PRINT SUMMARY REPORT***;
ods noresults;
ods listing close;
data _null_;
	set sashelp.vslib;
	if libname="%upcase(&interimLib.)" then do;
		call symput("outPath",trim(path));
	end;
run;
ods rtf file = "&outPath.\CSSP Sequential Testing - Look &this_look..rtf";
footnote "CSSP Method - Sequential Testing.sas (&sysdate, &systime)";

data _NULL_;
	Length boundary $3.;
	boundary = &bnd_type.;
	call symput("bound",boundary);
run;

title "Method: CSSP - Look: &this_look. - Information Proportion: &info_prop - Boundary: &bound. - Alpha: &alpha."; 

*high level testing overview;
data condp_sum (keep = look evt_doi_cum p_value bndry_value signal);
	set &interimLib..condp (rename = (condp=p_value));
	bndry_value = 1-cutrank/&nboot.;
run;
data condp_sum;
	retain look evt_doi_cum p_value bndry_value signal;
	set condp_sum;
run;
proc print data=condp_sum;
	title2 "Testing outcomes for all looks through current look (&this_look.)";
run;

*data counts documentation - used in analysis;
proc summary data= &interimLib..analysis_raw nway missing;
  class look stdt eddt;
  var subj: exp: evt:;
  format subj: exp: evt: comma12.0;
  output out=sum1( drop=_type_ _freq_) sum=;
run;
proc sort data =  sum1
          out = sum1d;
  by DESCENDING look;
run;
proc print data=sum1d;
	title2 "Analysis for all looks through current look (&this_look.) based on aggregate counts:";
run;

*data counts documentation - provided this look;
proc summary data= &look_summary. nway missing;
  class look stdt eddt;
  var subj: exp: evt:;
  format subj: exp: evt: comma12.0;
  output out=sum2( drop=_type_ _freq_) sum=;
run;
proc sort data =  sum2
          out = sum2d;
  by DESCENDING look;
run;
proc print data=sum2d;
	title2 "Aggregate counts provided at current look (&this_look.).";
run;

data diff (drop=subj_c subj_d evt_c evt_d exp_c exp_d);
	merge sum1 (IN=ana) sum2 (IN=new rename = (subj_cmp=subj_c subj_doi=subj_d exp_cmp=exp_c exp_doi=exp_d evt_cmp=evt_c evt_doi=evt_d));
	by look;
	subj_cmp=-1*(subj_cmp-subj_c);
	subj_doi=-1*(subj_doi-subj_d);
	exp_cmp=-1*(exp_cmp-exp_c);
	exp_doi=-1*(exp_doi-exp_d);
	evt_cmp=-1*(evt_cmp-evt_c);
	evt_doi=-1*(evt_doi-evt_d);
run;

proc print data=diff;
	title2 "Discrepancies between analysis base and aggregate counts provided at cuurent look (&this_look.):";
run;
*optional summaries by strata;
%IF &detail.=Y %THEN %DO;
proc sort data =  &interimLib..analysis_stats
          out = &interimLib..analysis_stats;
  by descending look;
run;
proc print data=&interimLib..analysis_stats;
	title2 "Analysis for current look (&this_look.) based on following data statistics:";
	title3 "";
run;

proc sort data =  &interimLib..analysis_raw
          out = &interimLib..analysis_raw;
  by descending look strata;
run;
proc print data=&interimLib..analysis_raw;
	title2 "Analysis for all looks through current look (&this_look.) based on following counts:";
	title3 "";
run;
%END;
ods rtf close;
ods listing;
ods results;
%mend;