tests/testthat/helper-wang2007.R

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

# Helper functions for Wang2007 test suite
# These functions support the split test files for parallel execution

################################################################################
# Data Setup Functions
################################################################################

# Function to setup Wang2007 base dataset
getWang2007BaseData <- function() {
  dat <- Wang2007
  dat$DV <- dat$Y
  dat
}

# Function to setup Wang2007 dataset with dosing (dat2)
getWang2007DoseData <- function() {
  dat <- Wang2007
  dat$DV <- dat$Y
  dat2 <- dat[dat$Time == 0, ]
  dat2$EVID <- 101
  dat2$AMT <- 10
  dat2 <- rbind(dat2, data.frame(dat, EVID = 0, AMT = 0))
  dat2 <- dat2[(order(dat2$ID, -dat2$EVID, dat2$Time)), ]
  dat2
}

# Function to setup log-transformed data
getWang2007LogData <- function() {
  dat <- Wang2007
  dat$DV <- dat$Y
  datl <- dat
  datl$DV <- log(datl$DV)
  datl
}

# Function to setup log-transformed data with dosing
getWang2007LogDoseData <- function() {
  dat2 <- getWang2007DoseData()
  datl2 <- dat2
  datl2$DV <- log(datl2$DV)
  datl2
}

################################################################################
# Model Definition Functions
################################################################################

# Base PK model without ODE
getWang2007BaseModel <- function() {
  f <- function() {
    ini({
      tke <- 0.5
      eta.ke ~ 0.04
      prop.sd <- sqrt(0.1)
    })
    model({
      ke <- tke * exp(eta.ke)
      ipre <- 10 * exp(-ke * t)
      f2 <- ipre / (ipre + 5)
      f3 <- f2 * 3
      lipre <- log(ipre)
      ipre ~ prop(prop.sd)
    })
  }
  .nlmixr(f)
}

# Base PK model with ODE
getWang2007OdeModel <- function() {
  fo <- function() {
    ini({
      tke <- 0.5
      eta.ke ~ 0.04
      prop.sd <- sqrt(0.1)
    })
    model({
      ke <- tke * exp(eta.ke)
      d/dt(ipre) <- -ke * ipre
      f2 <- ipre / (ipre + 5)
      f3 <- f2 * 3
      lipre <- log(ipre)
      ipre ~ prop(prop.sd)
    })
  }
  .nlmixr(fo)
}

################################################################################
# Core Testing Function
################################################################################

testWang2007ErrorModel <- function(type, fun, val = rep(NA_real_, 10), addProp = 2, log=FALSE) {
  ## message(type)
  valName <- as.character(substitute(val))

  # Get base models
  f <- getWang2007BaseModel()
  fo <- getWang2007OdeModel()

  .f <- fun(f)
  .fo <- fun(fo)
  .dode <- getWang2007DoseData()
  .d <- getWang2007BaseData()
  if (log) {
    .d <- getWang2007LogData()
    .dode <- getWang2007LogDoseData()
  }

  # Only the objective values (and that the fit prints) are checked below, so
  # skip residual/table computation (calcTables = FALSE).  This leaves the
  # objective identical to NONMEM while cutting each fit's time by ~10x.
  .combined <- paste0("combined", addProp)

  fit1 <- .nlmixr(.fo, .dode, "focei",
                  control = foceiControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit2 <- .nlmixr(.f, .d, "focei",
                  control = foceiControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit3 <- .nlmixr(.fo, .dode, "foce",
                  control = foceiControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit4 <- .nlmixr(.f, .d, "foce",
                  control = foceiControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit5 <- .nlmixr(.fo, .dode, "fo",
                  control = foceiControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))
  setOfv(fit5, "fo")

  fit6 <- .nlmixr(.f, .dode, "fo", control = foceiControl(
    maxOuterIterations = 0, covMethod = "",
    addProp = .combined, calcTables = FALSE))
  setOfv(fit6, "fo")

  fit7 <- .nlmixr(.fo, .dode, "agq",
                  control = agqControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit8 <- .nlmixr(.f, .d, "agq",
                  control = agqControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit9 <- .nlmixr(.fo, .dode, "laplace",
                  control = laplaceControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  fit10 <- .nlmixr(.f, .d, "laplace",
                  control = laplaceControl(
                    maxOuterIterations = 0, covMethod = "",
                    addProp = .combined, calcTables = FALSE))

  .n <- paste(type, c("focei ode", "focei", "foce ode", "foce", "fo ode", "fo",
                      "agq ode", "agq", "laplace ode", "laplace"),
              paste0("combined", addProp))
  ret <- c(fit1$objective, fit2$objective, fit3$objective, fit4$objective,
           fit5$objective, fit6$objective, fit7$objective, fit8$objective,
           fit9$objective, fit10$objective)
  ret <- setNames(ret, .n)
  if (length(val) < length(.n)) {
    val <- c(val, rep(NA_real_, length(.n) - length(val)))
  }
  val <- setNames(val, .n)

  ## if (!identical(round(val, 3), round(ret, 3))) {
  ##   print("\n")
  ##   t <- try(print(str2lang(paste0(valName, " <- ", deparse1(round(setNames(ret, NULL), 3))))), silent=TRUE)
  ##   if (inherits(t, "try-error")) {
  ##     print(type)
  ##     try(print(str2lang(deparse1(round(setNames(ret, NULL), 3)))))
  ##   }
  ## }

  ## Now test
  if (!all(is.na(val))) {
    test_that(
      type,{
        expect_equal(setNames(ret, NULL), setNames(val, NULL), tolerance=1e-3)
      })
    test_that(paste0(type, " print"),{
      utils::capture.output(suppressMessages({
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit1), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit1), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit2), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit2), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit3), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit3), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit4), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit4), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit5), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit5), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit6), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit6), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit7), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit7), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit8), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit8), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit9), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit9), NA))
        withr::with_options(list(cli.unicode=TRUE),expect_error(print(fit10), NA))
        withr::with_options(list(cli.unicode=FALSE),expect_error(print(fit10), NA))
      }))
    })
  }
  invisible(ret)
}