R/frailty.R

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

#' Laplace Transform and Derivatives for Frailty Distributions
#'
#' Computes Laplace transform L(s), first derivative L'(s),
#' and frailty variance Var(W) for five frailty families:
#' \code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, and \code{"gl2"}.
#'
#' @param s Non-negative numeric vector.
#' @param frailty Character string specifying frailty family:
#'   \code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}.
#' @param par Numeric vector of parameters per family:
#'   \code{numeric(0)} for \code{"none"}, \code{c(theta)} for \code{"gamma"},
#'   \code{c(eta)} for \code{"ig"}, \code{c(eta, epsilon)} for \code{"gl1"},
#'   or \code{c(theta, mu)} for \code{"gl2"}.
#'
#' @return A named list with components:
#'   \item{L}{Laplace transform L(s).}
#'   \item{L1}{First derivative L'(s).}
#'   \item{Var}{Variance of frailty distribution Var(W), evaluated as L''(0) - 1.}
#'
#' @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
#' fl_gamma <- frailty_laplace(s = 0.5, frailty = "gamma", par = c(0.8))
#' fl_gl1   <- frailty_laplace(s = 0.5, frailty = "gl1", par = c(1.2, 0.5))
#' fl_gl2   <- frailty_laplace(s = 0.5, frailty = "gl2", par = c(1.5, 0.8))
frailty_laplace <- function(s, frailty = c("none", "gamma", "ig", "gl1", "gl2"), par = numeric(0)) {
  frailty <- match.arg(frailty)

  if (any(is.na(s))) {
    stop("Argument 's' cannot contain NA values.")
  }

  s <- pmax(s, 0.0)

  if (frailty == "none") {
    L <- exp(-s)
    L1 <- -L
    var_w <- 0
    return(list(L = L, L1 = L1, Var = var_w))
  } else if (frailty == "gamma") {
    if (length(par) != 1 || any(is.na(par)) || par[1] <= 0) {
      stop("Gamma frailty requires 'par = c(theta)' with theta > 0.")
    }
    theta <- par[1]
    log_base <- log1p(theta * s)
    log_L <- - (1 / theta) * log_base
    L <- exp(log_L)
    L1 <- - exp(- ((1 + theta) / theta) * log_base)
    var_w <- theta
    return(list(L = L, L1 = L1, Var = var_w))
  } else if (frailty == "ig") {
    if (length(par) != 1 || any(is.na(par)) || par[1] <= 0) {
      stop("Inverse Gaussian frailty requires 'par = c(eta)' with eta > 0.")
    }
    eta <- par[1]
    inner <- sqrt(1 + 2 * eta * s)
    log_L <- (1 - inner) / eta
    L <- exp(log_L)
    L1 <- - L / inner
    var_w <- eta
    return(list(L = L, L1 = L1, Var = var_w))
  } else if (frailty == "gl1") {
    if (length(par) != 2 || any(is.na(par)) || any(par <= 0)) {
      stop("GL Type 1 frailty requires 'par = c(eta, epsilon)' with positive values.")
    }
    eta <- par[1]
    epsilon <- par[2]

    log_denom <- log(eta + epsilon)

    term1_log <- log(eta) - (1 / eta) * log1p(s * eta) - log_denom
    term2_log <- log(epsilon) - (1 / epsilon) * log1p(s * epsilon) - log_denom

    L <- exp(log_sum_exp(term1_log, term2_log))

    l1_term1_log <- log(eta) - ((1 + eta) / eta) * log1p(s * eta) - log_denom
    l1_term2_log <- log(epsilon) - ((1 + epsilon) / epsilon) * log1p(s * epsilon) - log_denom

    L1 <- - exp(log_sum_exp(l1_term1_log, l1_term2_log))

    var_w <- (eta^2 + epsilon^2) / (eta + epsilon)
    return(list(L = L, L1 = L1, Var = var_w))
  } else if (frailty == "gl2") {
    if (length(par) != 2 || any(is.na(par)) || par[1] <= 0 || par[2] <= 0 || par[2] >= (1 + par[1])) {
      stop("GL Type 2 frailty requires 'par = c(theta, mu)' with theta > 0 and 0 < mu < 1 + theta.")
    }
    theta <- par[1]
    mu <- par[2]
    nu <- theta * (1 + theta - mu)

    log_denom <- log(1 + theta)
    log_stheta <- log(s + theta)

    term1_log <- (mu + 1) * log(theta) - mu * log_stheta - log_denom
    term2_log <- nu * log(theta) - nu * log_stheta - log_denom

    L <- exp(log_sum_exp(term1_log, term2_log))

    l1_term1_log <- log(mu) + (mu + 1) * log(theta) - (mu + 1) * log_stheta - log_denom
    l1_term2_log <- log(nu) + nu * log(theta) - (nu + 1) * log_stheta - log_denom

    L1 <- - exp(log_sum_exp(l1_term1_log, l1_term2_log))

    var_w <- ((mu - theta)^2 + 1) / theta
    return(list(L = L, L1 = L1, Var = var_w))
  }
}

#' Random Generation for Inverse Gaussian (IG) Frailty Distribution
#'
#' Generates random variates from the Inverse Gaussian frailty distribution with mean E[W] = 1.
#'
#' @param n Number of observations to generate.
#' @param eta Frailty variance parameter (eta > 0).
#'
#' @return Numeric vector of length \code{n}.
#' @references
#' Hougaard, P. (1984). Life table methods for heterogeneous populations: distributions of frailties. Biometrika, 71(1), 75-83.
#' @export
#' @examples
#' set.seed(123)
#' w_ig <- r_ig(100, eta = 0.5)
r_ig <- function(n, eta) {
  if (length(n) > 1) n <- length(n)
  if (n <= 0 || n != as.integer(n)) stop("'n' must be a positive integer.")
  if (eta <= 0) stop("Parameter 'eta' must be positive.")

  mu <- 1.0
  lambda <- 1.0 / eta

  Z <- stats::rnorm(n)
  Y <- Z^2
  X <- mu + (Y * mu^2) / (2 * lambda) - (mu / (2 * lambda)) * sqrt(4 * mu * lambda * Y + mu^2 * Y^2)

  U <- stats::runif(n)
  w <- ifelse(U <= mu / (mu + X), X, mu^2 / X)
  w
}

#' Random Generation for Generalized Lindley Type 1 (GL1) Frailty Distribution
#'
#' Generates random variates from GL Type 1 frailty distribution (two-component Gamma mixture).
#'
#' @param n Number of observations to generate.
#' @param eta Parameter eta > 0.
#' @param epsilon Parameter epsilon > 0.
#'
#' @return Numeric vector of length \code{n}.
#' @references
#' 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.
#' @export
#' @examples
#' set.seed(123)
#' w_gl1 <- r_gl1(100, eta = 1.2, epsilon = 0.5)
r_gl1 <- function(n, eta, epsilon) {
  if (length(n) > 1) n <- length(n)
  if (n <= 0 || n != as.integer(n)) stop("'n' must be a positive integer.")
  if (eta <= 0 || epsilon <= 0) stop("Parameters 'eta' and 'epsilon' must be positive.")

  p1 <- eta / (eta + epsilon)
  comp <- stats::runif(n) < p1

  w <- numeric(n)
  n1 <- sum(comp)
  n2 <- n - n1

  if (n1 > 0) w[comp] <- stats::rgamma(n1, shape = 1 / eta, scale = eta)
  if (n2 > 0) w[!comp] <- stats::rgamma(n2, shape = 1 / epsilon, scale = epsilon)

  w
}

#' Random Generation for Generalized Lindley Type 2 (GL2) Frailty Distribution
#'
#' Generates random variates from GL Type 2 frailty distribution (two-component Gamma mixture with common rate).
#'
#' @param n Number of observations to generate.
#' @param theta Parameter theta > 0.
#' @param mu Parameter mu in (0, 1 + theta).
#'
#' @return Numeric vector of length \code{n}.
#' @references
#' 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)
#' w_gl2 <- r_gl2(100, theta = 1.5, mu = 0.8)
r_gl2 <- function(n, theta, mu) {
  if (length(n) > 1) n <- length(n)
  if (n <= 0 || n != as.integer(n)) stop("'n' must be a positive integer.")
  if (theta <= 0 || mu <= 0 || mu >= (1 + theta)) stop("'theta' must be positive and 'mu' must be in (0, 1 + theta).")

  nu <- theta * (1 + theta - mu)
  p1 <- theta / (1 + theta)
  comp <- stats::runif(n) < p1

  w <- numeric(n)
  n1 <- sum(comp)
  n2 <- n - n1

  if (n1 > 0) w[comp] <- stats::rgamma(n1, shape = mu, rate = theta)
  if (n2 > 0) w[!comp] <- stats::rgamma(n2, shape = nu, rate = theta)

  w
}