tests/testthat/test-vi-stan.R

In nlmixr2est: Nonlinear Mixed Effects Models in Population PK/PD, Estimation Routines

## External cross-check of est="fbvi" against Stan's own ADVI (rstan::vb()).  Uses a
## simple linear random-intercept model
##   y_ij = theta + eta_i + err,  eta_i ~ N(0, omega^2), err ~ N(0, sigma^2)
## which has the SAME variational objective in both tools, so the population posterior
## means should agree.
##
## The Stan side is precomputed into tests/testthat/baselines/vi-stan-golden.rds
## (regenerated by tools/viGenStanFixture.R, which needs rstan + a C++ toolchain), so this
## test needs neither at run time.  Stan's answer for a fixed model, fixed data and fixed
## seed does not depend on anything in this package, so recompiling and re-fitting it on
## every run bought nothing -- and gating the whole test behind NLMIXR2_ADVI_STAN=true
## meant the cross-check never actually ran anywhere.

test_that("est='fbvi' population posterior agrees with rstan::vb()", {
  skip_on_cran()
  .gp <- test_path("baselines", "vi-stan-golden.rds")
  skip_if_not(file.exists(.gp), "run tools/viGenStanFixture.R to build the Stan reference")
  .g <- readRDS(.gp)

  ## simulate a linear random-intercept dataset -- must match tools/viGenStanFixture.R
  .testSeed(42)
  nsub <- 60L; nobs <- 6L
  thetaTrue <- 5; omegaTrue <- 1; sigmaTrue <- 0.7
  eta <- stats::rnorm(nsub, 0, omegaTrue)
  dat <- do.call(rbind, lapply(seq_len(nsub), function(i) {
    data.frame(ID = i, TIME = seq_len(nobs),
               DV = thetaTrue + eta[i] + stats::rnorm(nobs, 0, sigmaTrue),
               EVID = 0, AMT = 0)
  }))
  ## Fail safe: the reference is only meaningful for the data Stan actually saw.  If the
  ## simulation drifts (RNG change, edited constants), skip rather than assert against a
  ## reference built from different data.
  skip_if_not(isTRUE(all.equal(dat$DV, .g$dv)),
              "simulated data no longer matches the Stan reference; regenerate it")

  ## est="fbvi" (full-Bayes) on the same model
  linmod <- function() {
    ini({ theta <- 4; eta ~ 1; add.sd <- 1 })
    model({ pred <- theta + eta; pred ~ add(add.sd) })
  }
  fA <- suppressMessages(suppressWarnings(
    nlmixr2(linmod, dat, est = "fbvi",
            control = fbviControl(iters = 800L, print = 0L, returnVi = TRUE))))

  ## population posterior means agree (theta, residual sd, between-subject sd)
  expect_equal(unname(fA$theta[1]), .g$theta, tolerance = 0.15)
  expect_equal(unname(fA$theta[2]), .g$sigma, tolerance = 0.2)   # add.sd
  expect_equal(unname(sqrt(fA$popOmega[1])), .g$omega, tolerance = 0.25)
})