R/baseline.R

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

#' Baseline Hazard and Cumulative Hazard Evaluation
#'
#' Computes baseline hazard \code{phi0(t)} and cumulative baseline hazard \code{Phi0(t)}
#' for Weibull and Generalized Weibull (GW) distributions.
#'
#' @param t Positive numeric vector of time points.
#' @param baseline Character string specifying baseline distribution: \code{"weibull"} or \code{"gw"}.
#' @param par Numeric vector of parameters: \code{c(lambda, gamma)} for Weibull,
#'   or \code{c(delta, zeta, xi)} for GW.
#'
#' @return A named list with components:
#'   \item{H0}{Cumulative baseline hazard \code{Phi0(t)}.}
#'   \item{h0}{Baseline hazard \code{phi0(t)}.}
#'
#' @references
#' Mudholkar, G. S., & Srivastava, D. K. (1993). Exponentiated Weibull family for
#' analyzing bathtub failure-rate data. IEEE Transactions on Reliability, 42(2), 299-302.
#'
#' @export
#' @examples
#' bh_weib <- baseline_hazard(1:5, baseline = "weibull", par = c(2, 1.5))
#' bh_gw <- baseline_hazard(1:5, baseline = "gw", par = c(0.5, 1.2, 1.1))
baseline_hazard <- function(t, baseline = c("weibull", "gw"), par) {
  baseline <- match.arg(baseline)

  if (any(is.na(t)) || any(t <= 0)) {
    stop("All time points 't' must be positive and non-NA.")
  }

  if (baseline == "weibull") {
    if (length(par) != 2 || any(is.na(par)) || any(par <= 0)) {
      stop("Weibull baseline requires 'par = c(lambda, gamma)' with positive values.")
    }
    lambda <- par[1]
    gamma <- par[2]

    log_t_lam <- log(t) - log(lambda)
    H0 <- exp(gamma * log_t_lam)
    h0 <- (gamma / lambda) * exp((gamma - 1) * log_t_lam)
    return(list(H0 = H0, h0 = h0))
  } else if (baseline == "gw") {
    if (length(par) != 3 || any(is.na(par)) || any(par <= 0)) {
      stop("Generalized Weibull baseline requires 'par = c(delta, zeta, xi)' with positive values.")
    }
    delta <- par[1]
    zeta <- par[2]
    xi <- par[3]

    # Cap u to 500 to prevent floating point underflow/overflow at extreme t
    u <- pmin(delta * (t^xi), 500)
    log_1mexp_u <- log1mexp(u)
    v <- pmax(-zeta * log_1mexp_u, 1e-15)

    log_Lambda <- log1mexp(v)
    H0 <- -log_Lambda

    log_phi0 <- log(xi) + log(zeta) + log(delta) + (xi - 1) * log(t) - u + (zeta - 1) * log_1mexp_u - log_Lambda
    h0 <- exp(log_phi0)
    return(list(H0 = H0, h0 = h0))
  }
}

#' Random Generation for Generalized Weibull (GW) Baseline Distribution
#'
#' Generates random variates from the 3-parameter Generalized Weibull baseline distribution.
#'
#' @param n Number of observations to generate.
#' @param delta Scale parameter (delta > 0).
#' @param zeta Shape parameter (zeta > 0).
#' @param xi Shape parameter (xi > 0).
#'
#' @return Numeric vector of length \code{n} containing random samples.
#'
#' @references
#' Mudholkar, G. S., & Srivastava, D. K. (1993). Exponentiated Weibull family for
#' analyzing bathtub failure-rate data. IEEE Transactions on Reliability, 42(2), 299-302.
#'
#' @export
#' @examples
#' set.seed(123)
#' sim_data <- r_gw(100, delta = 0.5, zeta = 1.2, xi = 1.1)
r_gw <- function(n, delta, zeta, xi) {
  if (length(n) > 1) n <- length(n)
  if (n <= 0 || n != as.integer(n)) stop("'n' must be a positive integer.")
  if (delta <= 0 || zeta <= 0 || xi <= 0) stop("Parameters 'delta', 'zeta', 'xi' must be positive.")

  U <- stats::runif(n)
  v <- -(1 / zeta) * log(U)
  log_term <- log1mexp(v)
  t <- ((-1 / delta) * log_term)^(1 / xi)
  t
}