tests/testthat/test-Gibbs-sampler.R

In LTFGRS: Implementation of Several Phenotype-Based Family Genetic Risk Scores

################################################################################


skip_if_not_installed("tmvtnorm")

################################################################################

test_that("Gibbs sampler works", {

  cov <- matrix(c(1, 0.2, 0.2, 0.5), 2)
  lower <- c(-Inf, 0)
  upper <- c(0, Inf)

  true <- tmvtnorm::rtmvnorm(100e3, mean = rep(0, 2), H = solve(cov),
                             lower = lower, upper = upper)

  test <- rtmvnorm.gibbs(100e3, covmat = cov, lower = lower, upper = upper, out = 1:2)

  expect_true(all(test[, 1] <= 0))
  expect_true(all(test[, 2] >= 0))
  expect_equal(colMeans(test), colMeans(true), tolerance = 0.01)
  expect_equal(cov(test), cov(true), tolerance = 0.015)
})

################################################################################

test_that("Gibbs sampler works with any cov", {

  cov <- matrix(c(1, 0.2, 0.2, 0.5), 2) * 100
  lower <- c(-Inf, 10)
  upper <- c(0, 20)

  true <- tmvtnorm::rtmvnorm(100e3, mean = rep(0, 2), H = solve(cov),
                             lower = lower, upper = upper, algorithm = "gibbs")

  test <- rtmvnorm.gibbs(100e3, covmat = cov, lower = lower, upper = upper, out = 1:2)

  expect_true(all(test[, 1] <= 0))
  expect_true(all(test[, 2] >= 10))
  expect_true(all(test[, 2] <= 20))
  expect_equal(colMeans(test), colMeans(true), tolerance = 0.01)
  expect_equal(cov(test), cov(true), tolerance = 0.5)
})

################################################################################

test_that("Gibbs sampler with fixed values works", {
  set.seed(555)

  cov <- matrix(c(1, 0.2, 0.2, 0.5), 2)
  lower <- c(-Inf, 0)
  upper <- c(0, 1e-5)

  true <- tmvtnorm::rtmvnorm(100e3, mean = rep(0, 2), H = solve(cov),
                             lower = lower, upper = upper, algorithm = "gibbs")

  test <- rtmvnorm.gibbs(100e3, covmat = cov, out = 1:2,
                         lower = lower, upper = upper, fixed = c(FALSE, TRUE))

  expect_true(all(test[, 1] <= 0))
  val <- test[1, 2]
  expect_gt(val, 0)
  expect_lt(val, 1e-5)
  expect_true(all(test[, 2] == val))
  expect_equal(colMeans(test), colMeans(true), tolerance = 0.01)
  expect_equal(cov(test), cov(true), tolerance = 0.075)
})

################################################################################