R/loglik.R

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

#' Log-Likelihood Function for MultiFrailty Models
#'
#' Computes the log-likelihood for shared frailty models across all 10 baseline-frailty combinations
#' with support for right, exact, left, and interval censoring, plus optional progressive censoring.
#'
#' @param par_all Vector of all model parameters on estimation scale (log/logit transformed).
#' @param time Primary event/censoring time vector.
#' @param status Event status vector (0 = right-censored, 1 = exact event, 2 = left-censored, 3 = interval-censored).
#' @param x Matrix of covariates (n x p). Default is 0-column matrix.
#' @param baseline Baseline hazard distribution (\code{"weibull"} or \code{"gw"}).
#' @param frailty Frailty distribution (\code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}).
#' @param time2 Vector of upper interval bounds when \code{status == 3}. Default is NULL.
#' @param prog_cen Vector of progressive censoring counts R_i per observation. Default is NULL.
#'
#' @return Scalar log-likelihood value. Returns \code{-1e12} sentinel on numerical invalidity.
#'
#' @export
#' @examples
#' par_all <- c(log(2), log(1.5), log(0.8), 0.1) # Weibull + Gamma + 1 beta
#' time <- c(1, 2, 3, 4)
#' status <- c(1, 0, 1, 0)
#' x <- matrix(c(0.5, -0.2, 0.1, 0.8), ncol = 1)
#' ll <- loglik_frailty(par_all, time, status, x, baseline = "weibull", frailty = "gamma")
loglik_frailty <- function(par_all, time, status, x = matrix(nrow = length(time), ncol = 0),
                           baseline = c("weibull", "gw"), frailty = c("none", "gamma", "ig", "gl1", "gl2"),
                           time2 = NULL, prog_cen = NULL) {
  baseline <- match.arg(baseline)
  frailty <- match.arg(frailty)

  sentinel <- -1e12

  tryCatch({
    n_base <- if (baseline == "weibull") 2L else 3L

    if (frailty == "none") {
      n_frailty <- 0L
    } else if (frailty %in% c("gamma", "ig")) {
      n_frailty <- 1L
    } else if (frailty %in% c("gl1", "gl2")) {
      n_frailty <- 2L
    }

    n_cov <- ncol(x)
    n_expected <- n_base + n_frailty + n_cov

    if (length(par_all) != n_expected) {
      return(sentinel)
    }

    # Extract baseline parameters on natural scale
    bpar_raw <- par_all[1:n_base]
    bpar <- exp(bpar_raw)

    # Extract frailty parameters on natural scale
    if (frailty == "none") {
      fpar <- numeric(0)
    } else if (frailty == "gamma") {
      fpar <- exp(par_all[n_base + 1])
    } else if (frailty == "ig") {
      fpar <- exp(par_all[n_base + 1])
    } else if (frailty == "gl1") {
      fpar <- exp(par_all[(n_base + 1):(n_base + 2)])
    } else if (frailty == "gl2") {
      theta <- exp(par_all[n_base + 1])
      mu_raw <- par_all[n_base + 2]
      p_val <- pmin(pmax(stats::plogis(mu_raw), 1e-15), 1.0 - 1e-15)
      mu <- (1.0 + theta) * p_val
      fpar <- c(theta, mu)
    }

    # Extract regression coefficients
    if (n_cov > 0) {
      beta <- par_all[(n_base + n_frailty + 1):n_expected]
      rho <- as.vector(exp(x %*% beta))
    } else {
      beta <- numeric(0)
      rho <- rep(1.0, length(time))
    }

    # Evaluate survival and density functions
    ff <- frailty_functions(t = time, eta = rho, frailty = frailty, fpar = fpar,
                            baseline = baseline, bpar = bpar)
    S <- ff$S
    f <- ff$f

    if (any(!is.finite(S)) || any(!is.finite(f)) || any(S <= 0) || any(S >= 1) || any(f <= 0)) {
      return(sentinel)
    }

    log_lik_vec <- numeric(length(time))

    # Exact events (status == 1)
    idx1 <- (status == 1)
    if (any(idx1)) {
      log_lik_vec[idx1] <- log(f[idx1])
    }

    # Right-censored (status == 0)
    idx0 <- (status == 0)
    if (any(idx0)) {
      log_lik_vec[idx0] <- log(S[idx0])
    }

    # Left-censored (status == 2)
    idx2 <- (status == 2)
    if (any(idx2)) {
      S2 <- S[idx2]
      log_lik_vec[idx2] <- log1mexp(-log(S2))
    }

    # Interval-censored (status == 3)
    idx3 <- (status == 3)
    if (any(idx3)) {
      if (is.null(time2)) return(sentinel)
      ff2 <- frailty_functions(t = time2[idx3], eta = rho[idx3], frailty = frailty, fpar = fpar,
                               baseline = baseline, bpar = bpar)
      S2_val <- ff2$S
      diff_S <- S[idx3] - S2_val
      if (any(diff_S <= 0) || any(!is.finite(diff_S))) return(sentinel)
      log_lik_vec[idx3] <- log(diff_S)
    }

    # Progressive censoring extra terms
    if (!is.null(prog_cen) && length(prog_cen) == length(time)) {
      prog_idx <- (prog_cen > 0)
      if (any(prog_idx)) {
        log_lik_vec[prog_idx] <- log_lik_vec[prog_idx] + prog_cen[prog_idx] * log(S[prog_idx])
      }
    }

    val <- sum(log_lik_vec)
    if (!is.finite(val)) return(sentinel)
    return(val)
  }, error = function(e) {
    return(sentinel)
  })
}