R/sim_frailty_data.R

In intercure: Cure Rate Estimators for Interval Censored Data

#' Generates a interval censored dataset using frailty cure rate model
#'
#' \code{sim_frailty} 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 beta Two parameters associated with the hazard function.
#' @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).
#' @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_frailty(20)
#' @export
sim_frailty <- function(N, theta = c(-1,1,0),
                        beta = c(0,0.5), A = 5, B = 15,
                        prob = 0.5) {
  u <- stats::runif(N)
  a <- stats::runif(N)
  C <- cbind(A,a * B)
  C <- C[,1] * (C[,1] <= C[,2]) + C[,2] * (C[,1] > C[,2])
  intercept <- 1
  xi1 <- stats::rbinom(N,1,prob)
  xi2 <- stats::rnorm(N)
  cov_theta <- data.frame(intercept, xi1, xi2)
  cov_beta <- data.frame(xi1, xi2)
  eta <- exp(as.vector(theta %*% t(cov_theta)))
  K_vector <- stats::rpois(N, eta / 2)
  U_vector <- K_vector * NA
  for(i in 1:length(K_vector)) {
    if(K_vector[i] == 0) U_vector[i] <- 0
    else{
      U_vector[i] <- 0
      for (j in 1:K_vector[i])
        U_vector[i] <- U_vector[i] + stats::rchisq(1, 2, ncp = 0)
    }
  }
  beta_x <- as.vector(beta %*% t(cov_beta))
  exp_pred_beta <- exp(beta_x)
  num <- -2 * log(1 - u)
  den <- U_vector * exp_pred_beta
  tempos <- ifelse(U_vector != 0, sqrt(num / den), Inf)
  Z <- ifelse(tempos < C, tempos, C)
  delta <- ifelse(tempos < C, 1, 0)
  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)
  return(dados)
}