## Here all all the functions to define before you run the code for the analyses## covariance_matrix<-function(binarydata, nR, nN, nD){ correlation_matrixY<-matrix(nrow = nN, ncol=2*nR); correlation_matrixX<-matrix(nrow = nD, ncol=2*nR); for (j in 1:nR){ for (i in 1:nN){ correlation_matrixY[i,j]=mean(binarydata[binarydata[,3]==i & binarydata[,5]==j,1]); correlation_matrixY[i,(j + nR)]=mean(binarydata[binarydata[,3]==i & binarydata[,5]==j,2]); } for (i in 1:nD){ correlation_matrixX[i,j]= mean(binarydata[binarydata[,4]==(nN+i) & binarydata[,5]==j,1]); correlation_matrixX[i,(j + nR)]= mean(binarydata[(binarydata[,4]==(nN+i)) & (binarydata[,5]==j),2]); } } Theta_hat <-matrix(nrow = 1, ncol = 2*nR); for (j in 1:nR){ Theta_hat[j]=mean(binarydata[binarydata[,5]==j,1]); Theta_hat[j + nR]=mean(binarydata[binarydata[,5]==j,2]); } S10 = (1/(nD-1))*(t(correlation_matrixX)%*%correlation_matrixX - nD*t(Theta_hat)%*%Theta_hat); S01 = (1/(nN-1))*(t(correlation_matrixY)%*%correlation_matrixY - nN*t(Theta_hat)%*%Theta_hat); correlation_matrix<-(1/nD)*S10 + (1/nN)*S01; return(correlation_matrix); } gen.data<-function(Y, nN, nD, nR){ output<- matrix(nrow = nR*nN*nD, ncol = 5); num_pts = nN+nD; #Make the data set #fill in reader, normal case, and diseased case numbers for(i in 1:nR){ #go through the readers count = nN*nD*(i - 1);#for indexing subset<-Y[Y[,5]==i, ]; #data just for reader i #reader number in 5th column: output[(count+1):(count+ nN*nD), 5] = i; for (j in 1:nN){ count.2<-(j-1)*nD; #case number k in 3rd column, s in 4th column: output[(count+count.2+1):(count + count.2+ nD), 3] = j; output[(count+count.2+1):(count + count.2+nD), 4] = (nN+1):(nD+nN); #fill in indicators output[(count + count.2 + 1:nD), 1] = (subset[j, 1] < subset[(nN + 1:nD), 1]) + .5*(subset[j, 1] ==subset[(nN + 1:nD), 1]); output[(count + count.2 + 1:nD), 2] = (subset[j, 2] < subset[(nN + 1:nD), 2]) + .5*(subset[j, 2] ==subset[(nN + 1:nD), 2]); } } return(output) ; } r2sigmac_hat<-function(correlation_matrix, nR){ #different j correlation matrix1<-correlation_matrix[1:nR, 1:nR]; matrix2<-correlation_matrix[(nR+1):(2*nR), (nR+1):(2*nR)]; sum1 = 0; sum2 = 0; num = (nR-1)*nR/2; for (i in 1:(nR-1)){ sum1= sum1 + sum(matrix1[i, (i+1):nR]); sum2= sum2 + sum(matrix2[i, (i+1):nR]); } #corr1=(sum1 + sum2)/num return(c(sum1/num, sum2/num)); } r3sigmac_hat<-function(correlation_matrix, nR){ #different i, j correlation matrix1<-correlation_matrix[1:nR, (nR+1):(2*nR)]; sum = 0; num = (nR-1)*nR; for (i in 1:(nR-1)){ sum= sum + sum(matrix1[i, (i+1):nR]); } matrix2<-t(matrix1); for (i in 1:(nR-1)){ sum= sum + sum(matrix2[i, (i+1):nR]); } corr=sum/num; return(corr); } OR_MST<- function(PS, nR, num_pts){#PS is pseudovalues matrix MST = 0; AveALL = mean(c(PS[,1], PS[,2])); for (i in 1:2){ Ave = mean(PS[, i]); MST = MST + (Ave -AveALL)^2; } MST=MST*nR; return(MST) } OR_MSRboth<- function(PS, nR){#PS is pseudovalues matrix MST = 0; AveALL = mean(c(PS[,1], PS[,2])); for (i in 1:nR){ Ave = mean(c(PS[PS[,4]==i, 1], PS[PS[,4]==i, 2])); MST = MST + (Ave -AveALL)^2; } MST=MST*2/(nR-1); return(MST) } OR_MSRsingle<- function(PS, nR){#PS is pseudovalues matrix MSR = c(0,0); Ave1 = mean(PS[,1]); Ave2 = mean(PS[,2]); for (i in 1:nR){ Ave11 = mean(PS[PS[,4]==i, 1]); Ave22 = mean(PS[PS[,4]==i, 2]); MSR[1] = MSR[1] + (Ave11 -Ave1)^2; MSR[2] = MSR[2] + (Ave22 -Ave2)^2; } return(MSR/(nR-1)) } OR_MSTR<- function(PS, nR, num_pts){ #PS is pseudovalues matrix SUM = 0; AVE = mean(c(PS[,1], PS[,2])); #Y dot dot dot SumIJ = c(0,0); for (i in 1:2){ Ave.i = mean(PS[,i]); #Y i dot dot for(j in 1:nR){ Ave.ij = mean(PS[PS[,4]==j, i]); # Y i j dot Ave.j = mean(c(PS[PS[,4]==j, 1], PS[PS[,4]==j, 2])); #Y dot j dot SumIJ[i] = SumIJ[i] + (Ave.ij - Ave.i - Ave.j + AVE)^2; } } SUM = SumIJ[1] + SumIJ[2]; MSTR=SUM/(nR-1); return(MSTR); } DDFH<-function(mstr, mstc, mstrc, nR){ max<- max(mstc-mstrc, 0); h<-(max+mstr)^2; h<-h/((mstr^2)/(nR-1)); return(h); } DDFH_OR<-function(mstr, nR, cov2, cov3){ maximum<-max(nR*(cov2-cov3), 0); h<-(maximum+mstr)^2; h<-h/((mstr^2)/(nR-1)); return(h); } DDFH_single<-function(msr, cov2, nR){ output<-c(0,0) temp<-c(0,0) for(i in 1:2){ temp[i]<-(msr[i]+max(nR*cov2[i], 0))^2; otheri<-msr[i]^2 output[i]=(nR-1)*temp[i]/otheri } return(output) } #single_trt_CI_radius<-function(msr, r2sigc, nR, dof){ #sum1 = 0; sum2 = 0; #theta_hat1= mean(binarydata[,1]); theta_hat2= mean(binarydata[,2]); #for (j in 1:nR){ # theta_hat1j= mean(binarydata[binarydata[,5]==j,1]); # theta_hat2j= mean(binarydata[binarydata[,5]==j,2]); # sum1 = sum1 + (theta_hat1j - theta_hat1)^2; # sum2 = sum2 + (theta_hat2j - theta_hat2)^2; #} #sum1 = sum1/(nR*(nR-1)); sum2 = sum2/(nR*(nR-1)); #t = qt(1-(0.05/2), dof); #bottom<-c(0,0); #bottom[1]<-msr[1]+max(r2sigc[1]*nR,0); #bottom[2]<-msr[2]+max(r2sigc[2]*nR,0); #ciradius<-c(t[1]*sqrt(sum1 + max(r2sigc[1],0)),t[2]*sqrt(sum2 + max(r2sigc[2],0))); #ciradius<-c(t[1]*sqrt((1/nR)*bottom[1]),t[2]*sqrt((1/nR)*bottom[2])); #return(ciradius); #} single_trt_CI_radius<-function(msr, r2sigc, nR, dof){ t = qt(1-(0.05/2), dof); bottom<-c(0,0); bottom[1]<-msr[1]+max(r2sigc[1]*nR,0); bottom[2]<-msr[2]+max(r2sigc[2]*nR,0); ciradius<-c(t[1]*sqrt((1/nR)*bottom[1]),t[2]*sqrt((1/nR)*bottom[2])); return(ciradius); } trt_difference_CI_radius<-function(mstr, r2sig, r3sig, nR, dof){ #trt1 - trt2 #sum = 0; #theta_hat1= mean(binarydata[,1]); theta_hat2= mean(binarydata[,2]); #for (j in 1:nR){ # theta_hat1j= mean(binarydata[binarydata[,5]==j,1]); # theta_hat2j= mean(binarydata[binarydata[,5]==j,2]); # sum = sum + ((theta_hat1j - theta_hat2j)-(theta_hat1 - theta_hat2))^2; #} #sum = sum/(nR*(nR-1)); #sum = sum + 2*max((r2sigc-r3sigc), 0); #sum = sqrt(sum); t = qt(1-(0.05/2), dof); bottom<-(2/nR)*(mstr + max(nR*(r2sig-r3sig),0)) return(t*sqrt(bottom)) } ORFstat<-function(r2sigma, r3sigma, mst, mstr, nR, dof){ f = mst/(mstr + nR*max((r2sigma - r3sigma),0)); return((f > qf(.95,1,dof))); } gen.data.COUNT.TIES<-function(Y, nN, nD, nR){ output<- matrix(nrow = nR*nN*nD, ncol = 5); num_pts = nN+nD; #Make the data set #fill in reader, normal case, and diseased case numbers for(i in 1:nR){ #go through the readers count = nN*nD*(i - 1);#for indexing subset<-Y[Y[,5]==i, ]; #data just for reader i #reader number in 5th column: output[(count+1):(count+ nN*nD), 5] = i; for (j in 1:nN){ count.2<-(j-1)*nD; #case number k in 3rd column, s in 4th column: output[(count+count.2+1):(count + count.2+ nD), 3] = j; output[(count+count.2+1):(count + count.2+nD), 4] = (nN+1):(nD+nN); #fill in indicators for (m in 1:nD){ output[count + count.2 + m, 1] = as.numeric(subset[subset[,4]==output[count + count.2 + m, 3], 1] < subset[subset[,4]==output[count + count.2 + m, 4], 1]) + .5*as.numeric(subset[subset[,4]==output[count + count.2 + m, 3], 1] ==subset[subset[,4]==output[count + count.2 + m, 4], 1]); output[count + count.2 + m, 2] = as.numeric(subset[subset[,4]==output[count + count.2 + m, 3], 2] < subset[subset[,4]==output[count + count.2 + m, 4], 2]) + .5*as.numeric(subset[subset[,4]==output[count + count.2 + m, 3], 2] ==subset[subset[,4]==output[count + count.2 + m, 4], 2]); } } } return(output) ; } AUC_est <- function(Y, nN, nD){ est = matrix(nrow = 1, ncol = 2); for (j in 1:2){ est[1, j] = wilcox.test(Y[Y[,3]==1, j], Y[Y[,3]==0, j])$statistic; } est = est/(nN*nD); } pseudos3<- function(Y, nN, nD, nR){ num_pts = nN+nD; ps = matrix(nrow = nR*(num_pts), ncol = 4); ps[,3]=Y[,4]; # case ps[,4]=Y[,5]; #reader for (i in 1:nR){ k = num_pts*(i-1); temp1<-Y[Y[,5]==i, ]; # data for only reader i auc_i <- AUC_est(temp1, nN, nD); for (j in 1:num_pts){ temp2 <- temp1[temp1[,4]!=j,]; d = temp1[j,3]; #indicator of disease in jth patient n = 1-d; #indicator of normal ps[j+k, 1:2] = num_pts*auc_i - (num_pts-1)*AUC_est(temp2, nN-n, nD-d); } } return(ps); } getINF <-function(data, xbeta){ #dU<- matrix(nrow = nR*nN*nD, ncol = 5); dU<- matrix(nrow = dim(data)[1], ncol = 5); temp1<- matrix(nrow = dim(data)[1], ncol = 2); temp2<- matrix(nrow = dim(data)[1], ncol = 2); dU[,3:5]=data[,3:5]; # column 3-> reader, 4-> normal case, 5-> diseased case dnorm.1 = dnorm(xbeta[1]); dnorm.2 = dnorm(xbeta[2]); phi.1 = pnorm(xbeta[1]); phi.2 = pnorm(xbeta[2]); INF<-matrix(nrow = 2, ncol = 2); temp1[,1]<- (data[,1] - phi.1)/(phi.1*(1-phi.1)); temp1[,2]<- (data[,2] - phi.2)/(phi.2*(1-phi.2)); temp2[,1]<- (-1)*(temp1[,1]^2); temp2[,2]<- (-1)*(temp1[,2]^2); dU[,1]<- temp2[,1]*(dnorm.1^2) - temp1[,1]*xbeta[1]*dnorm.1; dU[,2]<- temp2[,2]*(dnorm.2^2) - temp1[,2]*xbeta[2]*dnorm.2; INF[1,1] = sum(dU[,1], na.rm = TRUE) + sum(dU[,2], na.rm = TRUE); INF[1,2] = sum(dU[,2], na.rm = TRUE); INF[2,1] = INF[1,2]; INF[2,2] = sum(dU[,2], na.rm = TRUE); return(INF); } getB <- function(data, xbeta, nN, nD, nR){ B <- matrix(c(0,0,0,0), 2, 2); u<- matrix(nrow = dim(data)[1], ncol = 5); u[,3:5]=data[,3:5]; # column 3-> normal case, 4-> diseased case, 5->reader phi.1 = pnorm(xbeta[1]); phi.2 = pnorm(xbeta[2]); dnorm.1 = dnorm(xbeta[1]); dnorm.2 = dnorm(xbeta[2]); u[,1] <-((data[,1]-phi.1)*dnorm.1)/(phi.1*(1-phi.1)); u[,2] <-((data[,2]-phi.2)*dnorm.2)/(phi.2*(1-phi.2)); #U<-matrix(nrow = nR*nN*nN*2, ncol =2); #U[1:(nN*nR*nD),1]<-u[,1]; #U[(nN*nR*nD=1):(nN*nR*nD*2), 1]<-u[,2]; #U[1:(nN*nR*nD),2]<-0; #U[(nN*nR*nD=1):(nN*nR*nD*2), 2]<-u[,2]; for (k in 1:nN){ temp <- u[u[,3]==k,]; len <- dim(temp)[1]; subset<-matrix(c(temp[,1], temp[,2], rep(0, len), temp[,2]), 2*len, 2); SI = kronecker(subset, subset); si = matrix(c(sum(SI[,1], na.rm = TRUE), sum(SI[,3], na.rm = TRUE), sum(SI[,2], na.rm = TRUE), sum(SI[,4], na.rm = TRUE)), 2, 2); B = B + si; } for (s in (nN + 1):(nN + nD)){ temp <- u[u[,4]==s,]; len <- dim(temp)[1]; subset<-matrix(c(temp[,1], temp[,2], rep(0, len), temp[,2]), 2*len, 2); SI = kronecker(subset, subset); si = matrix(c(sum(SI[,1], na.rm = TRUE), sum(SI[,3], na.rm = TRUE), sum(SI[,2], na.rm = TRUE), sum(SI[,4], na.rm = TRUE)), 2, 2); B = B + si; } for (k in 1:nN){ for (s in (nN + 1):(nN + nD)){ temp <- u[(u[,4]==s & u[,3]==k),]; len <- dim(temp)[1]; subset<-matrix(c(temp[,1], temp[,2], rep(0, len), temp[,2]), 2*len, 2); SI = kronecker(subset, subset); si = matrix(c(sum(SI[,1], na.rm = TRUE), sum(SI[,3], na.rm = TRUE), sum(SI[,2], na.rm = TRUE), sum(SI[,4], na.rm = TRUE)), 2, 2); B = B + si; } } return(B); } MSR<- function(PS, nR, num_pts){ MSR = c(0, 0); for (i in 1:2){ SUM =0; Ave = mean(PS[, i]); for (j in 1:nR){ SUM = SUM + (mean(PS[PS[,4]==j, i]) -Ave)^2; } MSR[i] = num_pts*SUM/(nR-1); } return(MSR) } MST<- function(PS, nR, num_pts){ MST = 0; AveALL = mean(c(PS[,1], PS[,2])); for (i in 1:2){ Ave = mean(PS[, i]); MST = MST + (Ave -AveALL)^2; } MST=MST*nR*num_pts; return(MST) } MSC<- function(PS, nR, num_pts){ #SUM = 0; MSC = c(0, 0); for (i in 1:2){ SUM= 0; Ave = mean(PS[, i]); #ave pseudovalue for treatment i for (j in 1:num_pts){ SUM = SUM + (mean(PS[PS[, 3]==j, i])-Ave)^2; #averaged over the readers for that case (trt i) } MSC[i] = nR*SUM/(num_pts-1); } return(MSC); } MSRC<- function(PS, nR, num_pts){ #SUM =0; MSRC = c(0, 0); for (i in 1:2){ SUM = 0; Ave = mean(PS[,i]); for(j in 1:nR){ R.Ave = mean(PS[PS[,4]==j, i]); for (k in 1:num_pts){ K.Ave = mean(PS[PS[,3]==k, i]); SUM = SUM + (PS[PS[, 3]==k & PS[,4]==j, i]-R.Ave-K.Ave+Ave)^2; } } MSRC[i] = SUM/((nR-1)*(num_pts-1)); } return(MSRC); } MSTR<- function(PS, nR, num_pts){ SUM = 0; AVE = mean(c(PS[,1], PS[,2])); #Y dot dot dot SumIJ = c(0,0); for (i in 1:2){ Ave.i = mean(PS[,i]); #Y i dot dot for(j in 1:nR){ Ave.ij = mean(PS[PS[,4]==j, i]); # Y i j dot Ave.j = mean(c(PS[PS[,4]==j, 1], PS[PS[,4]==j, 2])); #Y dot j dot SumIJ[i] = SumIJ[i] + (Ave.ij - Ave.i - Ave.j + AVE)^2; } } SUM = SumIJ[1] + SumIJ[2]; MSTR=num_pts*SUM/(nR-1); return(MSTR); } MSTC<- function(PS, nR, num_pts){ SUM =0; AVE=mean(c(PS[,1], PS[,2])); #Y dot dot dot for (i in 1:2){ Ave.i = mean(PS[,i]); #Y i dot dot for(j in 1:num_pts){ Ave.ik = mean(PS[PS[,3]==j, i]); # Y i dot k Ave.k = mean(c(PS[PS[,3]==j, 1], PS[PS[,3]==j, 2])); #Y dot dot k SUM = SUM + (Ave.ik - Ave.i - Ave.k + AVE)^2; } } MSTR=nR*SUM/(num_pts-1); return(MSTR); } MSTRC<- function(PS, nR, num_pts){ MSTRC =0; Ydotdotdot = mean(c(PS[,1], PS[,2])); for (i in 1:2){ Yidotdot = mean(PS[,i]); for(j in 1:nR){ Yijdot = mean(PS[PS[,4]==j, i]); Ydotjdot = mean(c(PS[PS[,4]==j, 1], PS[PS[,4]==j, 2])); for (k in 1:num_pts){ Yijk = PS[PS[,3]==k & PS[,4]==j, i]; Ydotjk = mean(c(PS[PS[,3]==k & PS[,4]==j, 1], PS[PS[,3]==k & PS[,4]==j, 2])); Yidotk = mean(PS[PS[,3]==k, i]); Ydotdotk = mean(c(PS[PS[,3]==k, 1], PS[PS[,3]==k, 2])); MSTRC = MSTRC+ (Yijk + -Yijdot - Yidotk - Ydotjk + Yidotdot + Ydotjdot + Ydotdotk - Ydotdotdot)^2; } #k } #j } #i MSTRC=MSTRC/((num_pts-1)*(nR-1)); #return(MSTR); } sigma_est<-function(msR, msC, msRC, nR, nC){ var = (msR+msC-msRC)/(nR*(nC[1]+nC[2])); sig = sqrt(var); return(sig); } ftest<-function(mst, mstr, mstc, mstrc, dof.diff){ f = mst/(mstr + mstc - mstrc); return((f > qf(.95,1,dof.diff))); } ############ FUnctions for bootstrap method############## Bcanon<- function(x,nboot,theta, Y, nR, nC){ n <- length(x) thetahat <- AUC_est_theta(Y, nR);thetahat=thetahat[1]-thetahat[2]; bootsam <- matrix(sample(x,size=n*nboot,replace=T),nrow=nboot) thetastar <- apply(bootsam,1,theta,Y, nR, nC) z0 <- qnorm(sum(thetastar