R/functions.R

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

#' Unified Survival, Density, Hazard, and Cumulative Hazard Combiner
#'
#' Evaluates unconditional survival function \code{S(t)}, probability density function \code{f(t)},
#' hazard function \code{h(t)}, and cumulative hazard function \code{H(t)} for any specified
#' combination of frailty distribution and baseline hazard.
#'
#' @param t Positive numeric vector of time points.
#' @param eta Linear predictor risk scores \code{rho = exp(X \%*\% beta)} (numeric vector of same length as \code{t} or scalar).
#' @param frailty Character string specifying frailty family: \code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}.
#' @param fpar Numeric vector of parameters for frailty distribution.
#' @param baseline Character string specifying baseline hazard: \code{"weibull"} or \code{"gw"}.
#' @param bpar Numeric vector of parameters for baseline hazard.
#'
#' @return A named list containing:
#'   \item{S}{Unconditional survival probabilities \code{S(t)}.}
#'   \item{f}{Unconditional density values \code{f(t)}.}
#'   \item{h}{Unconditional hazard values \code{h(t)}.}
#'   \item{H}{Unconditional cumulative hazard values \code{H(t)}.}
#'
#' @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.
#'
#' Pandey, A., & Tyagi, S. (2021). Comparison of Multiplicative Frailty Models Under Weibull Baseline Distribution. Lobachevskii Journal of Mathematics, 42(13), 3184-3195.
#'
#' @export
#' @examples
#' gf <- frailty_functions(1:5, eta = 1, frailty = "gl1", fpar = c(2.13, 1.53),
#'                         baseline = "gw", bpar = c(0.71, 2.13, 1.54))
frailty_functions <- function(t, eta = 1, frailty = c("none", "gamma", "ig", "gl1", "gl2"),
                              fpar = numeric(0), baseline = c("weibull", "gw"), bpar) {
  frailty <- match.arg(frailty)
  baseline <- match.arg(baseline)

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

  bh <- baseline_hazard(t = t, baseline = baseline, par = bpar)
  Phi0 <- bh$H0
  phi0 <- bh$h0

  s_val <- Phi0 * eta

  fl <- frailty_laplace(s = s_val, frailty = frailty, par = fpar)
  L_s <- fl$L
  L1_s <- fl$L1

  S <- pmax(pmin(L_s, 1.0 - 1e-16), 1e-300)
  f <- pmax(- L1_s * phi0 * eta, 1e-300)
  h <- f / S
  H <- - log(S)

  list(S = S, f = f, h = h, H = H)
}