tests/testthat/test-frailty-laplace.R

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

test_that("Laplace transform satisfies L(0)=1 and L'(0)=-1 (E[W]=1) for all frailties", {
  frailties <- list(
    list(name = "none", par = numeric(0), exp_var = 0.0),
    list(name = "gamma", par = c(0.8), exp_var = 0.8),
    list(name = "ig", par = c(0.7), exp_var = 0.7),
    list(name = "gl1", par = c(1.2, 0.5), exp_var = (1.2^2 + 0.5^2) / (1.2 + 0.5)),
    list(name = "gl2", par = c(1.5, 0.8), exp_var = ((0.8 - 1.5)^2 + 1) / 1.5)
  )

  for (item in frailties) {
    res0 <- frailty_laplace(s = 0, frailty = item$name, par = item$par)
    expect_equal(res0$L, 1.0, tolerance = 1e-8, label = paste("L(0) for", item$name))
    expect_equal(res0$L1, -1.0, tolerance = 1e-8, label = paste("L'(0) for", item$name))
    expect_equal(res0$Var, item$exp_var, tolerance = 1e-6, label = paste("Var formula for", item$name))
  }
})

test_that("Frailty samplers produce non-negative draws with E[W] approx 1", {
  set.seed(123)
  n_sim <- 20000

  w_ig <- r_ig(n_sim, eta = 0.5)
  expect_true(all(w_ig > 0))
  expect_equal(mean(w_ig), 1.0, tolerance = 0.05)

  w_gl1 <- r_gl1(n_sim, eta = 1.2, epsilon = 0.5)
  expect_true(all(w_gl1 > 0))
  expect_equal(mean(w_gl1), 1.0, tolerance = 0.05)

  w_gl2 <- r_gl2(n_sim, theta = 1.5, mu = 0.8)
  expect_true(all(w_gl2 > 0))
  expect_equal(mean(w_gl2), 1.0, tolerance = 0.05)
})