/*******************************************************************************************/
/*  MACRO to calculate family history scores using the methods by Feng et al. (2008)       */ 
/*  Feng R., McClure L., Tiwari H.K., Howard G. (2008). Statistics in medicine. In review. */
/*                                                                                         */
/*                    EXAMPLE OF USAGE                                                     */
/*   %include 'c:\local directory\SLFS.sas';                                               */
/*   data testset;                                                                         */
/*   infile 'C:\local directory\samplefamily.txt' delimiter='09'x  firstobs=2;             */   
/*   input familyID	membername $3.	event	futime;                                        */
/*   RUN;                                                                                  */
/*   %SLFS(testset, familyID, group, event,	futime, out=outputdata);                       */                                               
/*   with following arguments:                                                             */
/*  testset is the dataset you wish to analyze                                             */
/*  familyId is a variable to index families within the dataste                            */
/*  group is a variable to define a group (or pointer) within the dataset                  */
/*  event is a variable of the status of event: 0-censored and 1-observed event            */
/*  futime is a variable of age of onset or age at censor                                  */
/*  out= is an optional augument, results will be saved in "resultSLFS" if not specified   */
/*  You have to change the groupname1-groupname4 if the values in the group variable       */
/*  are not "Dad" "Mom" "Bro" and "Sis"; and you can change number of groups be adding or  */
/*  deleting groupnamexx, where xx is a number.                                            */
/*  This is trial version - use at your own risk.                                          */
/*  For questions and bugs, email Rui Feng at rfeng@ms.soph.uab.edu                        */  
/*******************************************************************************************/
%MACRO SLFS(dataset, familyID, group, age, status, out=resultSLFS);
options NOMLOGIC NOMPRINT NOMRECALL NOSYMBOLGEN;
%let groupname1=Dad;  /*can be changed to other group value*/
%let groupname2=Mom;  /*can be changed to other group value*/
%let groupname3=Bro;  /*can be changed to other group value*/
%let groupname4=Sis;  /*can be changed to other group value*/
data &dataset; set &dataset;
 id=_n_;
Run;
proc sort data=&dataset;
by &familyID;
run;
data uniquefamid;    /*store unique family IDs in order to match SFLS with each family*/
   set &dataset;
   by &familyID;
   if first.&familyID=1 then output;
   keep &familyID;
run; 

%do i=1 %to 4; /*should match with number of groups in the dataset*/
   data data&i; set &dataset;  /*get all individuals within the current group*/
       where &group="&groupname&i";
   run;
   data timetmp; set data&i;  /*next 12 lines - find distinct time points of observed events*/
    where &status=1;           
    keep &age;                
   run;
   proc sort data=timetmp;
   by &age;
   run;
   data timepoint&i;
   set timetmp;
   by &age;
   if first.&age=1 then output;
   run;
%end;

PROC IML;
use &dataset;
	read all var {&familyID} into vfamilyID; 
close &dataset;
n=nrow(vfamilyID);               *Total number of subjects;

use uniquefamid;
      read all var {&familyID} into vuniquefID;
close uniquefamid;;
totalfam=nrow(vuniquefID);      *Total number of families;

risk=J(n,1,.); *log-rank score for each individual;

%do i=1 %to 4; *for each group, the risk score is calculated seperately;
  use data&i;
	read all var {&familyID, &age, &status, id} into datamatrix; /*storing data in matrix for calculation*/
  close data&i; 
  nsub=nrow(datamatrix);
  use timepoint&i;
      read all var {&age} into period;    *define time period to be unique event time point ;
  close timepoint&i;
  nperiod=nrow(period);
  n0=J(nperiod+1,1, 0); * n0j is the number of right censored observations in the ith time interval;
  n1=J(nperiod+1,1, 0); * n1j is the number of observed events in the ith time interval;
  do j=1 to nsub;
     g=1;
	 if (datamatrix[j,2]>period[1]) then do; 
          k=2;
          do while ((datamatrix[j,2]>period[k]) && (k<nperiod) );
			 k=k+1;
	      end;
		  g=k+1;
     end;  
	 if (datamatrix[j,3]=0) then n0[g]=n0[g]+1;
     if (datamatrix[j,3]=1) then n1[g]=n1[g]+1; 
  end;

  a0=J(nperiod+1, 1,  0); *risk score for censored observations with age at censor within the ith time interval;
  a1=J(nperiod+1, 1, 0); *risk score for observed events with age-on-set within the ith time interval;
  do j=1 to nperiod;
     do k=1 to j;
		tmp1=n0[k:nperiod];
		tmp2=n1[k:nperiod];
        tmp=tmp1[+]+tmp2[+]; 
        if(tmp>0) then a0[j]=a0[j]-n1[k]/tmp;
     end;
     a1[j]=a0[j]+1;
  end;
  a1[nperiod+1]=a1[nperiod];
  a0[nperiod+1]=a1[nperiod];
  do j=1 to nsub;    /*for all individuals*/
     g=1;           /*next 8 lines- find the current individual's corresponding ageonset or age-at-censor group*/
	 if (datamatrix[j,2]>period[1]) then do; 
          k=2;
          do while ((datamatrix[j,2]>period[k]) && (k<nperiod) );
			 k=k+1;
	      end;
		  g=k+1;
     end; 
     currentid=datamatrix[j,4]; 
	 if (datamatrix[j,3]=0) then risk[currentid]=a0[g]; /*assign corresponding risk to a censored individual*/
     else if (datamatrix[j,3]=1) then risk[currentid]=a1[g]; /*assign risk to the individual with an observed event*/
  end;  
%end; 
est_risk1=J(totalfam, 1,  .); *a vector to store SLFS for all families;
k=1;
do j=1 to totalfam;
     tmpid=datamatrix[k,4];   /*get current individual id*/ 
     tmp=risk[tmpid];         /*current individuals' risk score*/
	 k=k+1;
     do while (datamatrix[k,4]=tmpid); /*if an individual is within the current family*/
        tmpid=datamatrix[k,4];  
        tmp=tmp//risk[tmpid]; /*store this individual's risk score in a vector "tmp" */
        k=k+1;
	 end;
     est_risk1[j]=tmp[:]; /*average of the risk scores of all individuals within j th family*/
end;
resultm=vuniquefID||est_risk1;
cname={"familyID" "SLFS"};
create &out from resultm[colname=cname];
   append from resultm;
close &out;
quit;
%MEND SLFS;
/* Try on a dataset imported from a text example file; 
data testset;
infile 'C:\Documents and Settings\rfeng\My Documents\2007\familyhistory\samplefamily.txt' delimiter='09'x  firstobs=2; 
input familyID	membername $	event	futime;
membername=substr(membername,1,3);
RUN;
%SLFS(testset, familyID, membername, futime, event);
*/