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