options ps=55 ls=80;

data test15;

* data from Table 1 page 50 of ;
* ;
* Hanley JA and Hajian-Tilaki KO. ;
* Sampling Variability of Nonparametric Estimates of the Areas ;
* under Receiver Operating Characteristic Curves: An Update ;
* Academic Radiology, 1997 4:49-58. ;

input pt_no status r1 r2;
constant = 2;
lines;
1 1 1 1
2 0 2 1
3 1 5 5
4 0 1 1
5 0 1 1
6 1 1 1
7 1 2 4
8 0 1 1
9 0 2 2
10 1 2 2
11 0 1 1
12 0 1 1
13 1 5 5
14 0 1 1
15 0 1 1
;

* updated April 6, 2002

written by, and comments to:

J Hanley, McGill University email: James.Hanley@McGill.CA ;

*** specify dataset, truth and rating variables, polarity of each scale,
and the value which denotes criterion is present ;

%let dataset = test15; * name of dataset;

%let d = status; * name of truth variable;

%let positiv = 1 ; * value of truth variable denoting criterion present;

%let rating1 = r1 ; * variable containing ratings for <<test>> 1 ;

%let polar1 = 1 ; * 'polarity' of test1 scale ;
* positive value: higher test1 values INcrease Prob[criterion present] ;
* negative value: higher test1 values DEcrease Prob[criterion present] ;

%let rating2 = r2 ; * variable containing ratings for <<test>> 2 ;

%let polar2 = 1 ; * 'polarity' of test2 scale ;
* positive value: higher test1 values INcrease Prob[criterion present] ;
* negative value: higher test1 values DEcrease Prob[criterion present] ;

* If wish to obtain SE for a single auc, or single c ;
* use the same variable for rating1 and rating2 ;

* If let rating1 = a variable that is actually constant, ;
* then AUC1 = 0.5 (SE1=0) and CI(diff) is for (AUC2 - 0.5) ;
* To get CI(AUC2), add 0.5 to l95_diff and u95_diff ;

* *******************************************************************;

data temp_a; keep d &rating1 &rating2;
set &dataset; d = 0;
if &d = &positiv then d = 1;
if &polar1 < 0 then &rating1 = 0 - &rating1;
if &polar2 < 0 then &rating2 = 0 - &rating2;

proc means noprint data=temp_a;
var d ; output out=temp_n n = n_tot sum = n_d1;
*proc print data=temp_n; run;

proc rank data = temp_a out = temp_b;
var &rating1 &rating2; ranks overall1 overall2;
*proc print; run;

proc sort data=temp_b; by d;
proc rank data=temp_b out = temp_c;
var &rating1 &rating2; ranks spcific1 spcific2; by d;
run;

data temp_d;
retain count 0 n_d0 n_d1;
if count = 0 then do;
set temp_n;
n_d0 = n_tot - n_d1;
count = 1;
end;
set temp_c;
if d=1 then do;
locatn1 = (overall1-spcific1)/n_d0;
locatn2 = (overall2-spcific2)/n_d0;
end;
if d = 0 then do;
locatn1 = 1-(overall1-spcific1)/n_d1;
locatn2 = 1-(overall2-spcific2)/n_d1;
end;
product = locatn1*locatn2;
*proc print data=temp_d;
run;

*proc sort data=temp_d;
*by d;

proc means noprint data=temp_d mean n var;
var locatn1 locatn2 product; by d;
output out = temp_var
mean = m1 m2 m12
var = v1 v2 v12
n = n1 n2 n12;

*proc print data=temp_var;

data temp_cov;
set temp_var;
cov12 = ( n1/(n1-1) ) * (m12 - m1*m2) ;
v1 = v1/n1; v2 = v2/n2; cov12 = cov12/n1;
run;
*proc print data=temp_cov;
proc means noprint data=temp_cov;
var m1 m2 v1 v2 cov12;
output out = temp_sum
mean=auc1 auc2 m mm mmm
sum = s sss v1 v2 cov12 ;
*proc print data=temp_sum;
run;

data temp_inf;
keep auc1 auc2 se1 se2 r_auc1_2 a2minus1 se_diff
l95_diff u95_diff z_diff p_diff;
set temp_sum;
auc1 = round(auc1,0.001);
auc2 = round(auc2,0.001);
se1 = round(sqrt(v1),0.001); se2 = round(sqrt(v2),0.001);
a2minus1 = round(auc2-auc1,0.001);
se_diff = round(sqrt(v1 + v2 - 2*cov12),0.001);
l95_diff = round(a2minus1 - 1.96*se_diff,0.001);
u95_diff = round(a2minus1 + 1.96*se_diff,0.001);
z_diff = round(a2minus1 / se_diff, 0.001);
p_diff = round(2 * ( 1-probnorm(abs(z_diff)) ),0.001) ;
r_auc1_2 = round(cov12 / (se1 * se2), 0.01);

put / ' NOTATION..... '
// ' auc1, auc2 : (trapezoidal) areas under curves '
// ' SE1 , SE2 : their standard errors '
// ' a2minus1 : auc2 minus auc1 '
// ' SE_diff : SE of this difference '
// ' l95_diff : lower limit of 95% CI for true AUC difference'
/ ' u95_diff : upper limit of 95% CI for true AUC difference'
// ' z_diff : z statistic = (auc2 - auc1)/SE_diff '
// ' p_diff : corresponding 2-sided p-value [H0: AUC2 = AUC1 ] '
// ' r_auc1_2 : correlation(auc1,auc2) '
//
;
;
proc print data =temp_inf noobs heading=h ;
run;