tests/testthat/test-functions.R

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

test_that("frailty_functions produces monotone decreasing S(t) and consistent h(t) = f(t)/S(t)", {
  t_vec <- seq(0.1, 10, length.out = 100)
  combos <- list(
    list(f = "none", fp = numeric(0), b = "weibull", bp = c(2, 1.5)),
    list(f = "gamma", fp = c(0.8), b = "weibull", bp = c(2, 1.5)),
    list(f = "ig", fp = c(0.7), b = "weibull", bp = c(2, 1.5)),
    list(f = "gl1", fp = c(1.2, 0.5), b = "weibull", bp = c(2, 1.5)),
    list(f = "gl2", fp = c(1.5, 0.8), b = "weibull", bp = c(2, 1.5)),
    list(f = "none", fp = numeric(0), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gamma", fp = c(0.8), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "ig", fp = c(0.7), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gl1", fp = c(1.2, 0.5), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gl2", fp = c(1.5, 0.8), b = "gw", bp = c(0.5, 1.2, 1.1))
  )

  for (item in combos) {
    res <- frailty_functions(t_vec, eta = 1.2, frailty = item$f, fpar = item$fp, baseline = item$b, bpar = item$bp)

    expect_true(all(diff(res$S) <= 1e-12), label = paste("Monotonicity for", item$f, item$b))
    expect_equal(res$h, res$f / res$S, tolerance = 1e-8, label = paste("h = f/S for", item$f, item$b))
    expect_equal(res$H, -log(res$S), tolerance = 1e-8, label = paste("H = -log(S) for", item$f, item$b))
  }
})

test_that("Density f(t) integrates to 1.0 for all 10 combinations", {
  combos <- list(
    list(f = "none", fp = numeric(0), b = "weibull", bp = c(2, 1.5)),
    list(f = "gamma", fp = c(0.8), b = "weibull", bp = c(2, 1.5)),
    list(f = "ig", fp = c(0.7), b = "weibull", bp = c(2, 1.5)),
    list(f = "gl1", fp = c(1.2, 0.5), b = "weibull", bp = c(2, 1.5)),
    list(f = "gl2", fp = c(1.5, 0.8), b = "weibull", bp = c(2, 1.5)),
    list(f = "none", fp = numeric(0), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gamma", fp = c(0.8), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "ig", fp = c(0.7), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gl1", fp = c(1.2, 0.5), b = "gw", bp = c(0.5, 1.2, 1.1)),
    list(f = "gl2", fp = c(1.5, 0.8), b = "gw", bp = c(0.5, 1.2, 1.1))
  )

  for (item in combos) {
    int_val <- stats::integrate(function(t) {
      frailty_functions(t, eta = 1.0, frailty = item$f, fpar = item$fp, baseline = item$b, bpar = item$bp)$f
    }, lower = 1e-6, upper = 10000)$value

    expect_equal(int_val, 1.0, tolerance = 0.05, label = paste("Integration for", item$f, item$b))
  }
})