R/sim_bch_data.R

In intercure: Cure Rate Estimators for Interval Censored Data

#' Generates a interval censored dataset using promotion time cure rate model
#'
#' \code{sim_frailty_data} returns a dataset generated by the cure rate frailty
#' model.
#' @param N Size of the sample to be generated.
#' @param theta Three parameters associated with the cure linear predictor.
#' @param A A positive number representing a fixed right censoring.
#' @param B A positive number which multiplies an uniform random variable,
#'   defining another right censoring case.
#' @param prob Probability that individual presents treatment T1 (baseline is
#'   T0).
#' @param lambda Rate parameter for the exponential distributed latent
#'   variables.
#' @return A generated dataset with columns: \code{Z}, the actual event time;
#'   \code{L}, the leftmost limit of the censored interval; \code{R}, the
#'   rightmost limit of the censored interval; \code{delta}, the failure
#'   indicator; \code{xi1}, the treatment covariate assuming 1 with probability
#'   \code{prob} and 0 otherwise; \code{xi2}, second variable generated by a
#'   standard normal distribution.
#' @examples
#' sim_bch(20)
#' @export
sim_bch <- function(N, theta = c(1, 0.5, 0),
                    lambda = 1, A = 5, B = 15, prob = 0.5){
  intercept <- 1
  xi1 <- stats::rbinom(N, 1, prob)
  xi2 <- stats::rnorm(N)
  X <- data.frame(intercept, xi1, xi2)
  lambda_pois <- as.numeric(exp(theta %*% t(X)))
  N_L <- stats::rpois(N, lambda_pois)
  a <- stats::runif(N)
  C <- cbind(A,a * B)
  C <- C[,1] * (C[,1] <= C[,2]) + C[,2] * (C[,1] > C[,2])
  T <- c(1:N) * NA
  for (i in 1:N) {
    T[i] <- ifelse(N_L[i] > 0, min(stats::rexp(N_L[i],
                                               lambda)), Inf)
  }
  delta <- ifelse(T <= C,1,0)
  Z <- ifelse(T <= C, T, C)
  L <- R <- Z * NA
  for (i in 1:N){
    if (delta[i] == 0){
      L[i] <- Z[i]
      R[i] <- Inf
    }
    else{
      L[i] <- 0
      add  <- stats::runif(1, 0.1, 0.5)
      R[i] <- add
      check <- (L[i] <= Z[i] & Z[i] < R[i])
      while(!check){
        L[i] <- L[i] + add
        add <- stats::runif(1, 0.1, 0.5)
        R[i] <- R[i] + add
        check <- (L[i] <= Z[i] & Z[i] < R[i])
      }
    }
  }
  dados <- data.frame(Z, L, R, delta, xi1, xi2, N_L)
  return(dados)
}