R/simulate.R

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

#' Monte Carlo Simulation Framework for MultiFrailty Models
#'
#' Evaluates frequentist Maximum Likelihood Estimation performance (bias, relative bias, MSE, empirical coverage)
#' across repeated Monte Carlo replications under user-specified baseline, frailty, and censoring schemes.
#'
#' @param n_sim Number of Monte Carlo simulation replicates. Default is 50.
#' @param n Sample size per replicate. Default is 100.
#' @param baseline Character string for baseline hazard (\code{"weibull"} or \code{"gw"}).
#' @param bpar True baseline parameter vector.
#' @param frailty Character string for frailty family (\code{"none"}, \code{"gamma"}, \code{"ig"}, \code{"gl1"}, or \code{"gl2"}).
#' @param fpar True frailty parameter vector.
#' @param beta True regression parameter. Default 0.5.
#' @param cen_type Censoring mechanism. Default \code{"right"}.
#' @param cen_rate Censoring rate. Default 0.2.
#'
#' @return A data frame summarizing parameter true values, mean estimates, bias, relative bias, MSE, and coverage.
#'
#' @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
#' \donttest{
#' set.seed(123)
#' sim_res <- simulate_mle_performance(n_sim = 10, n = 50, baseline = "weibull",
#'                                     bpar = c(2, 1.5), frailty = "gamma", fpar = c(0.8))
#' print(sim_res)
#' }
simulate_mle_performance <- function(n_sim = 50, n = 100, baseline = "weibull",
                                     bpar = c(2, 1.5), frailty = "gamma", fpar = c(0.8),
                                     beta = 0.5, cen_type = "right", cen_rate = 0.2) {
  n_base <- if (baseline == "weibull") 2L else 3L
  n_frail <- if (frailty == "none") 0L else if (frailty %in% c("gamma", "ig")) 1L else 2L
  n_cov <- if (length(beta) > 0) length(beta) else 0L

  true_pars <- c(bpar, fpar, beta)
  n_tot <- length(true_pars)

  est_matrix <- matrix(NA_real_, nrow = n_sim, ncol = n_tot)
  se_matrix <- matrix(NA_real_, nrow = n_sim, ncol = n_tot)
  conv_vec <- logical(n_sim)

  for (s in 1:n_sim) {
    x_mat <- if (n_cov > 0) matrix(stats::rnorm(n * n_cov), ncol = n_cov) else matrix(nrow = n, ncol = 0)
    dat <- r_frailty(n = n, baseline = baseline, bpar = bpar, frailty = frailty, fpar = fpar,
                     x = x_mat, beta = beta, cen_type = cen_type, cen_rate = cen_rate)

    fit <- tryCatch({
      fit_frailty(time = dat$time, status = dat$status, x = x_mat, baseline = baseline, frailty = frailty)
    }, error = function(e) NULL)

    if (!is.null(fit) && fit$converged) {
      est_matrix[s, ] <- fit$coefficients$Estimate
      se_matrix[s, ] <- fit$coefficients$StdErr
      conv_vec[s] <- TRUE
    }
  }

  valid_idx <- which(conv_vec)
  n_valid <- length(valid_idx)

  if (n_valid == 0) {
    stop("All Monte Carlo simulation runs failed to converge.")
  }

  mean_est <- colMeans(est_matrix[valid_idx, , drop = FALSE], na.rm = TRUE)
  bias_val <- mean_est - true_pars
  rel_bias <- bias_val / true_pars
  mse_val <- colMeans((est_matrix[valid_idx, , drop = FALSE] - matrix(true_pars, nrow = n_valid, ncol = n_tot, byrow = TRUE))^2, na.rm = TRUE)

  cov_count <- numeric(n_tot)
  for (j in 1:n_tot) {
    lowers <- est_matrix[valid_idx, j] - 1.96 * se_matrix[valid_idx, j]
    uppers <- est_matrix[valid_idx, j] + 1.96 * se_matrix[valid_idx, j]
    cov_count[j] <- mean(true_pars[j] >= lowers & true_pars[j] <= uppers, na.rm = TRUE)
  }

  res_df <- data.frame(
    Parameter = rownames(fit$coefficients),
    True_Value = true_pars,
    Mean_Estimate = mean_est,
    Bias = bias_val,
    Rel_Bias = rel_bias,
    MSE = mse_val,
    Coverage = cov_count,
    Convergence_Rate = n_valid / n_sim,
    stringsAsFactors = FALSE
  )
  rownames(res_df) <- NULL
  res_df
}

#' Run Comprehensive Simulation Study
#'
#' @param sample_sizes Vector of sample sizes. Default c(25, 50, 100).
#' @param n_sim Number of simulation replicates per setting. Default 20.
#' @param baseline Baseline hazard. Default \code{"weibull"}.
#' @param frailty Frailty distribution. Default \code{"gl2"}.
#' @return Master data frame of Monte Carlo simulation metrics across sample sizes.
#' @export
simulation_study <- function(sample_sizes = c(25, 50, 100), n_sim = 20,
                             baseline = "weibull", frailty = "gl2") {
  master_list <- list()

  bpar_def <- if (baseline == "weibull") c(48.0, 1.22) else c(0.02, 1.2, 1.1)
  fpar_def <- if (frailty == "gl2") c(5.0, 2.4) else if (frailty == "gl1") c(1.2, 0.5) else c(0.8)

  for (n_val in sample_sizes) {
    df_n <- simulate_mle_performance(n_sim = n_sim, n = n_val, baseline = baseline,
                                    bpar = bpar_def, frailty = frailty, fpar = fpar_def,
                                    beta = -0.009, cen_type = "right", cen_rate = 0.05)
    df_n$n <- n_val
    master_list[[as.character(n_val)]] <- df_n
  }

  do.call(rbind, master_list)
}