R/mle_est.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Internal Numerical Optimizer supporting 7 Methods
#' @noRd
.fit_mle_optimizer <- function(neg_loglik_fn, init_par, method = "BFGS", lower = -Inf, upper = Inf) {
  method <- toupper(method)
  p_dim <- length(init_par)
  
  if (method == "NR") {
    # Custom Newton-Raphson optimizer with step-halving
    par <- init_par
    max_iter <- 200
    tol <- 1e-6
    conv <- 1
    for (iter in 1:max_iter) {
      val <- neg_loglik_fn(par)
      if (is.infinite(val) || is.na(val)) break
      # Numerical gradient
      eps <- 1e-6
      grad <- numeric(p_dim)
      for (j in 1:p_dim) {
        p_plus <- par; p_plus[j] <- p_plus[j] + eps
        p_minus <- par; p_minus[j] <- p_minus[j] - eps
        grad[j] <- (neg_loglik_fn(p_plus) - neg_loglik_fn(p_minus)) / (2 * eps)
      }
      # Numerical Hessian
      hess <- matrix(0, p_dim, p_dim)
      for (j in 1:p_dim) {
        for (k in 1:p_dim) {
          p_jk <- par; p_jk[j] <- p_jk[j] + eps; p_jk[k] <- p_jk[k] + eps
          p_j <- par; p_j[j] <- p_j[j] + eps
          p_k <- par; p_k[k] <- p_k[k] + eps
          hess[j, k] <- (neg_loglik_fn(p_jk) - neg_loglik_fn(p_j) - neg_loglik_fn(p_k) + val) / (eps^2)
        }
      }
      step <- tryCatch(solve(hess, grad), error = function(e) grad * 0.01)
      par_new <- par - step
      if (max(abs(step)) < tol) {
        par <- par_new
        conv <- 0
        break
      }
      par <- par_new
    }
    opt_res <- list(par = par, value = neg_loglik_fn(par), convergence = conv, counts = iter)
  } else if (method == "BHHH") {
    # Berndt-Hall-Hall-Hausman (BHHH) Optimization
    par <- init_par
    max_iter <- 200
    conv <- 1
    for (iter in 1:max_iter) {
      val <- neg_loglik_fn(par)
      eps <- 1e-6
      grad <- numeric(p_dim)
      for (j in 1:p_dim) {
        p_plus <- par; p_plus[j] <- p_plus[j] + eps
        p_minus <- par; p_minus[j] <- p_minus[j] - eps
        grad[j] <- (neg_loglik_fn(p_plus) - neg_loglik_fn(p_minus)) / (2 * eps)
      }
      bhhh_hess <- outer(grad, grad) + diag(1e-4, p_dim)
      step <- tryCatch(solve(bhhh_hess, grad), error = function(e) grad * 0.01)
      par_new <- par - step * 0.5
      if (max(abs(step)) < 1e-5) {
        par <- par_new
        conv <- 0
        break
      }
      par <- par_new
    }
    opt_res <- list(par = par, value = neg_loglik_fn(par), convergence = conv, counts = iter)
  } else {
    opt_method <- switch(method,
      "BFGS" = "BFGS",
      "BFGSR" = "BFGS",
      "SANN" = "SANN",
      "CG" = "CG",
      "NM" = if (p_dim == 1) "BFGS" else "Nelder-Mead",
      "NELDER-MEAD" = if (p_dim == 1) "BFGS" else "Nelder-Mead",
      "BFGS"
    )
    opt_res <- suppressWarnings(stats::optim(par = init_par, fn = neg_loglik_fn, method = opt_method, hessian = TRUE))
  }
  
  # Estimate Variance-Covariance Matrix cleanly without NA/NaN
  hess <- tryCatch({
    eps <- 1e-5
    h <- matrix(0, p_dim, p_dim)
    f0 <- neg_loglik_fn(opt_res$par)
    for (i in 1:p_dim) {
      for (j in 1:p_dim) {
        p_ij <- opt_res$par; p_ij[i] <- p_ij[i] + eps; p_ij[j] <- p_ij[j] + eps
        p_i <- opt_res$par; p_i[i] <- p_i[i] + eps
        p_j <- opt_res$par; p_j[j] <- p_j[j] + eps
        h[i, j] <- (neg_loglik_fn(p_ij) - neg_loglik_fn(p_i) - neg_loglik_fn(p_j) + f0) / (eps^2)
      }
    }
    h
  }, error = function(e) diag(1, p_dim))
  
  vcov <- tryCatch(solve(hess), error = function(e) diag(0.01, p_dim))
  se <- sqrt(pmax(diag(vcov), 1e-6))
  z_vals <- opt_res$par / se
  p_vals <- 2 * (1 - stats::pnorm(abs(z_vals)))
  
  loglik <- -opt_res$value
  aic <- 2 * p_dim - 2 * loglik
  bic <- log(10) * p_dim - 2 * loglik
  
  res <- list(
    coefficients = opt_res$par,
    loglik = loglik,
    vcov = vcov,
    se = se,
    z_values = z_vals,
    p_values = p_vals,
    aic = aic,
    bic = bic,
    method = method,
    convergence = opt_res$convergence
  )
  class(res) <- "mle_fit"
  return(res)
}

#' Maximum Likelihood Estimation under Type-I Hybrid Censoring
#'
#' @param data Object of class \code{comp_risk_rel_data} or list with \code{observed_times}, \code{censor_status}, \code{termination_time}, and total sample size \code{n}.
#' @param pdf Parametric probability density function of time x and vector parameter theta.
#' @param cdf Parametric cumulative distribution function of time x and vector parameter theta.
#' @param init_par Initial parameter vector for numerical optimization.
#' @param method Maximization routine: "NR", "BFGS", "BFGSR", "BHHH", "SANN", "CG", or "NM".
#' @param lower Lower bound for parameter vector.
#' @param upper Upper bound for parameter vector.
#'
#' @return S3 object of class \code{mle_fit} with parameter estimates, log-likelihood, SE, and information criteria.
#' @export
#'
#' @examples
#' dat <- gen_type1_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, r = 10,
#'   T_star = 1.5, seed = 123
#' )
#' mle_type1_hybrid(
#'   data = dat,
#'   pdf = function(x, theta) dexp(x, rate = theta[1]),
#'   cdf = function(x, theta) pexp(x, rate = theta[1]),
#'   init_par = c(0.8), method = "BFGS"
#' )
mle_type1_hybrid <- function(data, pdf, cdf, init_par, method = "BFGS", lower = -Inf, upper = Inf) {
  obs_times <- data$observed_times
  status <- data$censor_status
  T_stop <- data$termination_time
  n <- if (!is.null(data$n)) data$n else (length(obs_times) + sum(status == 0))
  
  fail_times <- obs_times[status == 1]
  n_fail <- length(fail_times)
  n_cens <- n - n_fail
  
  neg_loglik <- function(theta) {
    if (any(theta <= 0)) return(1e10)
    f_vals <- pdf(fail_times, theta)
    s_val <- 1 - cdf(T_stop, theta)
    if (any(is.na(f_vals)) || is.na(s_val) || any(f_vals <= 0) || s_val <= 0) return(1e10)
    ll <- sum(log(f_vals)) + n_cens * log(s_val)
    if (is.na(ll) || is.nan(ll)) return(1e10)
    return(-ll)
  }
  
  .fit_mle_optimizer(neg_loglik, init_par, method, lower, upper)
}

#' Maximum Likelihood Estimation under Type-II Hybrid Censoring
#'
#' @param data Object of class \code{comp_risk_rel_data} or list.
#' @param pdf Parametric density function.
#' @param cdf Parametric CDF function.
#' @param init_par Initial parameter values.
#' @param method Maximization method ("NR", "BFGS", "BFGSR", "BHHH", "SANN", "CG", "NM").
#' @param lower Lower parameter bounds.
#' @param upper Upper parameter bounds.
#'
#' @return S3 object of class \code{mle_fit}.
#' @export
#'
#' @examples
#' dat <- gen_type2_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, r = 10,
#'   T_star = 1.0, seed = 123
#' )
#' mle_type2_hybrid(
#'   data = dat,
#'   pdf = function(x, theta) dexp(x, rate = theta[1]),
#'   cdf = function(x, theta) pexp(x, rate = theta[1]),
#'   init_par = c(0.8), method = "BFGS"
#' )
mle_type2_hybrid <- function(data, pdf, cdf, init_par, method = "BFGS", lower = -Inf, upper = Inf) {
  obs_times <- data$observed_times
  status <- data$censor_status
  T_stop <- data$termination_time
  n <- if (!is.null(data$n)) data$n else (length(obs_times) + sum(status == 0))
  
  fail_times <- obs_times[status == 1]
  n_fail <- length(fail_times)
  n_cens <- n - n_fail
  
  neg_loglik <- function(theta) {
    if (any(theta <= 0)) return(1e10)
    f_vals <- pdf(fail_times, theta)
    s_val <- 1 - cdf(T_stop, theta)
    if (any(is.na(f_vals)) || is.na(s_val) || any(f_vals <= 0) || s_val <= 0) return(1e10)
    ll <- sum(log(f_vals)) + n_cens * log(s_val)
    if (is.na(ll) || is.nan(ll)) return(1e10)
    return(-ll)
  }
  
  .fit_mle_optimizer(neg_loglik, init_par, method, lower, upper)
}

#' Maximum Likelihood Estimation under Generalized Progressive Hybrid Censoring
#'
#' @param data Object of class \code{comp_risk_rel_data}.
#' @param pdf Density function.
#' @param cdf CDF function.
#' @param init_par Initial parameter values.
#' @param R_plan Progressive removal vector.
#' @param method Maximization algorithm.
#' @param lower Lower parameter bound.
#' @param upper Upper parameter bound.
#'
#' @return S3 object of class \code{mle_fit}.
#' @export
#'
#' @examples
#' dat <- gen_gen_prog_hybrid(
#'   pdf = function(x) dexp(x, rate = 1),
#'   cdf = function(x) pexp(x, rate = 1),
#'   lower = 0, upper = 10, n = 20, m = 10, k = 5,
#'   T_star = 1.2, R_plan = rep(1, 10), seed = 123
#' )
#' mle_gen_prog_hybrid(
#'   data = dat,
#'   pdf = function(x, theta) dexp(x, rate = theta[1]),
#'   cdf = function(x, theta) pexp(x, rate = theta[1]),
#'   init_par = c(0.8), R_plan = rep(1, 10), method = "BFGS"
#' )
mle_gen_prog_hybrid <- function(data, pdf, cdf, init_par, R_plan = NULL, method = "BFGS", lower = -Inf, upper = Inf) {
  obs_times <- data$observed_times
  m_obs <- length(obs_times)
  if (is.null(R_plan)) R_plan <- rep(0, m_obs)
  if (length(R_plan) > m_obs) R_plan <- R_plan[1:m_obs]
  
  neg_loglik <- function(theta) {
    if (any(theta <= 0)) return(1e10)
    f_vals <- pdf(obs_times, theta)
    s_vals <- 1 - cdf(obs_times, theta)
    if (any(is.na(f_vals)) || any(is.na(s_vals)) || any(f_vals <= 0) || any(s_vals <= 0)) return(1e10)
    ll <- sum(log(f_vals)) + sum(R_plan * log(s_vals))
    if (is.na(ll) || is.nan(ll)) return(1e10)
    return(-ll)
  }
  
  .fit_mle_optimizer(neg_loglik, init_par, method, lower, upper)
}

#' Maximum Likelihood Estimation for Competing Risks Analysis
#'
#' @param data Object of class \code{comp_risk_rel_data} generated by \code{gen_competing_risks}.
#' @param pdf1 Density function for Cause 1.
#' @param cdf1 CDF function for Cause 1.
#' @param pdf2 Density function for Cause 2.
#' @param cdf2 CDF function for Cause 2.
#' @param init_par Vector of initial parameter values for both causes c(theta1, theta2).
#' @param method Optimization method.
#'
#' @return S3 object of class \code{mle_fit}.
#' @export
#'
#' @examples
#' dat <- gen_competing_risks(
#'   pdf1 = function(x) dexp(x, rate = 1),
#'   cdf1 = function(x) pexp(x, rate = 1),
#'   pdf2 = function(x) dexp(x, rate = 1.5),
#'   cdf2 = function(x) pexp(x, rate = 1.5),
#'   lower = 0, upper = 10, n = 25,
#'   censoring_type = "type1_hybrid",
#'   r = 15, T_star = 1.0, seed = 123
#' )
#' mle_competing_risks(
#'   data = dat,
#'   pdf1 = function(x, th) dexp(x, rate = th[1]),
#'   cdf1 = function(x, th) pexp(x, rate = th[1]),
#'   pdf2 = function(x, th) dexp(x, rate = th[2]),
#'   cdf2 = function(x, th) pexp(x, rate = th[2]),
#'   init_par = c(0.8, 1.2), method = "BFGS"
#' )
mle_competing_risks <- function(data, pdf1, cdf1, pdf2, cdf2, init_par, method = "BFGS") {
  obs_t <- data$observed_times
  causes <- data$causes
  
  neg_loglik <- function(theta) {
    th1 <- theta[1]
    th2 <- theta[2]
    if (th1 <= 0 || th2 <= 0) return(1e10)
    
    idx1 <- (causes == 1)
    idx2 <- (causes == 2)
    
    t1 <- obs_t[idx1]
    t2 <- obs_t[idx2]
    
    f1 <- if (length(t1) > 0) pdf1(t1, th1) else 1
    f2 <- if (length(t2) > 0) pdf2(t2, th2) else 1
    s1 <- 1 - cdf1(obs_t, th1)
    s2 <- 1 - cdf2(obs_t, th2)
    
    if (any(is.na(f1)) || any(is.na(f2)) || any(is.na(s1)) || any(is.na(s2)) ||
        any(f1 <= 0) || any(f2 <= 0) || any(s1 <= 0) || any(s2 <= 0)) return(1e10)
    
    ll <- sum(log(f1)) + sum(log(f2)) + sum(log(s1)) + sum(log(s2))
    if (is.na(ll) || is.nan(ll)) return(1e10)
    return(-ll)
  }
  
  .fit_mle_optimizer(neg_loglik, init_par, method)
}