tests/testthat/testTheo.R

In tdmore-dev/tdmore: In silico evaluation of precision dosing

library(tdmore)
library(nlmixr)
library(ggplot2)
library(testthat)

asPhi <- function(estimate) {
  par <- coef(estimate)
  df1 <- as.data.frame( t(par) )
  colnames(df1) <- paste0("ETA(", seq_along(par), ")")

  vcov <- vcov(estimate)
  n <- nrow(vcov)
  result <- c()
  for(i in 1:n) {
    for(j in 1:i) {
      value <- vcov[i, j]
      names(value) <- paste0("ETC(",i,",",j,")")
      result <- c(result, value)
    }
  }
  df2 <- data.frame( t(result), check.names = F)

  cbind(df1, df2)
}

## initial values are FOCEi results for estimation of THEO dataset (with fixed dose of 320)
modelCode <- function() {
  ini({
    THETA1 <- 1.49072E+00
    THETA2 <- 3.24470E+01
    THETA3 <- 8.72901E-02
    ETA1 + ETA2 + ETA3 ~ c(4.34717E-01,  5.69202E-02,  1.93995E-02, -6.48548E-03,  1.20830E-02,  2.02334E-02)
    EPS1 <- 0.1305320
    EPS2 <- 0.2878263
  })
  model({
    KA=THETA1*exp(ETA1)
    V=THETA2*exp(ETA2)
    K=THETA3*exp(ETA3)

    d/dt(abs) = -KA*abs
    d/dt(centr) = KA*abs - K*centr
    #TODO: we do not support algebraic equations?
    #CONC=DOSE / V * KA/(KA-K)*(exp(-K*t) - exp(-KA*t))
    CONC=centr/V
    CONC ~ prop(EPS1) + add(EPS2)
  })
}

m1 <- nlmixrUI(modelCode)
m1tdm <- tdmore(m1, atol=1e-12, rtol=1e-12)
regimen <- data.frame(TIME=0, AMT=320)

## THEO dataset from NONMEM distribution
input <- data.frame(
  ID = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
         2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
         3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L,
         5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
         6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L,
         8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L,
         9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
         10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
         12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L),
  AMT = c(4.02, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 4.4, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, 4.53, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, 4.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 5.86,
          NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 4, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, 4.95, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          4.53, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.1, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, 5.5, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, 4.92, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 5.3, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, NA),
  TIME = c(0, 0.25, 0.57, 1.12, 2.02, 3.82, 5.1, 7.03, 9.05, 12.12,
           24.37, 0, 0.27, 0.52, 1, 1.92, 3.5, 5.02, 7.03, 9, 12, 24.3,
           0, 0.27, 0.58, 1.02, 2.02, 3.62, 5.08, 7.07, 9, 12.15, 24.17, 0,
           0.35, 0.6, 1.07, 2.13, 3.5, 5.02, 7.02, 9.02, 11.98, 24.65, 0,
           0.3, 0.52, 1, 2.02, 3.5, 5.02, 7.02, 9.1, 12, 24.35, 0, 0.27,
           0.58, 1.15, 2.03, 3.57, 5, 7, 9.22, 12.1, 23.85, 0, 0.25, 0.5,
           1.02, 2.02, 3.48, 5, 6.98, 9, 12.05, 24.22, 0, 0.25, 0.52, 0.98,
           2.02, 3.53, 5.05, 7.15, 9.07, 12.1, 24.12, 0, 0.3, 0.63, 1.05,
           2.02, 3.53, 5.02, 7.17, 8.8, 11.6, 24.43, 0, 0.37, 0.77, 1.02, 2.05,
           3.55, 5.05, 7.08, 9.38, 12.1, 23.7, 0, 0.25, 0.5, 0.98, 1.98,
           3.6, 5.02, 7.03, 9.03, 12.12, 24.08, 0, 0.25, 0.5, 1, 2, 3.52,
           5.07, 7.07, 9.03, 12.05, 24.15),
  DV = c(0.74, 2.84, 6.57, 10.5, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94,
         3.28, 0, 1.72, 7.91, 8.31, 8.33, 6.85, 6.08, 5.4, 4.55, 3.01,
         0.9, 0, 4.4, 6.9, 8.2, 7.8, 7.5, 6.2, 5.3, 4.9, 3.7, 1.05, 0,
         1.89, 4.6, 8.6, 8.38, 7.54, 6.88, 5.78, 5.33, 4.19, 1.15, 0, 2.02,
         5.63, 11.4, 9.33, 8.74, 7.56, 7.09, 5.9, 4.37, 1.57, 0, 1.29,
         3.08, 6.44, 6.32, 5.53, 4.94, 4.02, 3.46, 2.78, 0.92, 0.15, 0.85,
         2.35, 5.02, 6.58, 7.09, 6.66, 5.25, 4.39, 3.53, 1.15, 0, 3.05,
         3.05, 7.31, 7.56, 6.59, 5.88, 4.73, 4.57, 3, 1.25, 0, 7.37,
         9.03, 7.14, 6.33, 5.66, 5.67, 4.24, 4.11, 3.16, 1.12, 0.24, 2.89,
         5.22, 6.41, 7.83, 10.21, 9.18, 8.02, 7.14, 5.68, 2.42, 0, 4.86,
         7.24, 8, 6.81, 5.87, 5.22, 4.45, 3.62, 2.69, 0.86, 0, 1.25, 3.96,
         7.82, 9.72, 9.75, 8.57, 6.59, 6.11, 4.57, 1.17),
  BWT = c(79.6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 72.4, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, 70.5, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, 72.7, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          54.6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 80, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, 64.6, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, 70.5, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 86.4, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, NA, 58.2, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, 65, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          60.5, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)
)

library(dplyr)
expect_runtime( {
  results_m1 <- input %>% group_by(ID) %>% do({
    observed <- .data
    obs <- observed[, c("TIME", "DV")]
    colnames(obs) <- c("TIME", "CONC")
    est <- estimate(m1tdm, obs, regimen)

    cbind(asPhi(est),
          data.frame(
            ID=observed$ID[1],
            OBJ=est$ofv
          ))
  }) %>% ungroup() %>% as.data.frame()
}, "estimate_theo.runtime.txt", upper=2.5)
#2.5 times is OK; to avoid errors in Docker Cloud

expected_results_m1 <- data.frame(
  `ETA(1)` = c(-0.104660664215003, 0.329025944024965, 0.411113690829706,
             -0.34417129545638, -0.0303059358208419, -0.471091905410399,
             -0.858244023045384, -0.0682371600428652, 1.33496175489568,
             -0.71317160561277, 0.865596068010048, -0.496172377531442),
  `ETA(2)` = c(-0.192539549970711, 0.0563206253944116, 0.0439044137447579,
             -0.0420778311793217, -0.108824692232917, 0.13798880959419,
             0.00291914243084671, 0.0708294139509278, 0.177451612093571,
             -0.214528821024506, 0.176336899728112, -0.134716528736458),
  `ETA(3)` = c(-0.279385491469617, 0.0440627231614929, -0.0129154414961791,
             0.0156405168994427, -0.0957211045422365, 0.232775943204004,
             0.146435435807427, 0.0792041339649327, -0.0126039120480183,
             -0.149184709693789, 0.0777919267358549, -0.029206882488009),
  `ETC(1,1)` = c(0.0210834141442305, 0.0287458887412714, 0.0298933751973137,
               0.019480143631608, 0.020703773520106, 0.021957863250446,
               0.0168548046713357, 0.0247346001427467, 0.068259135317654,
               0.0153118558270015, 0.0403415142765707, 0.0149459623625841),
  `ETC(2,1)` = c(0.00366572693314466, 0.00411898319533273, 0.00424709825254382,
               0.00354080102638969, 0.0033994255821157, 0.00414907513168687,
               0.00355442832126625, 0.00396088656906548, 0.00616909249653027,
               0.00310710577078467, 0.00503929584092165, 0.00285679898245941),
  `ETC(2,2)` = c(0.00190834017311547, 0.00189919169128497, 0.00190660480646897,
               0.00195164892426712, 0.00178642865769596, 0.00230468978730928,
               0.0021954984020324, 0.00206583339205648, 0.00189046725040563,
               0.00188226071970371, 0.00198631899649564, 0.00177270548882738),
  `ETC(3,1)` = c(-0.00205399844793196, -0.00285799735076052,
               -0.00301228511143907, -0.0021037969206494,
               -0.00214511418541834, -0.00214015237456608, -0.00167023723239374,
               -0.00250563915779664, -0.00604430914896841, -0.00180767831066223,
               -0.00369438417115206, -0.00173649804932396),
  `ETC(3,2)` = c(-0.000424878877753965, -0.000449896881796006,
               -0.000505294722018099, -0.000434323023749888,
               -0.000412883531992041, -0.000314851566735198, -0.000297007127463613,
               -0.000388750719403323, -0.000531060443697142, -0.000444448641386078,
               -0.000462263674984931, -0.000408881433073268),
  `ETC(3,3)` = c(0.00303833864865961, 0.00344265672464375, 0.00333370435936775,
               0.00328374877428826, 0.00333994994172282, 0.00351935833951445,
               0.00345425473899878, 0.00342739505745427, 0.00382690833031373,
               0.00322080957437884, 0.00356028530018839, 0.00335456279096149),
  ID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L),
  OBJ = c(29.1410690800101, 27.5313632957081, 8.5609851062512,
          18.7811120410482, 29.8966881792101, 13.8129616612949,
          10.6676450568203, 13.6414107810644, 14.3392116396854, 13.8981508342129,
          7.94278603869878, 20.9078897884101),
  check.names=F
)

expect_equal(results_m1, expected_results_m1, tolerance=1E-5)