R/rng.R

In MultiFrailty: Shared Frailty Regression Models with Inverse Gaussian, Generalized Lindley, and Gamma Frailty Distributions

#' Random Data Generation under Eight Censoring Schemes
#'
#' Generates survival times and censoring indicators for shared frailty regression models
#' across all 10 baseline-frailty combinations under eight distinct censoring mechanisms.
#'
#' @param n Number of observations to generate.
#' @param baseline Baseline hazard distribution (\code{"weibull"} or \code{"gw"}).
#' @param bpar Baseline parameter vector.
#' @param frailty Frailty distribution (\code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}).
#' @param fpar Frailty parameter vector.
#' @param x Matrix of covariates (n x p). Default is 0-column matrix.
#' @param beta Regression coefficient vector matching columns of \code{x}.
#' @param cen_type Censoring mechanism: \code{"none"}, \code{"right"}, \code{"left"}, \code{"interval"},
#'   \code{"type1"}, \code{"type2"}, \code{"progressive"}, or \code{"progressive_type1"}.
#' @param cen_rate Exponential rate for right-censoring time generation. Default is 0.2.
#' @param left_threshold Threshold for left censoring. Default is 20th percentile.
#' @param int_width Width of censoring window for interval censoring. Default is 20\% of mean time.
#' @param cen_time Fixed cutoff time for Type-I censoring. Default is 70th percentile.
#' @param r_failures Target number of failures for Type-II censoring. Default is \code{floor(0.7 * n)}.
#' @param prog_scheme Vector of progressive removal counts for progressive censoring.
#' @param prog_times Inspection time points for progressive Type-I censoring.
#'
#' @return A data frame containing generated \code{time}, \code{time2} (for interval), \code{status}, and covariates.
#'
#' @references
#' Hougaard, P. (1984). Life table methods for heterogeneous populations: distributions of frailties. Biometrika, 71(1), 75-83.
#'
#' Pandey, A., Hanagal, D. D., & Tyagi, S. (2022). Shared Frailty Models Based on Cancer Data. International Journal of Statistics and Reliability Engineering, 9(3), 461-474.
#'
#' Pandey, A., & Tyagi, S. (2021). Comparison of Multiplicative Frailty Models Under Weibull Baseline Distribution. Lobachevskii Journal of Mathematics, 42(13), 3184-3195.
#'
#' @export
#' @examples
#' set.seed(123)
#' dat <- r_frailty(n = 100, baseline = "weibull", bpar = c(2, 1.5),
#'                  frailty = "gl1", fpar = c(1.2, 0.5),
#'                  cen_type = "right", cen_rate = 0.1)
#' head(dat)
r_frailty <- function(n, baseline = c("weibull", "gw"), bpar,
                      frailty = c("none", "gamma", "ig", "gl1", "gl2"), fpar = numeric(0),
                      x = matrix(nrow = n, ncol = 0), beta = numeric(0),
                      cen_type = c("none", "right", "left", "interval", "type1", "type2",
                                   "progressive", "progressive_type1"),
                      cen_rate = 0.2, left_threshold = NULL, int_width = NULL,
                      cen_time = NULL, r_failures = NULL, prog_scheme = NULL, prog_times = NULL) {
  baseline <- match.arg(baseline)
  frailty <- match.arg(frailty)
  cen_type <- match.arg(cen_type)

  if (n <= 0 || n != as.integer(n)) stop("'n' must be a positive integer.")
  if (!is.matrix(x)) x <- as.matrix(x)
  if (nrow(x) != n && nrow(x) != 0) stop("Number of rows in 'x' must equal 'n'.")

  n_cov <- ncol(x)
  if (n_cov > 0) {
    if (length(beta) != n_cov) stop("Length of 'beta' must match number of columns in 'x'.")
    rho <- as.vector(exp(x %*% beta))
  } else {
    rho <- rep(1.0, n)
  }

  # 1. Draw frailties W
  if (frailty == "none") {
    w <- rep(1.0, n)
  } else if (frailty == "gamma") {
    theta <- fpar[1]
    w <- stats::rgamma(n, shape = 1 / theta, scale = theta)
  } else if (frailty == "ig") {
    w <- r_ig(n, eta = fpar[1])
  } else if (frailty == "gl1") {
    w <- r_gl1(n, eta = fpar[1], epsilon = fpar[2])
  } else if (frailty == "gl2") {
    w <- r_gl2(n, theta = fpar[1], mu = fpar[2])
  }

  # 2. Invert survival function to draw true event times T
  U <- stats::runif(n)
  v_val <- -log(U) / (w * rho) # Phi0(t) = v_val

  if (baseline == "weibull") {
    lambda <- bpar[1]
    gamma <- bpar[2]
    # Phi0(t) = (t/lambda)^gamma = v_val => t = lambda * (v_val)^(1/gamma)
    T_true <- lambda * (v_val^(1 / gamma))
  } else if (baseline == "gw") {
    delta <- bpar[1]
    zeta <- bpar[2]
    xi <- bpar[3]
    # Phi0(t) = -log1mexp(v_inner_gw) = v_val
    v_inner <- -(1 / zeta) * log1mexp(v_val)
    T_true <- ((-1 / delta) * log1mexp(v_inner))^(1 / xi)
  }

  # 3. Apply requested censoring scheme
  time <- T_true
  time2 <- rep(NA_real_, n)
  status <- rep(1L, n)
  prog_cen <- rep(0L, n)

  if (cen_type == "none") {
    # exact events
  } else if (cen_type == "right") {
    C <- stats::rexp(n, rate = cen_rate)
    status <- ifelse(T_true <= C, 1L, 0L)
    time <- pmin(T_true, C)
  } else if (cen_type == "left") {
    if (is.null(left_threshold)) left_threshold <- stats::quantile(T_true, 0.20)
    status <- ifelse(T_true <= left_threshold, 2L, 1L)
    time <- ifelse(status == 2L, left_threshold, T_true)
  } else if (cen_type == "interval") {
    if (is.null(int_width)) int_width <- mean(T_true) * 0.2
    is_cen <- stats::runif(n) < 0.5
    status <- ifelse(is_cen, 3L, 1L)
    time <- ifelse(status == 3L, pmax(0.001, T_true - int_width / 2), T_true)
    time2 <- ifelse(status == 3L, T_true + int_width / 2, T_true)
  } else if (cen_type == "type1") {
    if (is.null(cen_time)) cen_time <- stats::quantile(T_true, 0.70)
    status <- ifelse(T_true <= cen_time, 1L, 0L)
    time <- pmin(T_true, cen_time)
  } else if (cen_type == "type2") {
    if (is.null(r_failures)) r_failures <- floor(0.7 * n)
    r_failures <- pmin(pmax(1, r_failures), n)
    T_sort <- sort(T_true)
    c_cutoff <- T_sort[r_failures]
    status <- ifelse(T_true <= c_cutoff, 1L, 0L)
    time <- pmin(T_true, c_cutoff)
  } else if (cen_type == "progressive") {
    if (is.null(prog_scheme)) {
      m <- floor(0.7 * n)
      prog_scheme <- rep(0L, m)
      rem <- n - m
      if (rem > 0) prog_scheme[m] <- rem
    }
    # Progressive censoring simulation
    ord <- order(T_true)
    time <- T_true[ord]
    status <- rep(1L, n)
    prog_cen <- rep(0L, n)
    m <- length(prog_scheme)
    if (m < n) {
      status[(m + 1):n] <- 0L
      prog_cen[1:m] <- prog_scheme
    }
  } else if (cen_type == "progressive_type1") {
    if (is.null(prog_times)) prog_times <- stats::quantile(T_true, c(0.25, 0.50, 0.75))
    k_times <- length(prog_times)
    status <- rep(1L, n)
    for (j in 1:k_times) {
      t_inspect <- prog_times[j]
      idx <- (T_true > t_inspect) & (status == 1L)
      if (any(idx)) {
        # remove a fraction at inspection time
        rem_idx <- which(idx)[stats::runif(sum(idx)) < 0.2]
        if (length(rem_idx) > 0) {
          status[rem_idx] <- 0L
          time[rem_idx] <- t_inspect
          prog_cen[rem_idx] <- 1L
        }
      }
    }
  }

  res_df <- data.frame(
    time = time,
    time2 = time2,
    status = status,
    prog_cen = prog_cen
  )

  if (n_cov > 0) {
    colnames(x) <- if (is.null(colnames(x))) paste0("X", 1:n_cov) else colnames(x)
    res_df <- cbind(res_df, as.data.frame(x))
  }

  res_df
}