tests/testthat/test-nlm-cens.R

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

nmTest({
  # Tests for censoring support in the NLM estimation engine
  # These tests verify that NLM handles M2, M3, M4 censoring
  # matching the behavior used in FOCEI/SAEM

  one.cmt <- function() {
    ini({
      tka <- 0.45
      tcl <- log(c(0, 2.7, 100))
      tv <- 3.45
      add.sd <- 0.7
    })
    model({
      ka <- exp(tka)
      cl <- exp(tcl)
      v <- exp(tv)
      linCmt() ~ add(add.sd)
    })
  }

  # Base dataset (no censoring)
  .dat <- nlmixr2data::theo_sd

  # censOption is inert for NLM (its finite-difference outer Hessian already reflects
  # censoring), so the NLM censoring text stays PLAIN (no " (laplace)"/" (gauss)" suffix --
  # that is FOCEI/FOCE only).  The strip is a defensive no-op that keeps the checks robust.
  .censMethod <- function(x) sub(" \\((laplace|gauss)\\)$", "", as.character(x$censInformation))

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " works without censoring (no CENS column)"), {
      expect_s3_class(fit, paste0("nlmixr2.", meth))
      expect_equal(as.character(fit$censInformation), "No censoring")
    })
  }

  .dat0 <- .dat
  .dat0$CENS <- 0L

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit0 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat0, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit0 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat0, est = meth, list(print = 0))
      ))
    }
    test_that("nlm & related methods works with CENS column all zeros", {
      expect_s3_class(fit0, paste0("nlmixr2.", meth))
      expect_equal(as.character(fit0$censInformation), "No censoring")
    })
  }

  # Create censored datasets
  # Mark observations below a threshold as left-censored (M3)
  .LLOQ <- 2.0
  .datM3 <- .dat
  .datM3$CENS <- ifelse(.datM3$DV < .LLOQ & .datM3$EVID == 0, 1L, 0L)
  # For censored obs, set DV to LLOQ
  .datM3$DV[.datM3$CENS == 1] <- .LLOQ

  test_that("nls does not support censoring", {
    expect_error(
      suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = "nls", list(print = 0)))))
  })

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " accepts and processes M3 (left) censored data"), {
      expect_s3_class(fit_m3, paste0("nlmixr2.", meth))
      expect_equal(.censMethod(fit_m3), "M3 censoring")
      # NLM censoring text is plain (censOption inert) -- no laplace/gauss suffix
      expect_equal(as.character(fit_m3$censInformation), "M3 censoring")
    })
  }


  test_that("nls does not support censoring", {
    expect_error(
      suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = "nls", list(print = 0))))
    )
  })

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_base <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0, method="BFGS"))
      ))
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0, method="BFGS"))
      ))
      expect_false(isTRUE(all.equal(fit_base$objf, fit_m3$objf)))
    } else {
      fit_base <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0))
      ))
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0))
      ))
      test_that(paste0(meth, "M3 censoring changes the objective function vs no censoring"), {
        expect_false(isTRUE(all.equal(fit_base$objf, fit_m3$objf)))
      })
    }
  }

  # M2 censoring: CENS=0 with a finite LIMIT
  .datM2 <- .dat
  .datM2$CENS <- 0L
  .datM2$LIMIT <- 0  # interval censoring: all obs have a lower bound of 0

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_m2 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM2, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit_m2 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM2, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " accepts and processes M2 (interval) censored data"), {
      expect_s3_class(fit_m2, paste0("nlmixr2.", meth))
      expect_equal(.censMethod(fit_m2), "M2 censoring")
    })
  }

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_base <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0, method="BFGS"))
      ))
      fit_m2 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM2, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit_base <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .dat, est = meth, list(print = 0))
      ))
      fit_m2 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM2, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " method: M2 censoring changes the objective function vs no censoring"), {
      expect_false(isTRUE(all.equal(fit_base$objf, fit_m2$objf)))
    })
  }

  # M4 censoring: CENS!=0 with a finite LIMIT
  .datM4 <- .datM3
  .datM4$LIMIT <- 0  # add LIMIT for M4

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_m4 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM4, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit_m4 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM4, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " method accepts and processes M4 (interval-censored) data"), {
      expect_s3_class(fit_m4, paste0("nlmixr2.", meth))
      expect_equal(.censMethod(fit_m4), "M2 and M4 censoring")
    })
  }

  for (meth in c("nlm", "bobyqa", "lbfgsb3c", "n1qn1", "newuoa", "nlminb", "optim")) {
    if (meth == "optim") {
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0, method="BFGS"))
      ))
      fit_m4 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM4, est = meth, list(print = 0, method="BFGS"))
      ))
    } else {
      fit_m3 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM3, est = meth, list(print = 0))
      ))
      fit_m4 <- suppressMessages(suppressWarnings(
        .nlmixr(one.cmt, .datM4, est = meth, list(print = 0))
      ))
    }
    test_that(paste0(meth, " method: M3 and M4 give different results"), {
      expect_false(isTRUE(all.equal(fit_m3$objf, fit_m4$objf)))
    })
  }

})