No AUC or cut-off value using %GLIMMROC

Tbobbo — Thu, 21 Feb 2019 17:21:54 GMT

%inc 'C:\myfolders\glimmix.sas'; 


%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_list;
      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;

%glimmroc(y=var1, x_list=var2 , /*z_list= */,ID=id, c_s_d=vc, c_s_r=csh, dataset=mydataset);run;

Hi all,

I'm trying to create a ROC curve for repeated-measures data, following the procedure published in: Liu, H., & Wu, T. (2003). Estimating the area under a receiver operating characteristic (ROC) curve for repeated measures design. J Stat Softw, 8(12), 1-18.

My outcome is binary (0/1) and my predictor is a continuous trait, measured 1 to 9 times for each individual. I've been playing around with %GLIMMROC and I am able to get the ROC curve; however, I could not get the AUC (I obtained a ROC1=1, but an AUC of 1 does not corresponds to the area) neither the optimal cutoff for the predictor.

To note, I obtained only 178 specificity values (12 missing values)..do you think that this could affect the AUC calculation?

Any suggestion? I will include the SAS script I have been using.

Thank you!

Re: No AUC or cut-off value using %GLIMMROC

StatDave — Mon, 25 Feb 2019 16:34:10 GMT

This can be done directly using the GLIMMIX and LOGISTIC procedures as shown in this note. Note that another modeling approach for repeated binary data is to use a GEE model available in PROC GEE or GENMOD. An ROC analysis can be obtained in the same way for the GEE model.

Re: No AUC or cut-off value using %GLIMMROC

Tbobbo — Wed, 27 Feb 2019 08:47:24 GMT

Thanks for your reply, I will take a look to the suggested procedures.

topic No AUC or cut-off value using %GLIMMROC in Statistical Procedures

No AUC or cut-off value using %GLIMMROC

Re: No AUC or cut-off value using %GLIMMROC

Re: No AUC or cut-off value using %GLIMMROC