/*********************************************************************************
%GLIMMROC (y=,x_list=,z_list=,id=,c_s_d=, c_s_r=,weight=,dataset=);

Purpose: To generate receiver operating characteristic (ROC) curve and
         to calculate the area under the curve through generalized
         linear mixed modeling with a Wilcoxon non-parametric approach
         for repeated measures longitudinal design

y----------the variable name of the binary outcome
x_list-----contains the list of all independent variable for the fixed
           effects with space in between(e.g., age sex race)
z_list----contains the list of all independent variables for the random
          effects with space in between (e.g., time) 
ID---------the variable of patient ID which identifies observations
           within a patient
c_s_d------specify the covariance structure of matrix D (for the random
           part of GLMM and the default is simple, e.g., diagonal
           matrix.)
c_s_r------specify the covariance structure of matrix R (for the 
           error of the model and the default is compound symmetry)
weight-----weight variable (the default is 1, which means unweighted)
dataset----the name of the data set on which one will run the
           macro procedure.

Output:    ROC curve and the statistic of the area under the curve    
****************************************************************************/

%inc 'c:\tt\ROC\glmm800.sas'; /***need to invoke %glimmix macro from SAS. The actual subdirectory name will be up to specific user*/

%macro glimmroc(y=,x_list=,z_list=,id=,c_s_d=, c_s_r=,weight=,dataset=);
%******* get working dataset DATAH ready *************;
      %if %length(&dataset)>0 %then %do;
data datah; set &dataset; rename &id=id;/*rename &y=obs;*/
                                   %end;
%else %if %length(&dataset)=0 %then %do;
data datah; set _last_; rename &id=id; /*rename &y=obs;*/
                                   %end;

%******* get z_list *************;
      %if %length(&z_list)>0 %then %do;
%let z=intercept%str( )&z_list;
                                   %end;
%else %if %length(&z_list)=0 %then %do;
%let z=intercept;                           
                                   %end;

%******* get c_s_d *************;
      %if %length(&c_s_d)>0 %then %do;
%let csd=&c_s_d;
                                   %end;
%else %if %length(&c_s_d)=0 %then %do;
%let csd=simple;                           
                                   %end;
%******* get c_s_r *************;
      %if %length(&c_s_r)>0 %then %do;
%let csr=&c_s_r;
                                   %end;
%else %if %length(&c_s_r)=0 %then %do;
%let csr=cs;                           
                                   %end;

%let wt   = %qupcase(&weight);
%if %length(&weight)=0 %then %do;
%glimmix(data=datah,
   stmts=%str(
      class &id;
      model &y = &x_list;
      repeated/ type=&csr subject=&id;
	  random &z / type=&csd; 
	           ),
 error=binomial, options=noprint

) ;
%end;
%else %do;

%glimmix(data=datah,
   stmts=%str(
      class &id;
      model &y = &x;
      repeated/ type=&csr subject=&id;
	  random &z / type=&csd; 
	           ),
 error=binomial, weight=&wt, options=noprint

) ;
%end; 


%*******get the predicted probability****************;
data pred;
   set _ds;
   /*ptid=&id;*/
   obs=&y;
run;

data data0;
  set pred;
  if obs=0;
run;
data data1;
  set pred;
  if obs=1;
  Rsum1=0;
run;

data temp;
  set pred;
  ROC1=0;
run;
data out;
  set temp;
  keep ROC1;
run;



data _dataa_; set data0 (keep=id mu obs); run;
data _datab_; set data1 (keep=id mu obs Rsum1); run;
data _out_; set out; run;

/* start iml: interative matrix language evironment */
proc iml;                                          
     use _dataa_; read all into dataa;
     use _datab_; read all into datab;
     use _out_; read all into out;

	    start ROC(dataa,datab,out);
			nra = nrow(dataa);
			nrb = nrow(datab);
			Rsum1=0;

			if nrb < nra then do;
				do i=1 to nrb;
				temp = datab[i,2];
				n1 = ncol(loc(dataa[,2] < temp));
				n2 = ncol(loc(dataa[,2] = temp));
				Rsum1 = Rsum1 + n1 + .5*n2;
				end;
			end; 

            else do;
				do i=1 to nra;
				temp = dataa[i,2];
				n1 = ncol(loc(datab[,2] > temp));
				n2 = ncol(loc(datab[,2] = temp));
				Rsum1 = Rsum1 + n1 + .5*n2;
				end;
			end;

			roc1 = Rsum1/(nra*nrb);
			print roc1;
			out[nrow(out),1] = roc1;
			return(out);
		finish ROC;



R=ROC(dataa,datab,out);
varnames={'ROC1'};
create outdata from R [colname=varnames] ;
append from R;
quit;

%*********get ROC curve*************;
data yy;
   set pred;
   /*if pred^=.;*/
   do i=1 to 200;
   if mu>0.005*i then y=1;
   else if .z<mu<0.005*i then y=0;
   else y=.;
   output;
   end;
run;
proc sort data=yy; by i; run;


%********get sensitivity and specificity*************;
proc freq data=yy;
tables y*obs/noprint out=pct outpct;
by i;
run; 

data sen;
   set pct;
   if y=0 and obs=0;
    sensi=PCT_COL;
   cut=0.005*i;
run;


data spc;
   set pct;
   if y=1 and obs=1;
    speci=PCT_COL;
    cut=0.005*i;
run;


data curve1;
    merge sen spc;
    by cut;
    _spc=100-speci;
run; 
data curve;
  set curve1;
  sensi1=sensi/100;
  _spc1=_spc/100;
run;

proc sort data=curve;
by _spc1;
run;

title "ROC curve";
symbol1 i=join v=star line=3 c=red;
axis1   order=(0 to 1 by 0.1);
axis2   order=(0 to 1 by 0.1) label=(a=90);

proc gplot data=curve;
plot sensi1*_spc1=1/vaxis=axis2 haxis=axis1;
label sensi1='Sensitivity';
label _spc1='1-Specificity';
run;
quit;

%mend glimmroc;