#' Generates a interval censored dataset using the promotion time cure rate
#' model
#'
#' \code{sim_bch} returns a dataset generated by the promotion time cure rate
#' 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 exponential random variable
#' with mean 1, 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. \code{N_latent} is the number of remaining
#' tumorous cells post treatment from the promotion time motivation. Cured
#' individuals are represented by zero on this variable.
#' @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::rexp(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_latent = N_L)
return(dados)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.