R/formula.R

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

#' Formula Interface for Shared Frailty Regression Models
#'
#' Fits shared frailty regression models using a formula specifying survival response
#' (e.g., \code{survival::Surv(time, status) ~ x1 + x2}) and data frame.
#'
#' @param formula Model formula of the form \code{survival::Surv(time, status) ~ x1 + x2}.
#' @param data Data frame containing variables in \code{formula}.
#' @param baseline Character string specifying baseline hazard: \code{"weibull"} or \code{"gw"}.
#' @param frailty Character string specifying frailty family: \code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}.
#' @param method Optimization method for \code{maxLik}: \code{"NR"} (Newton-Raphson) or \code{"BFGS"}.
#' @param ... Additional arguments passed to \code{\link{fit_frailty}}.
#'
#' @return Object of class \code{c("multifrailty", "multifrailty_fit")}.
#'
#' @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
#' set.seed(123)
#' df <- data.frame(
#'   time = stats::rexp(50, rate = 0.1),
#'   status = sample(c(0, 1), 50, replace = TRUE),
#'   age = stats::rnorm(50, mean = 50, sd = 10),
#'   sex = sample(c(0, 1), 50, replace = TRUE)
#' )
#' fit <- multifrailty(survival::Surv(time, status) ~ age + sex, data = df,
#'                     baseline = "weibull", frailty = "gamma")
#' print(fit)
multifrailty <- function(formula, data, baseline = c("weibull", "gw"),
                         frailty = c("none", "gamma", "ig", "gl1", "gl2"),
                         method = "NR", ...) {
  baseline <- match.arg(baseline)
  frailty <- match.arg(frailty)

  if (missing(data) || is.null(data)) {
    data <- environment(formula)
  }

  mf <- stats::model.frame(formula = formula, data = data)
  mt <- attr(mf, "terms")
  y <- stats::model.response(mf, "any")

  if (!inherits(y, "Surv")) {
    stop("Response in formula must be a 'Surv' object (from package 'survival').")
  }

  time <- y[, 1]
  status <- y[, 2]

  x_mat <- stats::model.matrix(mt, mf)
  if ("(Intercept)" %in% colnames(x_mat)) {
    x_mat <- x_mat[, colnames(x_mat) != "(Intercept)", drop = FALSE]
  }

  fit_res <- fit_frailty(time = time, status = status, x = x_mat,
                         baseline = baseline, frailty = frailty,
                         method = method, ...)

  fit_res$call <- match.call()
  fit_res$formula <- formula
  fit_res$terms <- mt

  class(fit_res) <- c("multifrailty", class(fit_res))
  fit_res
}