%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!
