###################################################################	
### <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))
	}


