For all these analyses, the data should be in the same format. There are 5 columns. The first has the ratings or scores from the first treatment, the second is the second treatment. The third column is disease status, fourht is case, and 5th is reader Disease status, case, and reader should all be sorted low to high (eg, start with ratings from the first reader, non-diseased patients, case 1 to number non-diseased etc). Also, cases must be numbered starting from 1. The columns are: Treat1 Treat2 Disease Case Reader #######Mammogram data example########## #read in the data and set number of readers, patients jiang<-Jiang[,c(1,2,3,5,6)] nC<-c(58, 46) nR<-10; nN<-58; nD<-46 #############DBM analysis############ Y <- jiang PS <- pseudos3(Y, nC[1], nC[2], nR); mst<- MST(PS, nR, nC[1]+nC[2]); msr<- MSR(PS, nR, nC[1]+nC[2]); msc<- MSC(PS, nR, nC[1]+nC[2]); msrc<- MSRC(PS, nR, nC[1]+nC[2]); mstr<- MSTR(PS, nR, nC[1]+nC[2]); mstc<- MSTC(PS, nR, nC[1]+nC[2]); mstrc<- MSTRC(PS, nR, nC[1]+nC[2]); est<- c(mean(PS[,1]), mean(PS[,2]), mean(PS[,1])-mean(PS[,2])); #difference CI dof.diff<-DDFH(mstr, mstc, mstrc, nR); sigma.diff<-sqrt(2)*sigma_est(mstr, mstc, mstrc, nR, nC); CI.radius.hillis <- abs(sigma.diff*qt(.025, dof.diff)); CI.diff.hillis <- c(est[3] - CI.radius.hillis, est[3] + CI.radius.hillis); # ftest: fCoverage1[i] <- ftest(mst, mstr, mstc, mstrc, dof.diff); f = mst/(mstr + mstc - mstrc); ftest =(f > qf(.95,1,dof.diff)); > #Confidence intervals for treatment means: dof1<- DDFH(msr[1], msc[1], msrc[1], nR); sig1<-sqrt((msr[1]+msc[1]-msrc[1])/(nR*(nC[1]+nC[2]))); num1<- sig1*qt(.975, dof1); CI1 <- c(est[1] - num1, est[1] + num1); dof2<- DDFH(msr[2], msc[2], msrc[2], nR); sig2<-sqrt((msr[2]+msc[2]-msrc[2])/(nR*(nC[1]+nC[2]))); num2<- sig2*qt(.975, dof2); CI2 <- c(est[2] - num2, est[2] + num2); ###########Results:############## #Estimates (treatment 1 AUC, treatment2 AUC, treatment1-treatment2) est #CI's CI.diff.hillis #difference CI1 #treatment 1 CI2 #treatment 2 #F statistic and test f #f statistic ftest #reject; 1 = yes, 0 = no 1-pf(f,1,dof.diff) # p-value ###########OR analysis############### Y<-jiang BIdata<-gen.data(Y, nC[1], nC[2], nR) x<-covariance_matrix(BIdata, nR, nC[1], nC[2]) r2sig<-r2sigmac_hat(x, nR) r3sig<-r3sigmac_hat(x, nR) PS<-pseudos3(Y, nC[1], nC[2], nR); ORmst<-OR_MST(PS, nR, nC[1]+nC[2]); ORmstr<-OR_MSTR(PS, nR, nC[1]+nC[2]); ORmsrBOTH<-OR_MSRboth(PS, nR) ORmsrSINGLE<-OR_MSRsingle(PS, nR) dofDIFF<-DDFH_OR(ORmstr, nR, .5*(r2sig[1]+r2sig[2]), r3sig) dofSINGLE<-DDFH_single(ORmsrSINGLE, r2sig, nR); trtCIs<-single_trt_CI_radius(ORmsrSINGLE, r2sig, nR, dofSINGLE) diffCIradius<-trt_difference_CI_radius(ORmstr, .5*(r2sig[1]+r2sig[2]), r3sig, nR, dofDIFF) ORests<-c(mean(BIdata[,1]), mean(BIdata[,2]), mean(BIdata[,1]) - mean(BIdata[,2])); if(is.na(diffCIradius)){diffCIradius<-0}; if(is.na(trtCIs[1])){trtCIs[1]<-0}; if(is.na(trtCIs[2])){trtCIs[2]<-0}; OR_f = ORmst/(ORmstr + nR*max((.5*(r2sig[1]+r2sig[2]) - r3sig),0)) OR_Fstat = (OR_f > qf(.95,1,dofDIFF)); ##############Results################# #Estimates ORests #first entry is treatment 1, second is treatment 2, third is 1-2 #CI's #Treatment 1: c(ORests[1]-trtCIs[1], ORests[1]+trtCIs[1]) #Treatment 2: c(ORests[2]-trtCIs[2], ORests[2]+trtCIs[2]) #Difference: c(ORests[3]-diffCIradius, ORests[3]+diffCIradius) #F OR_f #statistic OR_Fstat # reject at 95% level; 1 = yes, 0 = no #p-value 1-pf(OR_f,1,dofDIFF) ###############Marginal model############# BIdata<-gen.data(Y, nC[1], nC[2], nR) BIdata[which(BIdata==0.5)]<-0 ####This line is only for continuous data, and should not make a difference if there are no ties response<-c(BIdata[,1], BIdata[,2]); len <-length(BIdata[,1]); predictor.1<-c(rep(0, len), rep(1, len)); lin.mod.1<-glm(response~predictor.1, family = binomial(link = "probit")); beta0 = lin.mod.1$coefficients[1]; beta1 = lin.mod.1$coefficients[2]; beta = c(beta0, beta1, beta0+beta1); xbeta <- c(beta0, beta0 + beta1); AUC = c(pnorm(beta0), pnorm(beta0 + beta1), 1); AUC[3] = AUC[1]-AUC[2]; se <- matrix(c(0, 0, 0),3,1); se_AUC<-matrix(c(0, 0, 0),3,1); temp1 = dnorm(beta0); temp2 = dnorm(beta0 + beta1); #compute A A = getINF(BIdata, xbeta); A = solve(A); # compute B B = getB(BIdata, xbeta, nN, nD, nR); V = A%*%B%*%t(A); D <- matrix(c(1,1), 1, 2); V_beta2 = D%*%V%*%t(D); D <- matrix(c(temp1, temp2, temp1 - temp2, 0, temp2, (-1)*temp2), 3, 2); V_AUC <- D%*%V%*%t(D); se = c(sqrt(abs(diag(V))), sqrt(abs(V_beta2))); se_AUC = sqrt(abs(diag(V_AUC))); AUC_low = c(pnorm(beta[1] - qnorm(0.975)*se[1]), pnorm(beta[3] - qnorm(0.975)*se[3]), (AUC[3] - qnorm(0.975)*se_AUC[3])); AUC_up = c(pnorm(beta[1] + qnorm(0.975)*se[1]), pnorm(beta[3] + qnorm(0.975)*se[3]), (AUC[3] + qnorm(0.975)*se_AUC[3])); ##Results############### AUC #estimates AUC_low #lower 95% CI AUC_up #upper 95% CI ####################Bootstrap################ Y<-jiang; nC<-c(58,46); bootnumber<-2000; n<-(nC[1]+nC[2]); CI<-Bcanon(1:n, bootnumber, theta, Y, nR, nC) ##CI for Treatment 1 AUC - treatment 2 AUC ## This outputs both the BCA and percentile CI's