********************** test.sas  ******************;
***  Program to test Fleming-Harrington Macro   ***;
***  Programmer: Andy Liljestrand  				***;
***			September 2003						***;
***************************************************;


options ls=121 ps=65 macrogen symbolgen mprint missing=' ' nonumber nodate;
/*ods trace on;*/
dm log 'clear'; 
dm output 'clear';
/*
footnote "&botline";
*/
%let debug=ON;

data example;
	input treat months censor;
	cards;

	1 1.5 0
	1 3.5 0
    1 4.5 0
	1 4.5 0
	1 5.5 0
	1 8.5 0
	1 8.5 0
	1 9.5 0
	1 10.5 0
    1 11.5 0
	1 15.5 0
	1 16.5 0
    1 18.5 0
	1 23.5 0
	1 26.5 0
	1 2.5 1
	1 2.5 1
	1 3.5 1
	1 3.5 1
	1 3.5 1
	1 4.5 1
	1 5.5 1
	1 6.5 1
	1 6.5 1
	1 7.5 1
	1 7.5 1
	1 7.5 1
	1 7.5 1
	1 8.5 1
	1 9.5 1
	1 10.5 1
	1 11.5 1
	1 12.5 1
	1 12.5 1
	1 13.5 1
	1 14.5 1
	1 14.5 1
	1 21.5 1
	1 21.5 1
	1 22.5 1
	1 22.5 1
	1 25.5 1
	1 27.5 1
	2 0.5 0 
	2 0.5 0
	2 0.5 0
	2 0.5 0
	2 0.5 0
	2 0.5 0
	2 2.5 0
	2 2.5 0
	2 3.5 0
	2 6.5 0
	2 15.5 0
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 0.5 1
	2 1.5 1
	2 1.5 1
	2 1.5 1
	2 1.5 1
	2 2.5 1
	2 2.5 1
	2 2.5 1
	2 2.5 1
	2 2.5 1
	2 3.5 1
	2 3.5 1
	2 3.5 1
	2 3.5 1
	2 3.5 1
	2 4.5 1
	2 4.5 1
	2 4.5 1
	2 5.5 1
	2 5.5 1
    2 5.5 1
	2 5.5 1
	2 5.5 1
	2 6.5 1
	2 7.5 1
	2 7.5 1
	2 7.5 1
	2 8.5 1
	2 8.5 1
	2 8.5 1
	2 9.5 1
	2 9.5 1
	2 10.5 1
	2 10.5 1
	2 10.5 1
	2 11.5 1
	2 11.5 1
	2 12.5 1
	2 12.5 1
	2 12.5 1
	2 12.5 1
	2 14.5 1
	2 14.5 1
	2 16.5 1
	2 16.5 1
	2 18.5 1
	2 19.5 1
	2 19.5 1
	2 19.5 1
	2 20.5 1
	2 22.5 1
	2 24.5 1
	2 25.5 1	
    2 26.5 1
	2 26.5 1
	2 28.5 1
	;
run;
ods select all;


proc print data=example noobs;
	title 'Example Data';
	title2 'Times to Infection (in months) of Kidney Dialysis Patients with Different Catheterization Procedures';
	title3 'From Klein and Moeschberger (1997)';
run;

ods rtf file='h:\temp\SAS\exampgraph.rtf';

proc lifetest data=example plots=(s) outsurv=out1 nocensplot censoredsymbol=' ';
    time months*censor(1);
   	strata treat;
   	label months='Survival Time in Days from Start of Treatment';
   	symbol1 H=1 W=1 L=1 Ci=blue V=square CV=blue;
	symbol2 H=1 W=1 L=1 Ci=red V=triangle CV=red;
	title 'Table 2';
	title2 'Survival Analysis for Times to Infection of Kidney Dialysis Patients';
	title3 'Kaplan-Meier Curves';
	ods output ProductLimitEstimates=est;
run;

ods rtf close;
*ods listing close;
*ods listing;

%include 'h:\statking\macros\flemharr.mac';

*%flemharr(odsname=est,time=months,censor=censor,numstrat=2,
		  stratless=1,dsd=example,p=0,q=0,tabno=1,Company=Test of F-H Test,
		  study=Test);


%macro none(odsname,time,censor,numstrat,stratless,dsd,p,q,tabno,company,study);

	%let currdate=%sysfunc(today(),worddate.);

	data est;
		set &odsname;
	run;

	/*proc print data=est; run;*/

	proc sort data=est; by stratum &time left; run;

	data counts; set est; by stratum &time left;
		if first.months;
	run;

	%do a=1 %to &numstrat;
		data counts&a; set counts; by stratum &time left;
			if stratum=&a;
			failed&a=failed;
			risk&a=left;
			if first.stratum then do;
				call symput("total&a",trim(left(risk&a)));
			end;
			drop censor left failed;
		run;

		proc sort data=counts&a; by stratum &time; run;

		data time&a; set counts&a; by stratum &time;
			if &time=0 then delete;
			count+1;
			keep stratum &time failed&a count;
		run;

		data risk&a; set counts&a; by stratum &time;
			count+1;
			keep stratum risk&a count;
		run;
	
		proc sort data=time&a; by count; run;
		proc sort data=risk&a; by count; run;

		data counts&a; merge time&a risk&a; by count;
			if months=. then delete;
			drop stratum count;
			flipvar=1;
		run;

		/*proc print data=counts&a; run;*/
		proc sort data=counts&a; by &time; run;
	%end;

	data counts; merge %do b=1 %to &numstrat; counts&b %end; %str(;) by &time;
	run; 

	proc sort data=counts; by flipvar; run;
	
	data counts; set counts; by flipvar;
		count+1;
		if first.flipvar then do;
			%do d=1 %to &numstrat;
				if risk&d=. then do;		
					risk&d=&&total&d;
					failed&d=0;
				end;
			%end;	
		end;
	run;

	data counts; set counts;
		call symput('count',trim(left(count)));
	run;

	proc sort data=counts; by flipvar; run;

	data risks; set counts; 
		keep &time risk1-risk&numstrat;
	run;

	proc sort data=est; by stratum &time; run;

	ods listing close;

	proc freq data=est;
		by stratum;
		tables &time*&censor;
		ods output CrossTabFreqs=failure;
	run;

	proc freq data=est;
		by stratum;
		tables &censor*&time*left;
		ods output CrossTabFreqs=size;
	run;

	ods listing;

	data failure; set failure;
		if &time=0 then frequency=0;
	run;
	/*proc print data=failure; run;*/
	proc sort data=failure; by &time; run;

	data failure1; set failure; by &time;
		if _type_=11;
		array counts fail1-fail&numstrat cens1-cens&numstrat;
		retain fail1-fail&numstrat cens1-cens&numstrat;
		if first.&time then do;
			do over counts; counts=.; end;
		end;

		%do g=1 %to &numstrat;		
			if (stratum=&g and censor=0) then fail&g=frequency;
			if (stratum=&g and censor=1) then cens&g=frequency;
		%end;

		keep &time fail1-fail&numstrat cens1-cens&numstrat; 
	run;
	proc sort data=failure1; by &time; run;

	data failure1; set failure1; by &time; 
		if last.&time;
	run;
	/*proc print data=failure1; run;*/
	
	proc sort data=risks; by &time; run;
	proc sort data=failure1; by &time; run;

	data total; merge failure1 risks; by &time;
		if &time=0 then delete;
	run;

	data total; set total; by &time;
		%do j=1 %to &numstrat;
			left&j=risk&j-(fail&j + cens&j);
		%end;
		count+1;
	run;

	data total1; set total;
		keep &time fail1-fail&numstrat cens1-cens&numstrat 
			 risk1-risk&numstrat count;
	run;

	data total2; set total;
		count=count+1;
		%do k=1 %to &numstrat;
			risk&k=left&k;
		%end;
		keep count risk1-risk&numstrat;
	run;

	proc sort data=total1; by count; run;
	proc sort data=total2; by count; run;

	data total; merge total1 total2; by count; 
		if &time=. then delete;
		flipvar=1;
		%do h=1 %to &numstrat;
			if fail&h=. then fail&h=0;
			if cens&h=. then cens&h=0;
		%end;
	run;

	proc sort data=total; by flipvar; run;
	/*proc print data=total; run;*/

	data risks; set total; by flipvar;
		array risks %do l=1 %to &numstrat; atrisk&l._1-atrisk&l._&count %end; %str(;)
		retain %do l=1 %to &numstrat; atrisk&l._1-atrisk&l._&count %end; %str(;)

		if first.flipvar then do;
			do over risks; risks=.; end;
		end;

		%do c=1 %to &count;
			if count=&c then do;
				%do m=1 %to &numstrat;
					atrisk&m._&c=risk&m;
				%end;
			end;
		%end;

		if last.flipvar then output;

		keep %do l=1 %to &numstrat; atrisk&l._1-atrisk&l._&count %end; %str(;)
	run;

	data risks; set risks; 
		%do e=2 %to &count;
			%let f=%eval(&e-1);
			%do n=1 %to &numstrat;
				if atrisk&n._&e=. then atrisk&n._&e=atrisk&n._&f;
			%end;
		%end;
	run;

	%do o=1 %to &numstrat;

	  data risks&o; set risks;
		array atrisks atrisk&o._1-atrisk&o._&count;
		retain atrisk&o._1-atrisk&o._&count;

		do over atrisks; count=_i_; atrisk&o=atrisks; output; end;
		keep count atrisk&o;
	  run;
	
	  proc sort data=risks&o; by count; run;

	%end;

	data risks; merge %do l=1 %to &numstrat; risks&l %end; %str(;); by count; run;

	proc sort data=risks; by count; run;
	proc sort data=total; by count; run;

	data total; merge total risks; by count;
		drop risk1-risk&numstrat count flipvar;
		totfail=sum(of fail1-fail&numstrat);
		totleft=sum(of atrisk1-atrisk&numstrat);
	run;

	/*proc print data=total; run;*/


	ods listing close;
	  
	proc lifetest data=&dsd outsurv=out2 censoredsymbol=' ';
		time &time*&censor(1);
		label &time='Survival Time in Days from Start of Treatment';
		title3 'Survival Analysis for Death from All Causes';
		ods output ProductLimitEstimates=est2;
	run;

	data est2;
		set est2;
		if &censor=0;
		if survival=. then delete;
		drop stratum left;
		count+1;
	run;

	data est2a; set est2;
		keep survival failure stderr failed count;
	run;

	data est2b; set est2;
	 	count=count-1;
		if count=0 the delete;
		keep &time censor count;
	run;

	proc sort data=est2a; by count; run;
	proc sort data=est2b; by count; run;

	data est2; merge est2a est2b; by count;
		if &time=. then delete;
	run;
	
	/*proc print data=est2; title 'Est 2'; run;*/
*ods listing;
	proc sort data=total; by &time; run;
	proc sort data=est2; by &time; run;
proc print data=total; title 'combine'; run;
proc print data=est2b; run;
proc print data=est2; title 'est2'; run;


	data newest;
		merge total est2; by &time;
	run;

	proc print data=newest;
		title 'New est';
		var &time fail1-fail&&numstrat atrisk1-atrisk&numstrat
			survival failure totfail totleft;
	run;
*ods listing close;
	data fleming;
		set newest;
		if &time>0;
		p=&p;
		q=&q;
		tiefix=((totleft-totfail)/(totleft-1))*totfail;
		%let counter2= ;
		%do n=1 %to &numstrat;
			expect&n=atrisk&n*(totfail/totleft);
			diff&n=fail&n-expect&n;
			percent&n=(atrisk&n/totleft);
			partvar&n=percent&n*(1-percent&n)*tiefix;
			if p=0 and q=0 then weight=1;
			else if p=0 and q>0 then weight=(failure)**q;
			else if p>0 and q=0 then weight=(survival)**p;
			else if p>0 and q>0 then weight=((survival)**p)*((failure)**q);
			wdiff&n=weight*diff&n;
			wvar&n=(weight**2)*partvar&n;
		%end;
		%do t=1 %to &numstrat;
			%do s=%eval(&t+1) %to &numstrat;
				partcov&t&s=(weight**2)*percent&t*percent&s*tiefix;
				%let counter2=%eval(&counter2+1);
				partcov&counter2=partcov&t&s;
			%end;
		%end;
	run;
ods listing;
	proc print data=fleming;
		title 'Weighting';
		var &time atrisk1-atrisk&numstrat totfail totleft expect1-expect&numstrat diff1-diff&numstrat 
			partvar1-partvar&numstrat weight wdiff1-wdiff&numstrat 
			wvar1-wvar&numstrat partcov1-partcov&counter2;
	run;
ods listing close;
	%let newcount2= ;
	%do r=1 %to &numstrat;
		proc means data=fleming noprint;
			var wdiff&r wvar&r;
			output out=try&r sum=z&r var&r;
		run;

		data try&r;
			set try&r;
			mervar=1;
		run;

		proc sort data=try&r; by mervar; run;

		%do u=%eval(&r+1) %to &numstrat;
			%let newcount2=%eval(&newcount2+1);
			proc means data=fleming noprint;
				var partcov&r&u;
				output out=harr&newcount2 sum=cov&r&u;
			run;
			data harr&newcount2;
				set harr&newcount2;
				cov&r&u=-cov&r&u;
				cov&u&r=cov&r&u;
			run;
		%end;
	%end;

	data tot;
		merge try1 try2; by mervar;
	run;
	
	proc sort data=tot; by mervar; run;

	%if &numstrat>2 %then %do;
		%do w=3 %to &numstrat;
			data tot;
				merge tot try&w; by mervar;
			run;
		%end;
	%end;
	
	proc print data=tot; title 'tot1'; run;

	%do v=1 %to &counter2;
		data tot;
			merge tot harr&v;
		run;
		/*proc print data=tot; title "tot2&v"; run;*/
	%end;

	data stats;
		set tot;
		keep z1-z&numstrat mervar;
	run;

	proc sort data=stats; by mervar; run;

	data stats;
		set stats; by mervar;
		array zs z1-z&numstrat;
		retain z1-z&numstrat;

		do over zs; treat=_i_; teststat=zs; output; end;
		drop z1-z&numstrat mervar;
	run;

	/*proc print data=stats; title 'stats'; run;*/

	data tot; set tot;
		%do mmm=1 %to &numstrat;
			cov&mmm&mmm=var&mmm;
		%end;
	run;

	%do nnn=1 %to &numstrat;
		data varcov&nnn;
			set tot;
			keep cov&nnn.1-cov&nnn&numstrat mervar;
		run;

		proc sort data=varcov&nnn; by mervar; run;

		data varcov&nnn;
			set varcov&nnn; by mervar;
			array covs cov&nnn.1-cov&nnn&numstrat;
			retain cov&nnn.1-cov&nnn&numstrat;

			do over covs; treat=_i_; cov&nnn=covs; output; end;
			drop cov&nnn.1-cov&nnn&numstrat mervar;
		run;
	%end;

	data varcov;
		merge varcov1 varcov2; by treat; 
	run;

	%if &numstrat>2 %then %do;
		%do rrr=3 %to &numstrat;
			data varcov;
				merge varcov varcov&rrr; by treat;
			run;
		%end;
	%end;

	data tot;
		set tot;
		drop _type_ _freq_ mervar;
	run;

/*	proc print data=tot;
		var treat z sigma typecov cov;
	run;*/

	data calc;
		set tot;
		%do x=1 %to &numstrat;
			call symput("var&x",var&x);
			call symput("z&x",z&x);
			%do y=%eval(&x+1) %to &numstrat;
				call symput("cov&x&y",cov&x&y);
				call symput("cov&y&x",cov&y&x); 
			%end;
		%end;
	run;

	proc print data=calc;
		title 'calc';
	run;

	%if &numstrat=2 %then %do;
		data one;
			set calc;
			zstat=z1/(var1**.5);
			chisq=(zstat**2);
			test='Fleming-Harrington';
			keep chisq test;
		run;
		/*proc print data=calc; run;*/
	%end;

	%else %if &numstrat>2 %then %do;
		proc iml;
			use calc;
			zs=j(1,&stratless,1);
			v=j(&stratless,&stratless,1);
			zs={&z1 %do aa=2 %to &stratless; &&z&aa %end;};
			v={&var1 %do bb=2 %to &stratless; &&cov1&bb %end; ,
				%if &stratless>2 %then %do;
					%do cc=2 %to %eval(&stratless-1);
						%let dd=%eval(&cc-1);
						%let ee=%eval(&cc+1);
						%do ff=1 %to ⅆ &&cov&cc&ff %end; 
						&&var&cc %do gg=&ee %to &stratless; &&cov&cc&gg %end;,
					%end;
				%end;
				%do hh=1 %to %eval(&stratless-1); &&cov&stratless&hh %end;
					&&var&stratless};
			invv=inv(v);
			chisq=zs*invv*zs`;
			create one var{chisq}; 
			append;
			title "Analysis for p=&p, q=&q";
			print zs, v, invv, chisq;
		quit;
		run;
	%end;

	/*
	proc print data=try1; run;
	*/
	
	data one;
		set one;
		test='Fleming-Harrington';
	run;

	data appear;
		merge stats varcov; by treat;
		treat2=put(treat,1.);
		%do ttt=1 %to &numstrat;
			covtrt&ttt=put(cov&ttt,8.5);
			drop cov&ttt;
		%end;
	run;

	/*proc print data=appear; title 'appear'; run;*/

	data appear;
		merge appear one;
		if chisq>0 then df=&stratless;
		else df=.;
		df2=put(df,6.);
		df2=trim(left(df2));
		pval=(1-probchi(chisq,&stratless));
		label teststat='Fleming-Harrington'
			  treat2='Stratum'
			  %do vv=1 %to &numstrat;
			  	covtrt&&vv="&&vv"
			  %end;
			  chisq='Chi-Square'
			  df='DF'
			  pval='Pr > Chi-Square'
			  test='Test';
	run;

	ods listing;
/*	footnote;

	proc print data=appear noobs label;
		title1 "&company";
		title2 "&study";
		title3 'Fleming-Harrington Test';
		title4 'For Differences in Survival Curves';
		title5 "Using p=&p, q=&q";
		title6 ' ';
		title7 ' '; 
		title8 ' ';
		title9 'Rank Statisics'; 
		var treat2 teststat;
	run;

	options formdlim=' ';
	
	proc print data=appear noobs label;
		title 'Covariance Matrix for the Fleming-Harrington Statistics';
		var treat2 covtrt1-covtrt&numstrat;
	run;

	proc print data=appear noobs label;
		title 'Test of Equality over Strata';
		var test chisq df pval;
	run;

	options formdlim='';
*/
/*	proc lifetest data=&dsd outs=stat;
    	time &time*&censor(1);
   		strata &treat;
   		label &time='Survival Time in Days from Start of Treatment';
		title1 "&company";
		title2 "&study";
   		title3 'Survival Analysis for Death from All Causes';
	run;
*/

	proc report data=appear nowindows headline split='*';
		columns (' Rank Statistics * ' treat2 teststat)
				(' Covariance Matrix * * * ' 
					treat2 covtrt1-covtrt&numstrat)
				(' Test of Equality over Strata * * ' chisq df2 pval);
		define treat2/ 'Stratum' /*page*/ width=7 center 
						/*format=&format*/;
		define teststat/ 'Fleming-*Harrington*Statistic' width=15 format=10.6 center;
		%do uuu=1 %to &numstrat;
			define covtrt&uuu/ "&uuu" width=11 center;
		%end;
		define chisq/ 'Chi-Square' width=11 center;
		define df2/ 'DF' width=7 center;
		define pval/ 'Pr >*Chi-Square' width=12 format=10.8 center;
		title1 "Table &tabno";
		title2 "&company";
		title3 "&study";
		title4 'Fleming-Harrington Test';
		title5 'For Differences in Survival Curves';
		title6 "Using p=&p, q=&q";
		title7 ' ';
		title8 ' '; 
		title9 ' ';
	run;


	%mend none;

%flemharr(odsname=est,time=months,censor=censor,numstrat=2,
		  dsd=example,p=0.5,q=2,tabno=1,Company=Example - Times to Infection of Kidney Dialysis Patients,
		  study=From Klein and Moeschberger (1997));