R/censor_gr_test.R

In LugsailGR: Generalized Gelman-Rubin Diagnostic and Effective Sample Size for MCMC

#' Gelman-Rubin Convergence Diagnostics for Censored Data Models
#'
#' Evaluates lugsail Gelman-Rubin convergence diagnostics, effective sample size, 
#' and parameter estimates for MCMC chains fitted to censored data models under 
#' various censoring schemes (right, left, interval, Type-I, Type-II, progressive, and hybrid).
#'
#' @param pdf_fn Function `function(x, par)` returning probability density values for variable `x` given parameters `par`.
#' @param cdf_fn Function `function(x, par)` returning cumulative distribution function values.
#' @param surv_fn Optional function `function(x, par)` returning survival function values. If `NULL`, computed as `1 - cdf_fn(x, par)`.
#' @param data Numeric vector or matrix of observed sample data.
#' @param start_par Initial parameter vector.
#' @param censor_type Character string specifying censoring scheme. Options include:
#'   \code{"right"} (Right censoring), \code{"left"} (Left censoring), \code{"interval"} (Interval censoring),
#'   \code{"type1"} (Type-I censoring), \code{"type2"} (Type-II censoring),
#'   \code{"progressive"} (Progressive Type-II censoring), or \code{"hybrid"} (Hybrid censoring).
#' @param status Numeric vector or matrix indicating event/censoring status (e.g. 1 = observed, 0 = censored).
#' @param censor_info Optional list specifying parameters for censoring schemes (e.g. cutoff time \code{T_censor}, number of failures \code{r_censor}, or removal pattern \code{R_pattern}).
#' @param n_iter Integer, total number of MCMC iterations per chain (default \code{2000}).
#' @param n_chains Integer, number of parallel chains (default \code{3}).
#' @param burn_in Numeric, burn-in fraction (default \code{0.50}).
#' @param scale Numeric or vector, proposal scale for random walk Metropolis-Hastings (default \code{0.5}).
#' @param alpha Numeric, significance level (default \code{0.05}).
#' @param epsilon Numeric, relative volume tolerance (default \code{0.10}).
#'
#' @return An object of class \code{"lugsail_gr"} containing Gelman-Rubin diagnostics, 
#'   effective sample sizes, cutoff thresholds, parameter estimates, and MCMC samples.
#'
#' @references
#' Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. 
#' \emph{Statistical Science}, 36(4), 518–529. \doi{10.1214/20-STS812}.
#'
#' @export
#'
#' @examples
#' # Example: Right-censored exponential distribution
#' set.seed(123)
#' n_obs <- 30
#' true_rate <- 0.5
#' times <- rexp(n_obs, rate = true_rate)
#' censor_time <- 2.0
#' obs_times <- pmin(times, censor_time)
#' status <- as.numeric(times <= censor_time)
#'
#' f_exp <- function(x, rate) dexp(x, rate = rate)
#' F_exp <- function(x, rate) pexp(x, rate = rate)
#' S_exp <- function(x, rate) 1 - pexp(x, rate = rate)
#'
#' fit <- censor_gr_test(pdf_fn = f_exp, cdf_fn = F_exp, surv_fn = S_exp,
#'                      data = obs_times, start_par = c(rate = 0.8),
#'                      censor_type = "right", status = status,
#'                      n_iter = 1000, n_chains = 3)
#' print(fit)
censor_gr_test <- function(pdf_fn, cdf_fn, surv_fn = NULL, data, start_par,
                           censor_type = c("right", "left", "interval", "type1", "type2", "progressive", "hybrid"),
                           status = NULL, censor_info = NULL, n_iter = 2000,
                           n_chains = 3, burn_in = 0.5, scale = 0.5,
                           alpha = 0.05, epsilon = 0.10) {
  censor_type <- match.arg(censor_type)

  if (is.null(surv_fn)) {
    surv_fn <- function(x, par) {
      pmax(1e-10, 1 - cdf_fn(x, par))
    }
  }

  # Build log-likelihood function for specified censoring scheme
  censored_log_post <- function(par, data) {
    if (any(is.na(par)) || any(is.nan(par))) return(-Inf)

    ll <- 0.0
    if (censor_type == "right") {
      if (is.null(status)) status_vec <- rep(1, length(data)) else status_vec <- status
      f_vals <- suppressWarnings(pdf_fn(data[status_vec == 1], par))
      s_vals <- suppressWarnings(surv_fn(data[status_vec == 0], par))
      f_vals <- pmax(1e-10, f_vals)
      s_vals <- pmax(1e-10, s_vals)
      ll <- sum(log(f_vals)) + sum(log(s_vals))
    } else if (censor_type == "left") {
      if (is.null(status)) status_vec <- rep(1, length(data)) else status_vec <- status
      f_vals <- suppressWarnings(pdf_fn(data[status_vec == 1], par))
      cdf_vals <- suppressWarnings(cdf_fn(data[status_vec == 0], par))
      f_vals <- pmax(1e-10, f_vals)
      cdf_vals <- pmax(1e-10, cdf_vals)
      ll <- sum(log(f_vals)) + sum(log(cdf_vals))
    } else if (censor_type == "interval") {
      if (is.matrix(data) || is.data.frame(data)) {
        L_vals <- data[, 1]
        U_vals <- data[, 2]
        diff_cdf <- suppressWarnings(cdf_fn(U_vals, par) - cdf_fn(L_vals, par))
        diff_cdf <- pmax(1e-10, diff_cdf)
        ll <- sum(log(diff_cdf))
      } else {
        f_vals <- suppressWarnings(pdf_fn(data, par))
        f_vals <- pmax(1e-10, f_vals)
        ll <- sum(log(f_vals))
      }
    } else if (censor_type %in% c("type1", "type2", "progressive", "hybrid")) {
      if (!is.null(status)) {
        f_vals <- suppressWarnings(pdf_fn(data[status == 1], par))
        s_vals <- suppressWarnings(surv_fn(data[status == 0], par))
        f_vals <- pmax(1e-10, f_vals)
        s_vals <- pmax(1e-10, s_vals)
        ll <- sum(log(f_vals)) + sum(log(s_vals))
      } else {
        f_vals <- suppressWarnings(pdf_fn(data, par))
        f_vals <- pmax(1e-10, f_vals)
        ll <- sum(log(f_vals))
      }
    }

    if (is.na(ll) || is.nan(ll)) return(-Inf)
    return(ll)
  }

  res <- mcmc_gr_test(target_pdf = censored_log_post, data = data,
                      start_par = start_par, n_iter = n_iter,
                      n_chains = n_chains, burn_in = burn_in,
                      scale = scale, alpha = alpha, epsilon = epsilon)
  res$censor_type <- censor_type
  return(res)
}