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)
Tang, L. and N. Balakrishnan. “A Random-Sum Wilcoxon Statistic with Applications to Various Types of Diagnostic Data." Journal of Statistical Planning and Inference. 141 (1), 335-344. R code for estimating LROC curve and ROC curve with the limit of detection:
#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)