Bootstrap t confidence interval for censored observations

Question

I have wrote a simulation code for censored observations to find bootstrap-t confidence interval. However, I encountered some problem where my 'btAlpha' and 'btLambda' cannot compute the correct answer hence I cannot go to the next step which is to calculate the total error probabilities.

This is my code :

#BOOTSTRAP-T (20%)

library(survival)
n <- 100
N <- 1000
alpha <- 1
lambda <- 0.5
alphaHat <- NULL
lambdaHat <- NULL
cp <- NULL
btAlpha <- matrix (NA, nrow=N, ncol=2)
btLambda <- matrix (NA, nrow=N, ncol=2)

for (i in 1:1000) {
  u <- runif(n)
  c1 <- rexp(n, 0.1)
  t1 <- -(log(1 - u^(1/alpha))/lambda)
  t <- pmin(t1, c1)
  ci <- 1*(t1 < c1) #censored data
  cp[i] <- length(ci[ci == 0])/n #censoring proportion
  
  #FUNCTION TO CALL OUT
  estBoot < -function(data, j) {
    dat <- data [j, ]
    data0 <- dat[which(dat$ci == 0), ] # right censored data
    data1 <- dat[which(dat$ci == 1), ] # uncensored data
    dat
    
    #MAXIMUM LIKELIHOOD ESTIMATION  
    library(maxLik)
    LLF <- function(para) {
      alpha <- para[1]
      lambda <- para[2]
      a <- sum(log(alpha*lambda*(1 - exp(-lambda*data1$t1))^(alpha - 1)*
                     exp(-lambda*data1$t1)))
      b <- sum(log(1 - (1 - exp(-lambda*data0$t1)^(alpha))))
      l <- a + b
      return(l)
    }
    mle <- maxLik(LLF, start=c(alpha=1, lambda=0.5))
    
    alphaHat <- mle$estimate[1]
    lambdaHat <- mle$estimate[2]
    
    observedDi <- solve(-mle$hessian)
    
    return(c(alphaHat, lambdaHat, observedDi[1, 1], observedDi[2, 2]))
  }
  
  library(boot)
  
  bt <- boot(dat, estBoot, R=1000)
  
  bootAlphaHat <- bt$t[, 1] #t is from bootstrap
  bootAlphaHat0 <- bt$t0[1] #t0 is from original set
  seAlphaHat <- sqrt(bt$t[, 2])
  seAlphaHat0 <- sqrt(bt$t0[2]) #same as 'original' in bt
  zAlpha <- (bootAlphaHat - bootAlphaHat0)/seAlphaHat
  kAlpha <- zAlpha[order(zAlpha)]
  ciAlpha <- c(kAlpha[25], kAlpha[975])
  btAlpha[i, ] <- rev(bootAlphaHat0 - ciAlpha*seAlphaHat0)
  
  bootLambdaHat <- bt$t[, 2]
  bootLambdaHat0 <- bt$t0[2]
  seLambdaHat <- sqrt(bt$t[, 4])
  seLambdaHat0 <- sqrt(bt$t0[4])
  zLambda <- (bootLambdaHat - bootLambdaHat0)/seLambdaHat
  kLambda <- zLambda[order(zLambda)]
  ciLambda <- c(kLambda[25], kLambda[975])
  btLambda[i, ] <- rev(bootLambdaHat0 - ciLambda*seLambdaHat0)
}

leftAlpha <- sum(btAlpha[, 1] > alpha)/N
rightAlpha <- sum(btAlpha[, 2] < alpha)/N
totalEAlpha <- leftAlpha + rightAlpha

leftLambda <- sum(btLambda[, 1] > lambda)/N
rightLambda <- sum(btLambda[, 2] < lambda)/N
totalELambda <- leftLambda + rightLambda

#alpha=0.05
sealphaHat <- sqrt(0.05*(1 - 0.05)/N)

antiAlpha <- totalEAlpha > (0.05 + 2.58*sealphaHat)
conAlpha <- totalEAlpha < (0.05 - 2.58*sealphaHat )
asymAlpha <- (max(leftAlpha, rightAlpha)/min(leftAlpha, rightAlpha)) > 1.5

antiLambda <- totalELambda > (0.05 + 2.58 *sealphaHat)
conLambda <- totalELambda < (0.05 - 2.58 *sealphaHat)
asymLambda <- (max(leftLambda, rightLambda)/min(leftLambda, rightLambda)) > 1.5

anti <- antiAlpha + antiLambda
con <- conAlpha + conLambda
asym <- asymAlpha + asymLambda

My 'btAlpha[i,]' and 'btLambda[i,]' is two matrix data frame and only computed NA values hence I cannot calculate the next step which is total error probabilities etc. It should be simulated 1000 values through specified formula but I didnt get the desired output. I have tried to run this without using loops and same problems encountered. Do you guys have any idea? I could really use and truly appreciate your help.

Bootstrap t confidence interval for censored observations

0 Answers