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
## ......