#####################################################################################
#Empirical Likelihood Approach for Two-Group Comparisons Given a U-Statistic Constraint
#####################################################################################

#####################################################################################
#Test H0: Pr(Y>X)=d0 vs Ha: not so
#optimize function in R is used.
#####################################################################################

#There are two variables X and Y with sample sizes of n1 and n2 respectively. 
#And m=n1*n2  
#Obtain estimator for d0 
  indicators=matrix(data=NA, n1, n2)
  for ( u in 1:n1) {
    for ( z in 1:n2) {
    indicators[u,z]=as.numeric(x[u]<y[z])
    }
  }
  indicators_v=c(indicators)
  centerval<-indicators_v-d0
  fn.Lamda<-function(lamval){
    returnval<-(sum(centerval/(m+lamval*centerval)))^2
    return(returnval)
  }

#obtain initial value of lagrange multiplier
  inilamval<-sum(centerval)*m/sum(centerval^2)

#optimize procedure using initial value of lagrange multiplier
  lamval<-optimize(fn.Lamda,interval=c((inilamval-1*abs(inilamval)),(inilamval+1*abs(inilamval))),maximum=FALSE)$minimum

#obtain weights under H0 and make sure weights are within (0,1)
  wijs<-1/(m+lamval*centerval)
  if(length(wijs[wijs<=0])>0){
    xxxx<-matrix(seq((lamval-1*abs(lamval)),(lamval+1*abs(lamval)),0.5),1)
    yyyy<-apply(xxxx,2,fn.Lamda)
    posit<-match(min(yyyy),yyyy)
    inilamval2<-xxxx[posit]
    lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
    wijs<-1/(m+lamval*centerval)
      if(length(wijs[wijs<=0])>0){
      xxxx<-matrix(seq((inilamval-1*abs(inilamval)),(inilamval),0.5),1) 
      yyyy<-apply(xxxx,2,fn.Lamda)
      posit<-match(min(yyyy),yyyy)
      inilamval2<-xxxx[posit]
      lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
      wijs<-1/(m+lamval*centerval)
      }
   }

#Function to compute log(EL)
  fn1.2=function(p){
    sump<-0
    for(j in 1:m){
      sump<-sump+log(p[j])
    }
  return(-sump)
  }

#Compute ELR
  lr=-2*(m*log(n1)+m*log(n2)-fn1.2(wijs))

#Compute variance of estimator for d0 
  auc=d0
  v10=apply(indicators,1,mean)
  s10=sum((v10-auc)^2)/(n1-1)
  v01=apply(indicators,2,mean)
  s01=sum((v01-auc)^2)/(n2-1)
  v=((n2*s10+n1*s01)/(n1+n2)/(n1*n2/(n1+n2))

#Compute test statistics for EL approach
  sum_phi=sum((indicators_v-d0)^2)
  factor1=(v*m^2)/sum_phi
  chis1=lr/factor1

#Compute p-value for hypothesis testing

  pvalue=1-pchisq(chis1, 1)



#####################################################################################
#Test H0: Pr(Y>X)-Pr(Z>W)=0(d0) vs Ha: not so

#optimize function in R is used.
#####################################################################################

#X and W are correlated measurements from same person on biomarker1,   
#Y and Z are correlated measurements from same person on biomarker2,
#with sample sizes of n1 and n2 respectively.
#  
#m=n1*n2

  
#Obtain estimator for Pr(Y>X)-Pr(Z>W)
  indicator_xy=matrix(data=NA, n1, n2)
  for ( a in 1:n1) {
    for ( b in 1:n2) {
      indicator_xy[a,b]=as.numeric(x[a]<y[b])
    }
  }
  indicator_wz=matrix(data=NA, n1, n2)
  for ( p in 1:n1) {
    for ( q in 1:n2) {
    indicator_wz[p,q]=as.numeric(w[p]<z[q])
    }
  }
  indicator_vxy=c(indicator_xy)
  indicator_vwz=c(indicator_wz)
  centerval<-indicator_vxy-indicator_vwz
  fn.Lamda<-function(lamval){
    returnval<-(sum(centerval/(m+lamval*centerval)))^2
    return(returnval)
  }

#compute weights under H0 and make sure weights are within (0,1)

  if (length(centerval[centerval==0])==m) {
    pvalue=1-pchisq(0, 1) } 
  else {
  inilamval<-sum(centerval)*m/sum(centerval^2)
    if (inilamval==0) {
    lamval_0<-optimize(fn.Lamda,interval=c(-100,+100),maximum=FALSE)$minimum
    wijs<-1/(m+lamval_0*centerval)
      if(length(wijs[wijs<=0])>0){
      xxxx<-matrix(seq((lamval_0-1*abs(lamval_0)),(lamval_0+1*abs(lamval_0)),0.1),1)
      yyyy<-apply(xxxx,2,fn.Lamda)
      posit<-match(min(yyyy),yyyy)
      inilamval2<-xxxx[posit]
      lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
      wijs<-1/(m+lamval*centerval)
        if(length(wijs[wijs<=0])>0){
        xxxx<-matrix(seq(-100,0,0.01),1) 
        yyyy<-apply(xxxx,2,fn.Lamda)
        posit<-match(min(yyyy),yyyy)
        inilamval2<-xxxx[posit]
        lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
        wijs<-1/(m+lamval*centerval)
        }
      }
    }
    else {
    lamval<-optimize(fn.Lamda,interval=c((inilamval-1*abs(inilamval)),(inilamval+1*abs(inilamval))),maximum=FALSE)$minimum
    wijs<-1/(m+lamval*centerval)
      if(length(wijs[wijs<=0])>0){
      xxxx<-matrix(seq((lamval-1*abs(lamval)),(lamval+1*abs(lamval)),0.1),1)
      yyyy<-apply(xxxx,2,fn.Lamda)
      posit<-match(min(yyyy),yyyy)
      inilamval2<-xxxx[posit]
      lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
      wijs<-1/(m+lamval*centerval)
        if(length(wijs[wijs<=0])>0){
        xxxx<-matrix(seq((inilamval-1*abs(inilamval)),(inilamval),0.01),1) 
        yyyy<-apply(xxxx,2,fn.Lamda)
        posit<-match(min(yyyy),yyyy)
        inilamval2<-xxxx[posit]
        lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
        wijs<-1/(m+lamval*centerval)
        }
      }
    }
  }

#Function to compute log(EL)
  fn1.2=function(p){
    sump<-0
    for(j in 1:m){
    sump<-sump+log(p[j])
    }
    return(-sump)
  }

#compute ELR
  lr=-2*(m*log(n1)+m*log(n2)-fn1.2(wijs))

#Compute variance of estimator for Pr(Y>X)-Pr(Z>W)
  auc_xy<-mean(indicator_vxy)
  auc_wz<-mean(indicator_vwz)
  v10_x=apply(indicator_xy,1,mean)
  v01_y=apply(indicator_xy,2,mean)
  v10_w=apply(indicator_wz,1,mean)
  v01_z=apply(indicator_wz,2,mean)
  s10_12=sum((v10_x-auc_xy)*(v10_w-auc_wz))/(n1-1)
  s10_11=sum((v10_x-auc_xy)^2)/(n1-1)
  s10_22=sum((v10_w-auc_wz)^2)/(n1-1)
  s01_12=sum((v01_y-auc_xy)*(v01_z-auc_wz))/(n2-1)
  s01_11=sum((v01_y-auc_xy)^2)/(n2-1)
  s01_22=sum((v01_z-auc_wz)^2)/(n2-1)
  s=matrix(c(s10_11,s10_12,s10_12,s10_22),nrow=2,ncol=2,byrow=T)/n1+matrix(c(s01_11,s01_12,s01_12,s01_22),nrow=2,ncol=2,byrow=T)/n2
  v=matrix(c(1,-1),nrow=1) %*% s %*% t(matrix(c(1,-1),nrow=1))

#Compute test statistics for EL approach
  sum_phi=sum((centerval)^2)
  factor1=(v*m^2)/sum_phi
  chis1=lr/factor1

#Compute p-value for hypothesis testing

  pvalue=1-pchisq(chis1, 1)


#####################################################################################
#Test H0: There are no difference between two survival population represented by two 
#         groups of survival censored data
#     Ha: not so.

#optimize function in R is used.
#####################################################################################

#cx and cy are two groups of survival data with censorx and censory indicating if data is
#censored or not. censorx/censory=1 if censored, 0 otherwise 
#n1 and n2 are sample sizes of cx and cy repectively.
#m=n1*n2

  
#Obtain estimator
  indicator1_xy=matrix(data=0, n1, n2)
  indicator2_xy=matrix(data=0, n1, n2)
  for ( a in 1:n1) {
    for ( b in 1:n2) {
    if ((cx[a]>cy[b] & censorx[a]==0 & censory[b]==0)|(cx[a]>=cy[b] & censorx[a]==1 & censory[b]==0)){
    indicator1_xy[a,b]=1} else if ((cx[a]<cy[b] & censorx[a]==0 & censory[b]==0)|(cx[a]<=cy[b] & censorx[a]==0 & censory[b]==1)) {
    indicator2_xy[a,b]=1} 
    }
  }
  indicator1_vxy=c(indicator1_xy)
  indicator2_vxy=c(indicator2_xy)
  centerval<-indicator1_vxy-indicator2_vxy
  fn.Lamda<-function(lamval){
  returnval<-(sum(centerval/(m+lamval*centerval)))^2
  return(returnval)
  }

###compute weights under H0 and make sure weights are within (0,1)
  if (length(centerval[centerval==0])==m) {
  pvalue=1-pchisq(0, 1) } 
  else {
  inilamval<-sum(centerval)*m/sum(centerval^2)
    if (inilamval==0) {
    lamval_0<-optimize(fn.Lamda,interval=c(-100,+100),maximum=FALSE)$minimum
    wijs<-1/(m+lamval_0*centerval)
      if(length(wijs[wijs<=0])>0){
      xxxx<-matrix(seq((lamval_0-1*abs(lamval_0)),(lamval_0+1*abs(lamval_0)),0.1),1)
      yyyy<-apply(xxxx,2,fn.Lamda)
      posit<-match(min(yyyy),yyyy)
      inilamval2<-xxxx[posit]
      lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
      wijs<-1/(m+lamval*centerval)
        if(length(wijs[wijs<=0])>0){
        xxxx<-matrix(seq(-100,0,0.01),1) 
        yyyy<-apply(xxxx,2,fn.Lamda)
        posit<-match(min(yyyy),yyyy)
        inilamval2<-xxxx[posit]
        lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
        wijs<-1/(m+lamval*centerval)
        }
      }
    } 
    else {
    lamval<-optimize(fn.Lamda,interval=c((inilamval-1*abs(inilamval)),(inilamval+1*abs(inilamval))),maximum=FALSE)$minimum
    wijs<-1/(m+lamval*centerval)
      if(length(wijs[wijs<=0])>0){
      xxxx<-matrix(seq((lamval-1*abs(lamval)),(lamval+1*abs(lamval)),0.1),1)
      yyyy<-apply(xxxx,2,fn.Lamda)
      posit<-match(min(yyyy),yyyy)
      inilamval2<-xxxx[posit]
      lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
      wijs<-1/(m+lamval*centerval)
        if(length(wijs[wijs<=0])>0){
          xxxx<-matrix(seq((inilamval-1*abs(inilamval)),(inilamval),0.01),1) 
          yyyy<-apply(xxxx,2,fn.Lamda)
          posit<-match(min(yyyy),yyyy)
          inilamval2<-xxxx[posit]
          lamval<-optimize(fn.Lamda,interval=c((inilamval2-0.02),(inilamval2+0.02)),maximum=FALSE)$minimum
          wijs<-1/(m+lamval*centerval)
        }
      }
    }
  }

#Function to compute log(EL)
  fn1.2=function(p){
  sump<-0
  for(j in 1:m){
    sump<-sump+log(p[j])}
  return(-sump)
  }

#compute ELR
  lr=-2*(m*log(n1)+m*log(n2)-fn1.2(wijs))

#compute variance of estimator
  g=c(cx,cy)
  s=c(censorx,censory)
  uijs=matrix(data=NA, (n1+n2), (n1+n2))
  for ( u in 1:(n1+n2)) {
   for ( z in 1:(n1+n2)) {
   if ((g[u] > g[z] & s[u]==0 & s[z]==0)|(g[u]>=g[z] & s[u]==1 & s[z]==0)) {
   uijs[u,z]=1} else if ((g[u] < g[z] & s[u]==0 & s[z]==0)|(g[u]<=g[z] & s[u]==0 & s[z]==1)) {
   uijs[u,z]=-1} else {
   uijs[u,z]=0}
   }
  }
  ui=apply(uijs,1,sum)
  ustat=sum(ui[1:n1])
  var=n1*n2/((n1+n2)*(n1+n2-1))*sum(ui^2)
  v=var/(m^2)

#Compute test statistics for EL approach
  sum_phi=sum((centerval)^2)
  factor1=(v*m^2)/sum_phi
  chis1=lr/factor1

#Compute p-value for hypothesis testing
  pvalue=1-pchisq(chis1, 1)


#####################################################################################
#Test H0: Pr(Y>X)=d1
#         Pr(Z>W)=d2 
#  vs Ha: not so

#optim function in R is used.
#####################################################################################

#(X, W)are correlated measurements,   
#(Y, Z)are correlated measurements,
#with sample sizes of n1 and n2 respectively.
#m=n1*n2

#Obtain estimator for Pr(Y>X) and Pr(Z>W)
  indicator_1=matrix(data=NA, n1, n2)
  for ( a in 1:n1) {
    for ( b in 1:n2) {
    indicator_1[a,b]=as.numeric(G11[a]<G21[b])  #G11=X, G21=Y
    }
  }
  indicator_2=matrix(data=NA, n1, n2)
  for ( p in 1:n1) {
    for ( q in 1:n2) {
    indicator_2[p,q]=as.numeric(G12[p]<G22[q])  #G12=W, G22=Z
    }
  }
  indicator_v1<-c(indicator_1)
  indicator_v2<-c(indicator_2)
  centerval1<-indicator_v1-d1
  centerval2<-indicator_v2-d2
  fn.Lamda<-function(lamvals){
  lamval1<-lamvals[1]
  lamval2<-lamvals[2]
  returnval1<-(sum(centerval1/(m+lamval1*centerval1+lamval2*centerval2)))^2  
  returnval2<-(sum(centerval2/(m+lamval1*centerval1+lamval2*centerval2)))^2
  return(returnval1+returnval2)
  }

#obtain initial value of vector of lagrange multiplier
  Hval<-cov(cbind(centerval1,centerval2))*(m-1)/m
  inilamval<-c(solve(Hval*m)%*%matrix(c(sum(centerval1),sum(centerval2)),2)*m) 
  awH0<-optim(inilamval,fn.Lamda)
  lam1<-awH0$par[1]
  lam2<-awH0$par[2]

###compute weights under H0 and make sure weights are within (0,1)
  wijs<-1/(m+lam1*centerval1+lam2*centerval2)
  if(length(wijs[wijs<=0])>0){
    inilamval<-c(inilamval[1]+sign(abs(inilamval[1])-abs(lam1))*inilamval[1]*0.3,
    inilamval[2]+sign(abs(inilamval[2])-abs(lam2))*inilamval[2]*0.3)
    awH0<-optim(inilamval,fn.Lamda)
    lam1<-awH0$par[1]
    lam2<-awH0$par[2]
    wijs<-1/(m+lam1*centerval1+lam2*centerval2)
   }
  if(length(wijs[wijs<=0])>0){
  inilamval<-c(inilamval[1]-sign(abs(inilamval[1])-abs(lam1))*inilamval[1]*0.3,
  inilamval[2]-sign(abs(inilamval[2])-abs(lam2))*inilamval[2]*0.3)
  awH0<-optim(inilamval,fn.Lamda)
  lam1<-awH0$par[1]
  lam2<-awH0$par[2]
  wijs<-1/(m+lam1*centerval1+lam2*centerval2)
  }
  if(length(wijs[wijs<=0])>0){
  inilamval<-c(inilamval[1]+sign(abs(inilamval[1])-abs(lam1))*inilamval[1]*0.3,
  inilamval[2]-sign(abs(inilamval[2])-abs(lam2))*inilamval[2]*0.3)
  awH0<-optim(inilamval,fn.Lamda)
  lam1<-awH0$par[1]
  lam2<-awH0$par[2]
  wijs<-1/(m+lam1*centerval1+lam2*centerval2)
  }
  if(length(wijs[wijs<=0])>0){
  inilamval<-c(inilamval[1]-sign(abs(inilamval[1])-abs(lam1))*inilamval[1]*0.3,
  inilamval[2]+sign(abs(inilamval[2])-abs(lam2))*inilamval[2]*0.3)
  awH0<-optim(inilamval,fn.Lamda)
  lam1<-awH0$par[1]
  lam2<-awH0$par[2]
  wijs<-1/(m+lam1*centerval1+lam2*centerval2)
  }
  if(length(wijs[wijs<=0])>0){
  iniseqlam1<-seq(0, (1.5*inilamval[1]), (inilamval[1]/5))
  iniseqlam2<-seq(0, (1.5*inilamval[2]), (inilamval[2]/5))
  itsOK<-0
  for(findi in 1:length(iniseqlam1)){ 
    for(findj in 1:length(iniseqlam2)){
      iniseqvals<-c(iniseqlam1[findi],iniseqlam2[findj])
      awH0<-optim(iniseqvals,fn.Lamda)
      lam1<-awH0$par[1]
      lam2<-awH0$par[2]
      wijs<-1/(m+lam1*centerval1+lam2*centerval2)
      if(length(wijs[wijs<=0])==0){itsOK<-1;break} 
      }
    }
  }
  if(length(wijs[wijs<=0])>0)break;

#Function to compute log(EL)
  fn1.2=function(p){
  sump<-0
    for(j in 1:m){
    sump<-sump+log(p[j])
    }
  return(-sump)
  }

#Compute ELR
  lr=-2*(m*log(n1)+m*log(n2)-fn1.2(wijs))

#compute variance matrix of estimators
  covars<-function(indicator1,indicator2){
  indicator_xy<-indicator1
  indicator_wz<-indicator2
  indicator_vxy=c(indicator1)
  indicator_vwz=c(indicator2)
  auc_xy<-mean(indicator_vxy)
  auc_wz<-mean(indicator_vwz)
  v10_x=apply(indicator_xy,1,mean)
  v01_y=apply(indicator_xy,2,mean)
  v10_w=apply(indicator_wz,1,mean)
  v01_z=apply(indicator_wz,2,mean)
  s10_12=sum((v10_x-auc_xy)*(v10_w-auc_wz))/(n1-1)
  s10_11=sum((v10_x-auc_xy)^2)/(n1-1)
  s10_22=sum((v10_w-auc_wz)^2)/(n1-1)
  s01_12=sum((v01_y-auc_xy)*(v01_z-auc_wz))/(n2-1)
  s01_11=sum((v01_y-auc_xy)^2)/(n2-1)
  s01_22=sum((v01_z-auc_wz)^2)/(n2-1)
  s=matrix(c(s10_11,s10_12,s10_12,s10_22),nrow=2,ncol=2,byrow=T)/n1+matrix(c(s01_11,s01_12,s01_12,s01_22),nrow=2,ncol=2,byrow=T)/n2
  return(s)
  }
  covs12<-covars(indicator_1,indicator_2)
  Sigmaval<-covs12
  Hval<-cov(cbind(centerval1,centerval2))*(m-1)/m
  prefactor<-solve(Hval)%*%Sigmaval*m
  cofs<-eigen(prefactor)$values
  cof1<-cofs[1]
  cof2<-cofs[2]
  
  
# compute p-value
  chi2=cof1*rchisq(50000,1)+cof2*rchisq(50000,1)
  pvalue=(mean(1*((lr)<chi2), na.rm=TRUE))