Search code examples
rsurvival-analysis

Finding the empirical likelihood ratio

I have a homework about an article and the result must be in interval (0.10024, 1.0917).

As a concrete example we took the data of remission times for solid tumor patients (n=10), which is a slightly modified (breaktie) version of Statistical Methods for Survival Data Analysis, Elisa T. Lee, 1992, Example 4.2:

3, 6.5, 6.51, 10, 12, 15, 8.4+, 4+,5.7+, and 10+.

Suppose we are interested in getting a 95% confidence interval for the cumulative hazard at the time t = 9.8, ∆o (9.8). Hence θo=∆o(9.8). In this case the function g is an indicator function: g(t)=I[t9.8].

The 95% confidence interval using the empirical likelihood ratio, -2logALR, for ∆o (9.8) is (0.10024, 1.0917)
please help me to get result above. Thank you. 

remissiontime<-(3,4,5.7,6.5,6.51,8.4,10,10,12,15) 
status <- c(1,0,0,1,1,0,1,0,1,1)

and my code is (actually I am not sure about this code)

library(survival)
library(emplik)
x1 = c(3,4,5.7,6.5,6.51,8.4,10,10,12,15) 
d1 = c(1,0,0,1,1,0,1,0,1,1)

KM0 <- survfit(Surv(x1,d1) ~ 1,  type="kaplan-meier", conf.type="log")
summary(KM0)

myfun <-function(t){as.numeric(t <=9.8)}

emplikH1.test(x=x1,d=d1,theta=-log(0.643),fun=myfun)


myULfun <-function(theta,x,d){

  emplikH1.test(x=x1,d=d1,theta=theta,fun=function(t){as.numeric(t <= 9.8)})}

findUL(fun=myULfun,MLE =-log(0.643),x=x1,d=d1)
Solution

  • The last line of the code you presented

    findUL(fun=myULfun,MLE =-log(0.643),x=x1,d=d1)

    throws an error:

    Error in emplikH1.test(x = x1, d = d1, theta = theta, fun = function(t) { : given theta is too far away from theta0

    Changing MLE argument of findUL from -log(0.643) to 1 cure the issue. Please see the below:

    findUL(fun = myULfun, MLE = 1, x = x1, d = d1)

    Output, which coincides with the one required:

    $`Low`
[1] 0.1002516

$Up
[1] 1.09165

$FstepL
[1] 6.103516e-05

$FstepU
[1] 6.103516e-05

$Lvalue
[1] 3.839313

$Uvalue
[1] 3.839971