SOME EXAMPLE R-code FROM THE BOOK Example 10 Two sample log rank test and more on page 41 ======================================================== library(emplik) data(smallcell) x1 <- smallcell[1:62, 3] x2 <- smallcell[63:121, 3] d1 <- smallcell[1:62, 4] d2 <- smallcell[63:121, 4] temp11 <- Wdataclean3(z=x1, d=rep(1, length(d1))) temp12 <- DnR(x=temp11\$value, d=temp11\$dd, w=temp11\$weight) TIME1 <- temp12\$times RISK1 <- temp12\$n.risk fR1 <- approxfun(x=TIME1, y=RISK1, method="constant", yright=0, rule=2, f=1) temp21 <- Wdataclean3(z=x2, d=rep(1, length(d2))) temp22 <- DnR(x=temp21\$value, d=temp21\$dd, w=temp21\$weight) TIME2 <- temp22\$times RISK2 <- temp22\$n.risk fR2 <- approxfun(x=TIME2, y=RISK2, method="constant", yright=0, rule=2, f=1) flogrank <- function(t){fR1(t)*fR2(t)/(fR1(t)+fR2(t))} emplikHs.test2(x1=x1, d1=d1, x2=x2, d2=d2, theta=0, fun1=flogrank, fun2=flogrank) ## \$`-2LLR` ## [1] 7.247221 ## ...... "All round log rank test" NN <- length(x1) + length(x2) myfun7 <- function(t) {as.numeric((fR1(t)+fR2(t)) > NN*0.6)} fWlogrank <- function(t) {myfun7(t)*flogrank(t)} fBOTH <- function(t) {cbind(flogrank(t), fWlogrank(t))} emplikHs.test2(x1=x1, d1=d1, x2=x2, d2=d2, theta=c(0,0), fun1=fBOTH, fun2=fBOTH) ## \$`-2LLR` ## [1] 19.29306 ## ... 1-pchisq(19.29306, df=2) ## [1] 6.464951e-05 ==================================================================== Example 15, Test and confidence interval for F1(t)-F2(t) on page 79 ==================================================================== P1P2test <- function(theta, x1, d1, x2, d2, t1=365.25, t2=365.25) { a <- theta Jointllr <- function(u, x1=x1, d1=d1, x2=x2, d2=d2, a=a) { temp1 <- el.cen.EM2(x=x1, d=d1, fun=function(x){as.numeric(x >= t1)}, mu = u)\$"-2LLR" temp2 <- el.cen.EM2(x=x2, d=d2, fun=function(x){as.numeric(x >= t2)}, mu = u-a)\$"-2LLR" return(temp1 + temp2) } upBD <- min(0.999999, 1+a) loBD <- max(0.000001, a) temp <- optimize(f = Jointllr, lower = loBD, upper = upBD, x1 = x1, d1 = d1, x2 = x2, d2 = d2, a = a) cstar <- temp\$minimum val <- temp\$objective pvalue <- 1 - pchisq( val, df=1) tempMLE1 <- WKM(x=x1,d=d1) tempMLE2 <- WKM(x=x2,d=d2) P1 <- sum(tempMLE1\$jump[(tempMLE1\$times >= t1)]) P2 <- sum(tempMLE2\$jump[(tempMLE2\$times >= t2)]) MLE <- P1 - P2 list(`-2LLR`=val, Dmle=MLE, ConstrMLE=cstar, Pval=pvalue) } findUL(step=0.1, fun=P1P2test, MLE=0.24, x1=x1, d1=d1, x2=x2, d2=d2) =================================================================== Example 16, Test and confidence interval for F1(t)/F2(t) on page 80 =================================================================== F1F2ratio <- function(theta, x1, d1, x2, d2, t1=365.25, t2=365.25) { a <- theta Jointllr <- function(u, x1=x1, d1=d1, x2=x2, d2=d2, a=a) { temp1 <- el.cen.EM2(x=x1, d=d1, fun=function(x){as.numeric(x <= t1)}, mu = u)\$"-2LLR" temp2 <- el.cen.EM2(x=x2, d=d2, fun=function(x){as.numeric(x <= t2)}, mu=u/a)\$"-2LLR" return(temp1 + temp2) } upBD <- min(0.999999, a) loBD <- max(0.000001) temp <- optimize(f = Jointllr, lower = loBD, upper = upBD, x1 = x1, d1 = d1, x2 = x2, d2 = d2, a = a) cstar <- temp\$minimum val <- temp\$objective pvalue <- 1 - pchisq(val, df=1) tempMLE1 <- WKM(x=x1,d=d1) tempMLE2 <- WKM(x=x2,d=d2) P1 <- sum(tempMLE1\$jump[(tempMLE1\$times <= t1)]) P2 <- sum(tempMLE2\$jump[(tempMLE2\$times <= t2)]) MLE <- P1/P2 list(`-2LLR`=val, Rmle=MLE, ConstrainMLE=cstar, Pval=pvalue) } findUL(step=0.1, fun=F1F2ratio, MLE=0.4187, x1=x1, d1=d1, x2=x2, d2=d2) ================================================================== Example 19, COnfidence Interval in the Cox model, involve both beta and baseline ================================================================== library(survival) library(ELYP) data(smallcell) names(smallcell) ## to see what variable the data include coxph(Surv(survival, indicator)~arm, data=smallcell) ## from here we see beta(MLE) is 0.5337653, with SD=0.203 myy <- smallcell\$survival myd <- smallcell\$indicator myZ <- smallcell\$arm myfun <- function(t){as.numeric(t <= 300)} CoxEL(y=myy, d=myd, Z=myZ, beta=0.5337653, lam=0, fun=myfun)\$logEmpLik ## [1] -500.8937 ## this gives the max that logEmpLik can achieve bvec <- 1:100/110 + 0.06 lvec <- (1:100)*1.3 - 50 myEtafun <- function(beta,theta){(exp(beta) - 1)*theta} LLout1 <- LLout2 <- LLout3 <- matrix(NA, 100,100) for(i in 1:100)for(j in 1:100) { temp <- CoxEL(y=myy, d=myd, Z=myZ, beta=bvec[i], lam=lvec[j],fun=myfun) LLout1[i,j] <- temp\$logEmpLik LLout2[i,j] <- temp\$mu LLout3[i,j] <- myEtafun(beta=bvec[i], theta=temp\$mu) } contour(bvec, lvec, z=LLout1, level=(-500.8937 - 3.84/2), xlab="beta", ylab="lambda", main="logEL and Eta") par(new=TRUE) contour(bvec, lvec, z=LLout3) myIndex <- (LLout1 >= (-500.8937 - 3.84/2)) max( LLout3[myIndex] ) ## [1] 0.240251 min( LLout3[myIndex] ) ## [1] 0.03054699 CoxFindL2(BetaMLE=0.5337653, StepSize=c(0.05, 3), Efun=myEtafun, Hfun=myfun, y=myy, d=myd, Z=myZ) ## ...... ## \$Lower ## [1] 0.03061114 ## ...... CoxFindU2(BetaMLE=0.5337653, StepSize=c(0.05, 3), Efun=myEtafun, Hfun=myfun, y=myy, d=myd, Z=myZ) ## ...... ## \$Upper ## [1] 0.2402764 ## ......