R/gen_data.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Generate Data under Type-I Hybrid Censoring Scheme
#'
#' @param pdf Probability density function.
#' @param cdf Cumulative distribution function.
#' @param lower Lower bound of the support of the distribution.
#' @param upper Upper bound of the support of the distribution.
#' @param n Total sample size.
#' @param r Minimum required number of failures.
#' @param T_star Fixed termination time.
#' @param seed Optional integer random seed for reproducibility.
#'
#' @return A list containing:
#' \item{observed_times}{Vector of observed failure/censoring times.}
#' \item{censor_status}{Binary indicator (1 for failure, 0 for censored at T_star).}
#' \item{n_failures}{Number of observed failures.}
#' \item{termination_time}{Effective termination time of the test.}
#' \item{scheme}{Character string indicating censoring scheme.}
#' @export
#'
#' @examples
#' gen_type1_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, r = 10,
#'   T_star = 1.5, seed = 123
#' )
gen_type1_hybrid <- function(pdf, cdf, lower = 0, upper = Inf, n, r, T_star, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  u <- stats::runif(n)
  inv_cdf <- function(p) {
    sapply(p, function(pv) {
      if (pv <= 0) return(lower)
      if (pv >= 1) return(upper)
      res <- tryCatch(
        stats::uniroot(function(x) cdf(x) - pv, lower = lower, upper = ifelse(is.infinite(upper), 1e5, upper))$root,
        error = function(e) NA
      )
      if (is.na(res)) lower else res
    })
  }
  x <- sort(inv_cdf(u))
  # Type-I Hybrid: test stops at T_stop = min(X_(r), T_star)
  x_r <- if (r <= n) x[r] else Inf
  T_stop <- min(x_r, T_star)
  
  obs_times <- pmin(x, T_stop)
  censor_status <- ifelse(x <= T_stop, 1, 0)
  obs_times <- obs_times[censor_status == 1 | obs_times <= T_stop]
  
  res <- list(
    observed_times = obs_times,
    censor_status = censor_status,
    n_failures = sum(censor_status),
    termination_time = T_stop,
    scheme = "Type-I Hybrid Censoring"
  )
  class(res) <- "comp_risk_rel_data"
  return(res)
}

#' Generate Data under Type-II Hybrid Censoring Scheme
#'
#' @param pdf Probability density function.
#' @param cdf Cumulative distribution function.
#' @param lower Lower bound of the support of the distribution.
#' @param upper Upper bound of the support of the distribution.
#' @param n Total sample size.
#' @param r Minimum required number of failures.
#' @param T_star Target observation time limit.
#' @param seed Optional integer random seed for reproducibility.
#'
#' @return A list containing observed failure times, censor statuses, and termination details.
#' @export
#'
#' @examples
#' gen_type2_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, r = 10,
#'   T_star = 1.0, seed = 123
#' )
gen_type2_hybrid <- function(pdf, cdf, lower = 0, upper = Inf, n, r, T_star, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  u <- stats::runif(n)
  inv_cdf <- function(p) {
    sapply(p, function(pv) {
      if (pv <= 0) return(lower)
      if (pv >= 1) return(upper)
      res <- tryCatch(
        stats::uniroot(function(x) cdf(x) - pv, lower = lower, upper = ifelse(is.infinite(upper), 1e5, upper))$root,
        error = function(e) NA
      )
      if (is.na(res)) lower else res
    })
  }
  x <- sort(inv_cdf(u))
  # Type-II Hybrid: test stops at T_stop = max(X_(r), T_star)
  x_r <- if (r <= n) x[r] else x[n]
  T_stop <- max(x_r, T_star)
  
  obs_times <- x[x <= T_stop]
  censor_status <- rep(1, length(obs_times))
  
  res <- list(
    observed_times = obs_times,
    censor_status = censor_status,
    n_failures = length(obs_times),
    termination_time = T_stop,
    scheme = "Type-II Hybrid Censoring"
  )
  class(res) <- "comp_risk_rel_data"
  return(res)
}

#' Generate Data under Generalized Progressive Hybrid Censoring Scheme
#'
#' @param pdf Probability density function.
#' @param cdf Cumulative distribution function.
#' @param lower Lower bound of support.
#' @param upper Upper bound of support.
#' @param n Total sample size.
#' @param m Target number of failures.
#' @param k Threshold failure count.
#' @param T_star Termination time limit.
#' @param R_plan Progressive censoring plan vector of length m.
#' @param seed Optional integer random seed for reproducibility.
#'
#' @return A list containing observed failure times, progressive removal plan, and censoring details.
#' @export
#'
#' @examples
#' gen_gen_prog_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, m = 10,
#'   k = 5, T_star = 1.2, R_plan = rep(1, 10),
#'   seed = 123
#' )
gen_gen_prog_hybrid <- function(pdf, cdf, lower = 0, upper = Inf, n, m, k, T_star, R_plan, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  if (length(R_plan) != m) R_plan <- rep(floor((n - m) / m), m)
  
  u <- stats::runif(n)
  inv_cdf <- function(p) {
    sapply(p, function(pv) {
      if (pv <= 0) return(lower)
      if (pv >= 1) return(upper)
      res <- tryCatch(
        stats::uniroot(function(x) cdf(x) - pv, lower = lower, upper = ifelse(is.infinite(upper), 1e5, upper))$root,
        error = function(e) NA
      )
      if (is.na(res)) lower else res
    })
  }
  x <- sort(inv_cdf(u))
  
  # Progressive Type-II sampling simulation
  obs_times <- numeric(m)
  rem_pool <- x
  for (i in 1:m) {
    if (length(rem_pool) == 0) break
    obs_times[i] <- rem_pool[1]
    rem_pool <- rem_pool[-1]
    r_i <- R_plan[i]
    if (r_i > 0 && length(rem_pool) >= r_i) {
      drop_idx <- sample(seq_along(rem_pool), r_i)
      rem_pool <- rem_pool[-drop_idx]
    }
  }
  
  # Apply time cutoff T_star according to generalized progressive hybrid scheme
  D <- sum(obs_times <= T_star)
  if (D >= m) {
    final_obs <- obs_times[1:m]
    censor_status <- rep(1, m)
    T_stop <- obs_times[m]
  } else if (D >= k && D < m) {
    final_obs <- obs_times[1:D]
    censor_status <- rep(1, D)
    T_stop <- T_star
  } else {
    eff_m <- min(k, length(obs_times))
    final_obs <- obs_times[1:eff_m]
    censor_status <- rep(1, eff_m)
    T_stop <- obs_times[eff_m]
  }
  
  res <- list(
    observed_times = final_obs,
    censor_status = censor_status,
    n_failures = length(final_obs),
    termination_time = T_stop,
    scheme = "Generalized Progressive Hybrid Censoring"
  )
  class(res) <- "comp_risk_rel_data"
  return(res)
}

#' Generate Data for Competing Risks Analysis
#'
#' @param pdf1 Density function for Cause 1.
#' @param cdf1 Cumulative distribution function for Cause 1.
#' @param pdf2 Density function for Cause 2.
#' @param cdf2 Cumulative distribution function for Cause 2.
#' @param lower Lower bound of support.
#' @param upper Upper bound of support.
#' @param n Total sample size.
#' @param censoring_type Type of censoring ("type1_hybrid", "type2_hybrid", "prog_hybrid").
#' @param r Number of failures required.
#' @param T_star Time limit.
#' @param R_plan Progressive plan (if applicable).
#' @param seed Optional random seed.
#'
#' @return A list containing observed failure times, causes of failure (1 or 2), and censoring flags.
#' @export
#'
#' @examples
#' gen_competing_risks(
#'   pdf1 = function(x) dexp(x, rate = 1),
#'   cdf1 = function(x) pexp(x, rate = 1),
#'   pdf2 = function(x) dexp(x, rate = 1.5),
#'   cdf2 = function(x) pexp(x, rate = 1.5),
#'   lower = 0, upper = 10, n = 25,
#'   censoring_type = "type1_hybrid",
#'   r = 15, T_star = 1.0, seed = 123
#' )
gen_competing_risks <- function(pdf1, cdf1, pdf2, cdf2, lower = 0, upper = Inf, n,
                                censoring_type = c("type1_hybrid", "type2_hybrid", "prog_hybrid"),
                                r = NULL, T_star = NULL, R_plan = NULL, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  censoring_type <- match.arg(censoring_type)
  
  # Generate latent failure times X1 and X2 under independent competing risks
  inv_cdf1 <- function(p) sapply(p, function(pv) stats::uniroot(function(x) cdf1(x) - pv, lower = lower, upper = ifelse(is.infinite(upper), 1e5, upper))$root)
  inv_cdf2 <- function(p) sapply(p, function(pv) stats::uniroot(function(x) cdf2(x) - pv, lower = lower, upper = ifelse(is.infinite(upper), 1e5, upper))$root)
  
  x1 <- inv_cdf1(stats::runif(n))
  x2 <- inv_cdf2(stats::runif(n))
  
  x_min <- pmin(x1, x2)
  cause <- ifelse(x1 <= x2, 1, 2)
  
  ord <- order(x_min)
  x_min <- x_min[ord]
  cause <- cause[ord]
  
  if (censoring_type == "type1_hybrid") {
    r_val <- if (is.null(r)) floor(n / 2) else r
    t_val <- if (is.null(T_star)) stats::median(x_min) else T_star
    T_stop <- min(x_min[r_val], t_val)
    idx <- which(x_min <= T_stop)
    obs_t <- x_min[idx]
    obs_cause <- cause[idx]
  } else if (censoring_type == "type2_hybrid") {
    r_val <- if (is.null(r)) floor(n / 2) else r
    t_val <- if (is.null(T_star)) stats::median(x_min) else T_star
    T_stop <- max(x_min[r_val], t_val)
    idx <- which(x_min <= T_stop)
    obs_t <- x_min[idx]
    obs_cause <- cause[idx]
  } else {
    m_val <- if (is.null(r)) floor(n / 2) else r
    obs_t <- x_min[1:m_val]
    obs_cause <- cause[1:m_val]
    T_stop <- x_min[m_val]
  }
  
  res <- list(
    observed_times = obs_t,
    causes = obs_cause,
    n_failures = length(obs_t),
    termination_time = T_stop,
    scheme = paste("Competing Risks -", censoring_type)
  )
  class(res) <- "comp_risk_rel_data"
  return(res)
}

#' Generate Data for Stress-Strength Reliability Models
#'
#' @param pdf_X Density function for stress X.
#' @param cdf_X Cumulative distribution function for stress X.
#' @param pdf_Y Density function for strength Y.
#' @param cdf_Y Cumulative distribution function for strength Y.
#' @param lower_X Lower bound for X.
#' @param upper_X Upper bound for X.
#' @param lower_Y Lower bound for Y.
#' @param upper_Y Upper bound for Y.
#' @param n_X Sample size for stress X.
#' @param n_Y Sample size for strength Y.
#' @param censoring_type Censoring scheme.
#' @param r_X Minimum required failures for X.
#' @param T_X Time cutoff for X.
#' @param r_Y Minimum required failures for Y.
#' @param T_Y Time cutoff for Y.
#' @param seed Optional random seed.
#'
#' @return A list containing generated stress dataset X and strength dataset Y.
#' @export
#'
#' @examples
#' gen_stress_strength(
#'   pdf_X = function(x) dexp(x, rate = 1.2),
#'   cdf_X = function(x) pexp(x, rate = 1.2),
#'   pdf_Y = function(y) dexp(y, rate = 0.8),
#'   cdf_Y = function(y) pexp(y, rate = 0.8),
#'   n_X = 20, n_Y = 20,
#'   censoring_type = "type1_hybrid",
#'   r_X = 12, T_X = 1.2, r_Y = 12, T_Y = 1.5, seed = 123
#' )
gen_stress_strength <- function(pdf_X, cdf_X, pdf_Y, cdf_Y,
                                lower_X = 0, upper_X = Inf, lower_Y = 0, upper_Y = Inf,
                                n_X, n_Y, censoring_type = c("type1_hybrid", "type2_hybrid", "prog_hybrid"),
                                r_X = NULL, T_X = NULL, r_Y = NULL, T_Y = NULL, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  censoring_type <- match.arg(censoring_type)
  
  rx <- if (is.null(r_X)) floor(n_X * 0.7) else r_X
  tx <- if (is.null(T_X)) 1.5 else T_X
  ry <- if (is.null(r_Y)) floor(n_Y * 0.7) else r_Y
  ty <- if (is.null(T_Y)) 1.5 else T_Y
  
  data_X <- if (censoring_type == "type1_hybrid") {
    gen_type1_hybrid(pdf_X, cdf_X, lower_X, upper_X, n_X, rx, tx)
  } else {
    gen_type2_hybrid(pdf_X, cdf_X, lower_X, upper_X, n_X, rx, tx)
  }
  
  data_Y <- if (censoring_type == "type1_hybrid") {
    gen_type1_hybrid(pdf_Y, cdf_Y, lower_Y, upper_Y, n_Y, ry, ty)
  } else {
    gen_type2_hybrid(pdf_Y, cdf_Y, lower_Y, upper_Y, n_Y, ry, ty)
  }
  
  res <- list(
    data_X = data_X,
    data_Y = data_Y,
    scheme = paste("Stress-Strength Models -", censoring_type)
  )
  class(res) <- "comp_risk_rel_ss_data"
  return(res)
}