tests/testthat/test-impmap.R

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

nmTest({
  test_that("impmapControl option sanity and IS defaults", {
    expect_error(impmapControl(), NA)
    .ctl <- impmapControl()
    expect_s3_class(.ctl, "impmapControl")
    # inherits the FOCEI/MAP option surface
    expect_false(is.null(.ctl$maxOuterIterations))
    # mu-referencing is forced on for the MAP proposal center
    expect_identical(.ctl$muModel, "lin")
    # importance-sampling / EM defaults
    expect_identical(.ctl$isample, 300L)
    expect_identical(.ctl$nIter, 100L)
    expect_identical(.ctl$mapIter, 1L)
    expect_identical(.ctl$gamma, 1.0)
    expect_identical(.ctl$iscaleMin, 0.1)
    expect_identical(.ctl$iscaleMax, 10.0)
    expect_identical(.ctl$iaccept, 0.4)
    expect_null(.ctl$ctol)
    # nConvWindow is a RULE-DEPENDENT default: the "target" rule (now the default)
    # tracks a Monte-Carlo statistic and needs a longer window to average it out.
    expect_identical(.ctl$nConvWindow, 20L)
    expect_identical(impmapControl(gammaRule = "floor")$nConvWindow, 10L)
    expect_identical(.ctl$impSeed, 42L)

    # round-trips through do.call (used by getValidNlmixrCtl / .foceiFamilyControl)
    expect_error(do.call(impmapControl, .ctl), NA)
  })

  test_that("impmapControl overrides and focei passthrough", {
    .ctl <- impmapControl(isample = 50L, gamma = 2, impSeed = 7L,
                          maxOuterIterations = 3L)
    expect_identical(.ctl$isample, 50L)
    expect_identical(.ctl$gamma, 2.0)
    expect_identical(.ctl$impSeed, 7L)
    # forwarded to foceiControl via ...
    expect_identical(.ctl$maxOuterIterations, 3L)
  })

  test_that("down-conversion to foceiControl strips IS-only names", {
    .env <- new.env()
    .env$impmapControl <- impmapControl()
    .fc <- nlmixr2est:::.impmapControlToFoceiControl(.env, assign = FALSE)
    expect_s3_class(.fc, "foceiControl")
    # IS/EM-only names must not leak into the plain foceiControl
    for (.n in nlmixr2est:::.impmapIsControlNames) {
      expect_null(.fc[[.n]])
    }
    # MAP-relevant focei options are preserved
    expect_identical(.fc$muModel, "lin")
  })

  test_that("impmap dispatch is registered and discoverable", {
    expect_true("impmap" %in% nlmixr2AllEst())
    expect_true(is.function(getS3method("nlmixr2Est", "impmap")))
    # mu-hook activation gate is a control-dependent predicate (mfocei-style)
    expect_true(is.function(attr(nlmixr2Est.impmap, "mu")))
  })

  test_that("getValidNlmixrCtl.impmap yields a default impmapControl", {
    expect_s3_class(getValidNlmixrCtl.impmap(list(NULL)), "impmapControl")
  })

  test_that("M1: impmap MAP pass exposes per-subject mode and Hessian", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .dat <- nlmixr2data::theo_sd
    # One EM iteration is enough to exercise the MAP + per-subject Hessian.
    .imp <- suppressWarnings(
      nlmixr2(one.cmt, .dat, "impmap", impmapControl(print = 0L, nIter = 1L)))

    expect_true(inherits(.imp, "nlmixr2FitCore"))
    # The MAP pass stashes each subject's mode + eta Hessian; check the Hessian is
    # present, square, symmetric, and positive-definite for subject 1.
    .env <- .imp$env
    expect_true(is.matrix(.env$impEtaMode) &&
                  nrow(.env$impEtaMode) == length(unique(.dat$ID)))
    .H <- .env$impEtaHess
    expect_true(is.list(.H) && length(.H) == length(unique(.dat$ID)))
    .H1 <- .H[[1]]
    expect_true(is.matrix(.H1) && all(dim(.H1) == c(2, 2)))
    expect_equal(.H1, t(.H1), tolerance = 1e-6)
    expect_true(all(eigen(.H1, symmetric = TRUE, only.values = TRUE)$values > 0))
  })

  test_that("M2: threefry proposal sampler matches N(mode, gamma*H^-1)", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .dat <- nlmixr2data::theo_sd
    .nsub <- length(unique(.dat$ID))
    .gamma <- 1.5
    rxode2::rxSetSeed(42)
    .f <- suppressWarnings(
      nlmixr2(one.cmt, .dat, "impmap",
              impmapControl(print = 0L, nIter = 1L, isample = 4000L, gamma = .gamma)))
    .e <- .f$env
    expect_identical(.e$impNsample, 4000L)
    expect_equal(.e$impGammaUsed, .gamma)
    .S <- .e$impSamples
    expect_true(is.list(.S) && length(.S) == .nsub)
    expect_true(all(dim(.S[[1]]) == c(4000L, 2L)))
    # subject 1: empirical mean ~ mode, empirical cov ~ gamma * H^-1
    .m1 <- as.numeric(.e$impEtaMode[1, ])
    expect_equal(colMeans(.S[[1]]), .m1, tolerance = 0.03)
    .covTarget <- .gamma * solve(.e$impEtaHess[[1]])
    expect_equal(unname(cov(.S[[1]])), unname(.covTarget), tolerance = 0.02)
  })

  test_that("M2: sampler is thread-count independent (D6)", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .dat <- nlmixr2data::theo_sd
    .thr0 <- rxode2::getRxThreads()
    on.exit(rxode2::setRxThreads(.thr0), add = TRUE)
    # Perturb the ambient RNG differently before each run and do NOT set a seed
    # ourselves: the fit seeds its own E-step from impmapControl(impSeed=), so the
    # samples must be bit-identical regardless of thread count or ambient state.
    .run <- function(nthr) {
      rxode2::setRxThreads(nthr)
      rxode2::rxSetSeed(sample.int(9999L, 1L)); stats::runif(sample.int(50L, 1L))
      suppressWarnings(
        nlmixr2(one.cmt, .dat, "impmap",
                impmapControl(print = 0L, nIter = 1L, isample = 100L)))$env$impSamples
    }
    .s1 <- .run(1L)
    .s4 <- .run(4L)
    for (.i in seq_along(.s1)) {
      expect_identical(.s1[[.i]], .s4[[.i]])
    }
  })

  test_that("M3: importance weights recover the conditional mean/variance", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .dat <- nlmixr2data::theo_sd
    .gamma <- 2
    rxode2::rxSetSeed(42)
    .f <- suppressWarnings(
      nlmixr2(one.cmt, .dat, "impmap",
              impmapControl(print = 0L, nIter = 1L, isample = 6000L, gamma = .gamma)))
    .e <- .f$env
    # E-step outputs present and well-formed
    expect_true(is.numeric(.e$impObj) && is.finite(.e$impObj))
    expect_equal(nrow(.e$impCondMean), length(.e$impNeff))
    expect_true(all(is.finite(.e$impLi)))
    # effective sample size strictly between 1 and isample
    expect_true(all(.e$impNeff > 1 & .e$impNeff <= 6000))

    .k <- 1L
    .mode1 <- as.numeric(.e$impEtaMode[.k, ])
    .H1 <- .e$impEtaHess[[.k]]
    # conditional mean tracks the mode for a near-Gaussian posterior
    expect_equal(as.numeric(.e$impCondMean[.k, ]), .mode1, tolerance = 0.03)
    # KEY: the importance weights reweight samples drawn from the inflated
    # proposal (cov = gamma*H^-1) back to the TRUE posterior covariance ~ H^-1,
    # NOT the proposal covariance.  This is what validates the weights.
    .B1 <- .e$impCondVar[[.k]]
    .postCov <- solve(.H1)
    .propCov <- .gamma * .postCov
    expect_equal(unname(.B1), unname(.postCov), tolerance = 0.02)
    # and B is clearly closer to H^-1 than to the proposal covariance
    expect_lt(max(abs(.B1 - .postCov)), max(abs(.B1 - .propCov)))
  })

  test_that("M4: EM converges to FOCEI on the mu-referenced params and Omega", {
    # Non-mu parameters (tv, add.sd) are held fixed so this isolates the EM update
    # of the mu-referenced thetas and Omega (the non-mu FD updates are a later
    # module); impmap should then match FOCEI on tka/tcl and the Omega diagonal.
    mfix <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- fix(3.45)
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- fix(0.7)
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(mfix, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(mfix, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 40L, isample = 300L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    expect_equal(fixef(.fi)[c("tka", "tcl")], fixef(.ff)[c("tka", "tcl")], tolerance = 0.05)
    expect_equal(unname(diag(.fi$omega)), unname(diag(.ff$omega)), tolerance = 0.1)
  })

  test_that("M4: mu-referenced covariate (updateMuGroups) is estimated", {
    # cl.wt is a mu-referenced covariate effect -- handled by the covariate
    # regression update (updateMuGroups), which impmap must drive.  The estimate
    # should match FOCEI, and the fit should report a nonzero mu covariate group.
    mcov <- function() {
      ini({
        tka <- 0.45; tcl <- 1; cl.wt <- 0.75; tv <- fix(3.45)
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- fix(0.7)
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + cl.wt * log(WT / 70) + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(mcov, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(mcov, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 40L, isample = 300L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    # a mu covariate group was actually set up and driven
    expect_true(.fi$env$impMuGroupN >= 1L)
    # the covariate coefficient matches FOCEI
    expect_equal(unname(fixef(.fi)["cl.wt"]), unname(fixef(.ff)["cl.wt"]), tolerance = 0.03)
    expect_equal(fixef(.fi)[c("tka", "tcl")], fixef(.ff)[c("tka", "tcl")], tolerance = 0.05)
  })

  test_that("M5: non-mu structural theta converges to FOCEI (symbolic sensitivity Newton step)", {
    # tv is a non-mu structural theta -- estimated by the M-step Newton update on
    # the IS-weighted score / Gauss-Newton Hessian built from the symbolic
    # d(f)/d(theta) sensitivity model.  add.sd (sigma) stays fixed (a later module).
    mstr <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- fix(0.7)
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        d/dt(depot) <- -ka * depot
        d/dt(central) <- ka * depot - cl / v * central
        cp <- central / v
        cp ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(mstr, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(mstr, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 30L, isample = 300L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    # the structural theta tv was actually estimated (moved off its start toward FOCEI)
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
    expect_equal(fixef(.fi)[c("tka", "tcl")], fixef(.ff)[c("tka", "tcl")], tolerance = 0.05)
    expect_equal(unname(diag(.fi$omega)), unname(diag(.ff$omega)), tolerance = 0.1)
    # the fixed residual-error theta stays put (not swept into the Newton step)
    expect_equal(unname(fixef(.fi)["add.sd"]), 0.7)
  })

  test_that("M6: residual-error sigma converges to FOCEI (symbolic d(V)/d(sigma))", {
    # add.sd (additive residual error) is estimated by the same M-step Newton update,
    # now including the variance-part score/Fisher from the symbolic d(V)/d(sigma).
    madd <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        d/dt(depot) <- -ka * depot
        d/dt(central) <- ka * depot - cl / v * central
        cp <- central / v
        cp ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(madd, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(madd, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 30L, isample = 300L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.03)
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
  })

  test_that("M6: combined additive+proportional error converges to FOCEI", {
    # proportional error makes V depend on the prediction f, so d(V)/d(theta)
    # couples through the structural theta as well -- the general sensitivity path
    # (d(f)/d(theta) and d(V)/d(theta)) handles both add.sd and prop.sd.
    mcomb <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.5; prop.sd <- 0.1
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        d/dt(depot) <- -ka * depot
        d/dt(central) <- ka * depot - cl / v * central
        cp <- central / v
        cp ~ add(add.sd) + prop(prop.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(mcomb, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(mcomb, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 30L, isample = 300L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
    expect_equal(unname(fixef(.fi)["prop.sd"]), unname(fixef(.ff)["prop.sd"]), tolerance = 0.02)
  })

  test_that("M7: multiple endpoints with more structural thetas than etas (pool sized for theta-sens)", {
    # 2-endpoint PK/PD (indirect response).  Only eta.cl is random, so the inner
    # model has few states while the theta-sensitivity model (tka, tv, tec50, tkout,
    # te0) has many -- exercising the pool-sized-for-the-larger-structure path where
    # the inner MAP runs with ind->neqOverride.  All thetas + both sigmas should
    # converge to FOCEI, and the fit must not crash.
    skip_on_cran()
    mpkpd <- function() {
      ini({
        tka <- 0.5; tcl <- -3.2; tv <- -0.7; tec50 <- 2; tkout <- -2; te0 <- 4.6
        eta.cl ~ 0.09
        add.sd <- 0.4; pdadd.sd <- 2
      })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        ec50 <- exp(tec50); kout <- exp(tkout); e0 <- exp(te0)
        d/dt(depot) <- -ka * depot
        d/dt(center) <- ka * depot - cl / v * center
        cp <- center / v
        effect(0) <- e0
        d/dt(effect) <- kout * (e0 * (1 - cp / (ec50 + cp)) - effect)
        cp ~ add(add.sd) | center
        effect ~ add(pdadd.sd) | effect
      })
    }
    # simulate a 2-endpoint dataset from the model
    .testSeed(1); rxode2::rxSetSeed(1)
    .ev <- rxode2::et(amt = 100, cmt = "depot", id = 1:12)
    .ev <- rxode2::et(.ev, seq(0.5, 24, by = 3), cmt = "center")
    .ev <- rxode2::et(.ev, seq(0.5, 24, by = 3), cmt = "effect")
    .d <- as.data.frame(rxode2::rxSolve(mpkpd, .ev, addDosing = TRUE))
    .dose <- .d[.d$evid != 0, c("id", "time", "CMT", "amt", "evid")]
    .dose$dv <- NA_real_
    names(.dose)[names(.dose) == "CMT"] <- "cmt"
    .obs <- .d[.d$evid == 0, c("id", "time", "CMT", "sim")]
    .obs$amt <- 0; .obs$evid <- 0
    names(.obs)[names(.obs) == "CMT"] <- "cmt"
    names(.obs)[names(.obs) == "sim"] <- "dv"
    .dat <- rbind(.dose[, c("id", "time", "dv", "cmt", "amt", "evid")],
                  .obs[, c("id", "time", "dv", "cmt", "amt", "evid")])
    .dat <- .dat[order(.dat$id, .dat$time, -.dat$evid), ]
    # this ill-conditioned 2-endpoint model (more structural thetas than etas) has a
    # residual sigma sensitive to solve accuracy, so pin sigdig=6 on both fits to
    # compare the methods at a converged tolerance (the sigdig=4 default leaves the
    # FOCEI add.sd ~25% off, which is solve noise, not a method difference).
    rxode2::rxSetSeed(42)
    .ff <- suppressWarnings(nlmixr2(mpkpd, .dat, "focei",
                                    foceiControl(print = 0L, covMethod = "", sigdig = 6)))
    rxode2::rxSetSeed(42)
    .fi <- suppressWarnings(nlmixr2(mpkpd, .dat, "impmap",
                                    impmapControl(print = 0L, nIter = 20L, isample = 300L, sigdig = 6)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    # This fit runs after a prior parallel FOCEI fit.  Two bugs used to degrade it
    # non-deterministically: the within-fit pool-sizing thread race (fixed by forcing
    # this path serial) and, underneath, rxode2 LSODA work memory left uninitialised
    # (malloc) that the extreme-eta E-step solves read (fixed by zeroing it, calloc).
    # With both fixed the fit is deterministic and every subject keeps a healthy
    # effective sample size regardless of thread count.
    .neffFrac <- .fi$env$impNeff / .fi$env$impNsample
    expect_false(anyNA(.neffFrac))
    # The default rule ("target") drives xi onto iaccept, which deliberately
    # widens the proposal and so LOWERS the Kish effective-sample fraction -- that
    # is the trade it makes to bound the weights, not a defect.  The bound below
    # is a sanity floor; pin gammaRule = "floor" if you want the old >0.9.
    expect_true(min(.neffFrac) > 0.4)
    # PD structural thetas (in the higher-state theta-sensitivity model) match FOCEI
    expect_equal(fixef(.fi)[c("tec50", "tkout", "te0")],
                 fixef(.ff)[c("tec50", "tkout", "te0")], tolerance = 0.05)
    # both endpoints' residual-error sigmas match FOCEI (the E-step is seeded from
    # impmapControl(impSeed=) so the fit is reproducible and thread-count independent)
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
    expect_equal(unname(fixef(.fi)["pdadd.sd"]), unname(fixef(.ff)["pdadd.sd"]), tolerance = 0.05)
  })

  test_that("M8: windowed convergence controller stops early and adapts gamma", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.cl ~ 0.1
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(one.cmt, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    .fi <- suppressWarnings(nlmixr2(one.cmt, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 100L, isample = 300L,
                                                  nConvWindow = 10L)))
    .E <- .fi$env
    # the windowed criterion should trip well before the nIter cap on this
    # well-behaved (near-Gaussian) problem
    expect_true(isTRUE(.E$impConverged))
    expect_true(.E$impIter < 100L)
    # per-iteration diagnostics are recorded and internally consistent
    expect_length(.E$impObjTrace, .E$impIter)
    expect_length(.E$impGammaTrace, .E$impIter)
    expect_length(.E$impNeffFrac, .E$impIter)
    # gamma stays within the ISCALE bounds
    expect_true(all(.E$impGammaTrace >= 0.1 - 1e-8 & .E$impGammaTrace <= 10 + 1e-8))
    # "a well-covered proposal is not inflated" is the FLOOR rule's semantics, and
    # this fit runs the default "target" rule, which adjusts gamma both ways until
    # xi approximates iaccept -- so assert THAT instead.  The floor behaviour is
    # covered by the gammaRule test in test-imp-xi-gamma.R.
    expect_equal(unname(tail(.E$impXiTrace, 1)),
                 .E$impmapControl$iaccept, tolerance = 0.1)
    # and the early-stopped fit still matches FOCEI
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
  })

  test_that("C2: MC covariance (covMethod='imp') is the default; theta SEs match FOCEI |r|", {
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6; eta.cl ~ 0.3
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    # covMethod="" turns the covariance step off: no covariance stash
    .f0 <- suppressWarnings(nlmixr2(one.cmt, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 1L, covMethod = "")))
    expect_null(.f0$env$impSe)
    # default covMethod="imp": the full (theta + Omega) MC covariance
    .fi <- suppressWarnings(nlmixr2(one.cmt, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 40L, isample = 500L)))
    .se <- as.numeric(.fi$env$impSe)
    .nth <- .fi$env$impCovThetaN
    expect_true(all(is.finite(.se) & .se > 0))
    # the full covariance is symmetric positive-definite
    .cov <- .fi$env$impCov
    expect_equal(.cov, t(.cov), tolerance = 1e-8)
    expect_true(all(eigen(.cov, symmetric = TRUE, only.values = TRUE)$values > 0))
    # published as the fit covariance so standard errors show in the parameter table
    expect_false(is.null(.fi$cov))
    expect_true(all(is.finite(.fi$parFixedDf[["SE"]][seq_len(.nth)])))
    # theta and Omega rows/columns of vcov() are both labelled
    expect_true(all(c("tka", "tcl", "tv", "add.sd", "om.eta.ka", "om.eta.cl") %in%
                      dimnames(.fi$cov)[[1]]))
    # theta SEs match the Hessian-based FOCEI covariance (|r|)
    .ff <- suppressWarnings(nlmixr2(one.cmt, .d, "focei",
                                    foceiControl(print = 0L, covMethod = "r")))
    skip_if(is.null(.ff$cov), "FOCEI |r| covariance unavailable")
    .fse <- sqrt(diag(.ff$cov))[seq_len(.nth)]
    expect_equal(.se[seq_len(.nth)], unname(.fse), tolerance = 0.1)
  })

  test_that("Censoring: the M-step gradient uses the analytic censored score (matches FOCEI)", {
    m <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.cl ~ 0.1
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    # Left-censor (BLQ / M3): observations below the LOQ carry DV = LOQ, CENS = 1.
    .d <- nlmixr2data::theo_sd
    .loq <- 2.5
    .d$CENS <- 0L
    .blq <- .d$EVID == 0 & .d$DV < .loq & .d$DV > 0
    .d$CENS[.blq] <- 1L
    .d$DV[.blq] <- .loq
    expect_true(sum(.blq) > 5) # the censored branch is actually exercised
    .ff <- suppressWarnings(nlmixr2(m, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    .fi <- suppressWarnings(nlmixr2(m, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 40L, isample = 500L)))
    # the residual sigma is the parameter most sensitive to censoring; the non-mu
    # structural theta goes through the same M-step gradient
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
  })

  test_that("Mixture: recovers the sub-population proportion and the two clearance groups", {
    # two well-separated clearance groups (CL ~ 3 and ~ 9) with true proportion 0.6
    .mkg <- function(cl0, ids) {
      ka <- 1.5; v <- 8
      do.call(rbind, lapply(ids, function(id) {
        cli <- cl0 * exp(stats::rnorm(1, 0, 0.25))
        tt <- seq(0.25, 24, by = 2)
        cp <- (100 * ka / (v * (ka - cli / v))) * (exp(-cli / v * tt) - exp(-ka * tt))
        cp <- pmax(cp, 1e-3) * exp(stats::rnorm(length(tt), 0, 0.12))
        rbind(data.frame(id = id, time = 0, dv = NA_real_, amt = 100, evid = 1, cmt = "depot"),
              data.frame(id = id, time = tt, dv = cp, amt = 0, evid = 0, cmt = "cen"))
      }))
    }
    .testSeed(11); rxode2::rxSetSeed(11)
    .d <- rbind(.mkg(3.0, 1:30), .mkg(9.0, 31:50))
    .d <- .d[order(.d$id, .d$time, -.d$evid), ]
    m <- function() {
      ini({
        tka <- log(1.5); tcl1 <- log(2.5); tcl2 <- log(8); tv <- log(8)
        p1 <- 0.5
        eta.cl ~ 0.2
        add.sd <- 0.3
      })
      model({
        ka <- exp(tka)
        cl <- mix(exp(tcl1 + eta.cl), p1, exp(tcl2 + eta.cl))
        v <- exp(tv)
        d/dt(depot) <- -ka * depot
        d/dt(cen) <- ka * depot - cl / v * cen
        cp <- cen / v
        cp ~ add(add.sd)
      })
    }
    .fi <- suppressWarnings(nlmixr2(m, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 20L, isample = 200L)))
    .cl1 <- exp(unname(fixef(.fi)["tcl1"]))
    .cl2 <- exp(unname(fixef(.fi)["tcl2"]))
    .p1 <- unname(fixef(.fi)["p1"])
    # the mean-posterior proportion update recovers roughly the true 0.6 and, above
    # all, does NOT collapse to a single component (collapse would drive p1 to 0/1
    # and Omega to 0).  The estimate legitimately varies ~0.6-0.76 run to run, so the
    # band is generous -- the separation + non-collapse checks below carry the weight.
    expect_true(.p1 > 0.4 && .p1 < 0.85)
    # the two clearance groups stay separated (comp 1 low, comp 2 high)
    expect_true(.cl2 > 1.8 * .cl1)
    # Omega did not collapse
    expect_true(.fi$omega[1, 1] > 0.01)
  })

  test_that("IOV: an inter-occasion-variability model fits (BSV eta on a param without IOV)", {
    # This is the case that exercised the rxFromSE Subs/relational path: a between-
    # subject eta (eta.ka) on a parameter that does NOT carry the IOV (iov.cl).
    one.cmt <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.ka ~ 0.6
        eta.cl ~ 0.3
        iov.cl ~ 0.1 | occ
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl + iov.cl)
        v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_md
    .d$occ <- 1L
    .d$occ[.d$TIME >= 144] <- 2L
    .fi <- suppressWarnings(nlmixr2(one.cmt, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 20L, isample = 200L)))
    # the fit completes and reports the per-occasion IOV estimates
    expect_true("iov.cl" %in% names(.fi))
    expect_true(is.finite(.fi$objDf$OBJF[1]))
    # the structural population parameters come back finite and sensible
    expect_true(all(is.finite(fixef(.fi))))
    expect_true(fixef(.fi)["tcl"] > 0 && fixef(.fi)["tcl"] < 2)
  })

  test_that("General likelihood: a custom ll() model fits and matches FOCEI", {
    skip_if_not(rxode2hasLlik(), "rxode2 build has no llik support")
    # A general (custom) likelihood written out as EXACTLY the normal log-density,
    # so the reliable ground truth is the equivalent `add()` normal model (the
    # FOCEI outer optimizer does not itself estimate a non-`err` sigma living inside
    # a raw `ll()` -- a separate limitation).  impmap's non-mu M-step reads the
    # analytic d(ll)/d(theta) directly, so it estimates the raw-`ll()` add.sd and
    # must recover the same value as the equivalent normal fit.
    m <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.cl ~ 0.1
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        cp <- linCmt()
        ll(cp) ~ -0.5 * log(2 * pi) - log(add.sd) - 0.5 * ((DV - cp) / add.sd)^2
      })
    }
    mNorm <- function() {
      ini({
        tka <- 0.45; tcl <- 1; tv <- 3.45
        eta.cl ~ 0.1
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        cp <- linCmt(); cp ~ add(add.sd)
      })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(mNorm, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    .fi <- suppressWarnings(nlmixr2(m, .d, "impmap",
                                    impmapControl(print = 0L, nIter = 30L, isample = 300L)))
    expect_true(all(is.finite(fixef(.fi))))
    # the residual sigma and structural non-mu theta go through the M-step gradient
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
  })

  test_that("General likelihood: a + dnorm() (Laplace-form) endpoint fits", {
    skip_if_not(rxode2hasLlik(), "rxode2 build has no llik support")
    # dnorm() switches an otherwise-normal endpoint to the exact-likelihood
    # (Laplace) form; it drives the predMinusDv residual path where the fit's
    # foceiControl is reconstructed, so the impmap-internal M-step control fields
    # must be stripped from that down-conversion (else do.call(foceiControl, .)
    # errors "unused argument: impMuThetaIdx, ...").
    m <- function() {
      ini({ tka <- 0.45; tcl <- 1; tv <- 3.45; eta.cl ~ 0.1; add.sd <- 0.7 })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        linCmt() ~ add(add.sd) + dnorm()
      })
    }
    .fi <- suppressWarnings(nlmixr2(m, nlmixr2data::theo_sd, "impmap",
                                    impmapControl(print = 0L, nIter = 20L, isample = 200L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    expect_true("CWRES" %in% names(.fi))
    expect_true(all(is.finite(fixef(.fi))))
    # the down-converted foceiControl is valid (the failure mode was it not being one)
    expect_true(inherits(.fi$foceiControl, "foceiControl"))
  })

  test_that("est=\"imp\": importance-sampling EM without the MAP search", {
    # imp shares the impmap kernel but skips the per-iteration MAP search: the
    # proposal is centered at the running conditional mean with covariance gamma*V
    # (the previous conditional variance).  The mu-referenced, non-mu structural, and
    # residual-error thetas + Omega should all still converge to FOCEI.
    m <- function() {
      ini({ tka <- 0.45; tcl <- 1; tv <- 3.45; eta.cl ~ 0.1; add.sd <- 0.7 })
      model({ ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv); linCmt() ~ add(add.sd) })
    }
    .d <- nlmixr2data::theo_sd
    .ff <- suppressWarnings(nlmixr2(m, .d, "focei", foceiControl(print = 0L, covMethod = "")))
    .fi <- suppressWarnings(nlmixr2(m, .d, "imp",
                                    impControl(print = 0L, nIter = 40L, isample = 500L)))
    expect_true(inherits(.fi, "nlmixr2FitCore"))
    expect_equal(.fi$env$impmapControl$mapIter, 0L)   # imp == no MAP re-centering
    expect_equal(.fi$env$method, "imp")               # print header labels the method
    expect_true(all(is.finite(fixef(.fi))))
    # the non-mu structural theta + residual sigma move off their initial values and
    # match FOCEI (this is the path that needs the theta-sensitivity model built)
    expect_true(abs(unname(fixef(.fi)["add.sd"]) - 0.7) > 0.05)
    expect_equal(unname(fixef(.fi)["tcl"]), unname(fixef(.ff)["tcl"]), tolerance = 0.05)
    expect_equal(unname(fixef(.fi)["tv"]), unname(fixef(.ff)["tv"]), tolerance = 0.03)
    expect_equal(unname(fixef(.fi)["add.sd"]), unname(fixef(.ff)["add.sd"]), tolerance = 0.05)
    # the published objective is the converged FOCEi evaluation, on the same scale
    # as focei (a stale initial-MAP objective used to leak through the finalize and
    # get the n*log(2*pi) adjustment applied twice)
    expect_equal(.fi$objf, .ff$objf, tolerance = 0.01)
    expect_equal(as.numeric(logLik(.fi)), as.numeric(logLik(.ff)), tolerance = 0.01)
  })

  test_that("Bounds and fixed parameters are respected (theta bounds, fix(theta), fix(omega))", {
    # bounded tcl (log(c(lower, est, upper))), a fixed theta (tv), and a fixed
    # Omega diagonal (eta.cl) -- the EM M-step must hold the fixed values and keep
    # the bounded theta inside its limits.
    m <- function() {
      ini({
        tka <- 0.45
        tcl <- log(c(0, 2.7, 100))
        tv <- fix(3.45)
        eta.cl ~ fix(0.1)
        add.sd <- 0.7
      })
      model({
        ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv)
        linCmt() ~ add(add.sd)
      })
    }
    .fi <- suppressWarnings(nlmixr2(m, nlmixr2data::theo_sd, "impmap",
                                    impmapControl(print = 0L, nIter = 20L, isample = 200L)))
    # the fixed theta is held exactly
    expect_equal(unname(fixef(.fi)["tv"]), 3.45)
    expect_true(.fi$iniDf[.fi$iniDf$name == "tv", "fix"])
    # the fixed Omega diagonal is held exactly
    expect_equal(unname(.fi$omega[1, 1]), 0.1)
    # the bounded theta stayed inside (0, 100) on the natural scale
    expect_true(exp(unname(fixef(.fi)["tcl"])) > 0 && exp(unname(fixef(.fi)["tcl"])) < 100)
  })

  test_that("Parameter history is captured as the standard $parHist", {
    one.cmt <- function() {
      ini({ tka <- 0.45; tcl <- 1; tv <- 3.45; eta.cl ~ 0.1; add.sd <- 0.7 })
      model({ ka <- exp(tka); cl <- exp(tcl + eta.cl); v <- exp(tv); linCmt() ~ add(add.sd) })
    }
    .n <- 8L
    .fi <- suppressWarnings(nlmixr2(one.cmt, nlmixr2data::theo_sd, "impmap",
                                    impmapControl(print = 0L, nIter = .n, isample = 150L)))
    .ph <- .fi$parHist
    expect_true(is.data.frame(.ph))
    # one row per EM iteration, standard iter/objf columns + the estimated params
    expect_equal(nrow(.ph), .n)
    expect_true(all(c("iter", "objf", "tka", "tcl", "tv", "add.sd") %in% names(.ph)))
    expect_equal(.ph$iter, seq_len(.n))
    # the objective and the estimates are finite and moving toward convergence
    expect_true(all(is.finite(.ph$objf)))
    expect_true(.ph$objf[.n] < .ph$objf[1])
    # the EM walk is unscaled, so no redundant "Unscaled" rows and no gradient rows
    .types <- table(.fi$parHistData$type)
    expect_equal(unname(.types[["Unscaled"]]), 0L)
    expect_true(all(.types[grepl("Gradient|Difference|Sensitivity", names(.types))] == 0L))
    # the print header identifies the method (not FOCEi's outer-optimizer label)
    expect_equal(.fi$env$method, "impmap")
    expect_equal(.fi$env$extra, "")
  })
})