################################################################### ### <June 2001>Hemant Ishwaran ################################################################### ### Test fixed linear constrasts for dependent U-statistics ### Based on DeLong, DeLong, and Pearson (1988). Biometrics. ### ### Works for multivariate ROC data; ie the setting where ### each individual is ranked by M different raters with ### no missing data [for case of missing data a bootstrap ### procedure is also available; please contact me] ### ### ### Function returns an object containing: ### ------------------------------------- ### (1) area: ROC areas for each of M raters (Mx1) ### (2) test.statistic: test statistic value (lx1) ### (3) chi.square: chisquare value for overall test, ### degree of freedom, ### and p-value ### (4) Var: variance matrix for test statistic ### ### where contrast is a l x M contrast matrix representing ### the desired test (rank is <= M) ### ### Example Degree of suspicion of stage C/D prostate ### ------- cancer on ordinal scale {1,2,...100}. Each ### patient ranked by 5 different radiologists. ### ### -- true disease status for each patient (99 values): ### ### group <- c( ### 1,0,1,0,1,1,0,1,1,1,1,0,1,1,0,0,1,1,1,0,1,1,0,0,1,0,1,0,1,0,0,1,0,1,0, ### 0,0,1,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,1,1,0,1,1,0,1,1, ### 1,1,0,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,0,1,0,1,1,1,1,1,1,1) ### ### -- ordinal data for five readers (99x5): ### ### ordinal <- matrix(c( ### 99,90,100,100,100,5,5,10,25,49,40,25,95,75,65,20,20,20,15,50,10,35,50, ### 20,55,30,20,20,30,55,5,10,10,10,20,50,80,95,80,95,5,15,20,70,60,25, ### 35,60,100,80,50,5,35,50,30,5,40,10,15,10,60,75,90,75,95,10,30,10,60, ### 55,5,40,10,15,5,5,50,5,65,25,15,40,70,35,60,15,20,25,50,30,5,25,50, ### 20,10,5,5,10,10,10,40,95,100,100,95,60,30,60,30,51,5,0,5,5,1,5,10,25, ### 10,20,20,10,60,60,45,5,90,10,25,10,60,25,50,25,50,5,20,10,30,10,15, ### 60,85,25,45,5,20,40,60,30,5,25,10,60,55,30,10,10,30,10,25,20,50,25, ### 20,5,50,70,50,75,5,20,10,25,5,5,25,10,25,10,5,25,10,25,10,60,35,50, ### 35,65,40,20,20,20,30,5,80,10,40,20,5,25,20,25,5,10,90,50,20,65,10,15, ### 10,10,20,5,5,10,75,15,0,5,10,20,1,40,80,85,100,95,5,40,10,25,10,50, ### 25,25,25,40,5,5,25,0,15,5,5,25,25,10,5,60,20,60,35,70,80,100,100,45, ### 5,40,25,50,30,50,100,95,50,85,0,10,10,50,5,5,20,10,35,20,30,25,50,15, ### 55,10,30,30,75,10,20,20,10,75,10,60,80,60,100,5,5,25,10,30,25,10,35, ### 30,20,5,5,50,50,60,35,10,80,10,50,51,5,10,10,60,10,70,100,100,100,90, ### 10,20,10,70,5,15,20,25,100,90,40,60,25,100,55,10,10,10,10,55,80,75, ### 95,100,100,40,20,15,55,5,5,60,5,100,40,60,20,90,100,65,5,25,25,20,35, ### 60,50,30,40,50,5,10,50,10,65,10,25,10,20,25,10,5,15,50,10,50,60,25, ### 80,51,15,89,85,50,75,5,10,60,30,60,60,60,100,70,75,20,65,30,100,5,70, ### 80,85,70,85,10,25,10,100,35,5,20,10,10,10,30,70,50,100,30,10,20,5,50, ### 65,50,70,10,30,35,10,20,30,50,10,60,60,10,40,35,5,60,30,10,55,40,50, ### 20,100,65,50,50,95,70,90,80,80,100,100,100,5,20,40,15,35,15,65,45,30, ### 45,25,50,70,75,75),ncol=5,byrow=T) ### ### -- compute ROC areas and standard errors ### ### contrast.mat <- diag(5) ### prostate.out <- u.contrast(group,ordinal,contrast.mat) ### readers.area <- prostate.out$test ### readers.se <- sqrt(diag(prostate.out$Var)) ### ### -- test if readers 1,2 and 3 are different than 4 and 5 ### -- also test if readers 4 and 5 are different from one another ### ### contrast.mat2 <- rbind(c(1/3,1/3,1/3,-1/2,-1/2),c(0,0,0,1/2,-1/2)) ### prostate.out2 <- u.contrast(group,ordinal,contrast.mat2) ### prostate.test <- prostate.out2$test ### prostate.pv <- prostate.out2$chi[3] ### ### -- from prostate.pv the overall p.value is 0.036 ### -- signifying that we reject the null (no differences) ### -- also look at prostate.test for individual test statistics ### ### ### ------------------------------------------------------------- ### Hemant Ishwaran ishwaran@bio.ri.ccf.org ### Dept of Biostatistics/Wb4 phone: 216-444-9932 ### Cleveland Clinic Foundation fax: 216-445-2781 ### 9500 Euclid Avenue ### Cleveland, OH 44195 ### ### http://www.bio.ri.ccf.org/Resume/Pages/Ishwaran ################################################################### ### Function which computes U-statistics. Assign weight of ### 0.5 for ties. ustat.con_function(neg,pos) { n.neg_length(neg) n.pos_length(pos) if (n.neg==1 | n.pos==1){ #labour saving device sum((neg<pos)+0.5*(neg==pos))/(n.neg*n.pos) } else{ ranked.cut_sort(unique(c(neg,pos))) n.ranked_length(ranked.cut) n.atoms_n.ranked+1 fp_rep(0,n.atoms) tp_rep(0,n.atoms) for (i in 1:n.ranked){ fp[i]_sum(neg>=ranked.cut[i])/n.neg tp[i]_sum(pos>=ranked.cut[i])/n.pos } fp.diff_-diff(fp) tp.mean_(tp[1:(n.atoms-1)] + tp[2:(n.atoms)])/2 sum(fp.diff*tp.mean) } } u.contrast_function(group,ordinal,contrast) { M_dim(ordinal)[2] group.uniq_unique(sort(group)) contrast_as.matrix(contrast) n.neg_length(group[group==group.uniq[1]]) n.pos_length(group)-n.neg Y_ordinal[(group==group.uniq[1]),c(1:M)] # Nondiseased group X_ordinal[(group==group.uniq[2]),c(1:M)] # Diseased group ### Compute all areas (there are M of these) theta_rep(0,M) for (i in 1:M){ theta[i]_ustat.con(Y[,i],X[,i]) } ### Compute X,Y variance matrices V01_matrix(0,n.neg,M,byrow=T) V10_matrix(0,n.pos,M,byrow=T) for (row in 1:n.neg){ for (col in 1:M){ V01[row,col]_ustat.con(Y[row,col],X[,col]) } } for (row in 1:n.pos){ for (col in 1:M){ V10[row,col]_ustat.con(Y[,col],X[row,col]) } } S01_(t(V01)%*%V01-n.neg*outer(theta,theta))/(n.neg-1) S10_(t(V10)%*%V10-n.pos*outer(theta,theta))/(n.pos-1) Svar_S10/n.pos+S01/n.neg Var_contrast%*%Svar%*%t(contrast) test.L_contrast%*%theta #Test value test_t(test.L)%*%solve(Var)%*%test.L #Chisquare value deg.f_min(dim(contrast)) chi.p_1-pchisq(test,deg.f) ### Return promised object: return(list(area=theta,test.statistic=c(test.L), chi.square=c(test,deg.f,chi.p), Var=Var)) }