Software

The projects on ROC and likelihood ratio are currently supported by Department of Justice (NIJ grant 2018-DU-BX-0228, Tang as contact PI, Danica Ommen at Iowa State University as Joint PI, Chris Saunders as subcontract PI at South Dakota University; NIJ grant 2019-DU-BX-0011, Tang as PI).

  • R codes and data used in “Order-Restricted Inference for Correlated and Clustered ROC Data with Application to Fingerprint Data” by Zhang, Tang, Li, Tabassi, and Lee. (supported by NIJ grant 2018-DU-BX-0228)

#Both sample sizes are random


library(MASS)

# modified from Delong's SAS code

mwcomp<-function(z)

{

rz<-rank(z[,1]) #average ranks

nx<-sum(z[,2]) #num. of x's

ny<-nrow(z)-nx #num. of y's

loc<-which(z[,2] == 1, arr.ind=TRUE) #x indexes

psi<-matrix(0,nrow(z),1) # creates a matrix of identical values

psi[loc]<-(rz[loc]-rank(z[loc,1]))/ny # x components

loc<-which(z[,2] == 0, arr.ind=TRUE) # y indexes

psi[loc]<-(nx+rank(z[loc,1])-rz[loc])/nx # y components

return(psi)

}

mwvarp23<-function(z)

{

ind<-z[,2] #indicator values

who<-which(is.na(z[,1] *z[,1] ) != TRUE, arr.ind=TRUE) #nonmissing pairs

psii<- mwcomp(cbind(z ,ind)[who,]) #components V^{i}

inow<-ind[who] #x's and y's

tinow<-t(inow)

inow<-t(tinow)

m<-sum(inow) #current x's

n<-nrow(psii)-m #current y's

mi<-sum(psii*inow)/m #means

psii<-psii #center

p2<-sum(psii*psii*inow)/ (m-1)

p3<-sum(psii*psii*(1-inow))/ (n-1)

return(c(p1=mi,p2, p3))

}

## generate ordinal data

ordscale <- function(q0, val) {

valn <- c()

for (i in 1:length(val))

for(j in 1:(length(q0)-1)){

if (val[i] >=q0[j]&&val[i] < q0[j+1]) valn[i] <- j

}

return(valn)

}

# read in the data in this step, the last column indicates the disease status,

# 1 for diseased, 0 for non-diseased

# I used the simulated data here

# sample sizes 20 for diseased, 30 for nondiseased, 7 biomarkers,

cancerdata<-matrix(0, nrow=20+30, ncol=7+1)

cancerdata[,8]<-c(rep(1,20), rep(0,30))

for(i in 1:7) cancerdata[,i]<-c(rnorm(20 ,i ,1),rnorm(30))

m<-nrow(cancerdata[cancerdata[,8]==1,])

n<-nrow(cancerdata[cancerdata[,8]==0,])

N=m+n

auclod<-rep(0,7)

aucvar<-rep(0,7)

for(i in 1:7){

datai<-cancerdata[,c(i,8)]

# Suppose LOD= 0.1

data.temp<-datai[datai[,1] >0.01,]

X<-data.temp[data.temp[,2]==1,]

Y<-data.temp[data.temp[,2]==0,]

Y[,2]<-0

k1<-length(X[,1])

k2<-length(Y[,1])

# create z data matrix

z<-rbind( X , Y )

p1.hat<-k1/m

p2.hat<-n/n

p123<- mwvarp23(z)

p1<-p123[1]

emp.R <- p1*length(X[,1])/m*length(Y[,1])/n

p2<-p123[2]

p3<-p123[3]

auclod[i]<-(1-k1/m)*(1/2+k2/(2*n))+emp.R

aucvar[i]<- 1/4*(p1.hat*(1-p1.hat)*(1+p2.hat)^2/m+p2.hat*(1-p2.hat)*(1-p1.hat)^2)+ p1.hat*p2.hat^2*(p2-p1^2)/m+p1.hat^2*p2.hat*(p3-p1^2)/n+( m*n^2*p1.hat*(1-p1.hat)*p2.hat^2+m^2*n*p1.hat^2*p2.hat*(1-p2.hat)) *p1^2/(m*n)^2

}

# AUCs and their variances

result<-cbind(auclod, aucvar)