focei estimation fails with max() and abs()

[billdenney]

I'm wanting to estimate the continuous number of transit compartments for an oncology model. But, I get an error when trying to estimate (that is not a problem when simulating). I'm pretty sure that my hat function with max() and abs() is the problem.

library(nlmixr2)
#> Loading required package: nlmixr2data

oncology_sdm_lobo_2002 <- function() {
  description <- "Signal transduction model for delayed concentration effects on cancer cell growth"
  reference <- "Lobo ED, Balthasar JP. Pharmacodynamic modeling of chemotherapeutic effects: Application of a transit compartment model to characterize methotrexate effects in vitro. AAPS J. 2002;4(4):212-222. doi:10.1208/ps040442"
  # Values for lkng, ltau, lec50, and kmax are for methotrexate from Lobo 2002,
  # Table 2.  propErr and addErr are added as reasonable values though not from
  # Lobo 2002 where no value is apparent in the paper.
  ini({
    lkng <- log(0.02) ; label("Cell net growth rate (growth minus death) (1/hr)")
    ltau <- log(c(1, 34.1, 500)) ; label("Mean transit time of each transit compartment (hr)")
    lec50 <- log(c(1, 50, 2000)) ; label("Drug concentration reducing the cell growth by 50% (ug/mL)")
    kmax <- 0.01 ; label("Maximum drug-related reduction in cell growth (1/hr)")
    ncmt <- c(1, 4, 10) ; label("Number of transit compartments")

    propErr <- c(0, 0.3) ; label("Proportional residual error (fraction)")
    addErr <- c(0, 50, 1000) ; label("Additive residual error (tumor volume units)")
  })
  model({
    # cp is the drug concentration
    kng <- exp(lkng)
    tau <- exp(ltau)
    taulast <- tau
    ec50 <- exp(lec50)

    edrug <- kmax*cp/(ec50 + cp)

    tumor(0) <- tumor0
    d/dt(transit1) <- (edrug - transit1)/tau
    d/dt(transit2) <- (transit1 - transit2)/tau
    d/dt(transit3) <- (transit2 - transit3)/tau
    d/dt(transit4) <- (transit3 - transit4)/tau
    d/dt(transit5) <- (transit4 - transit5)/tau
    d/dt(transit6) <- (transit5 - transit6)/tau
    d/dt(transit7) <- (transit6 - transit7)/tau
    d/dt(transit8) <- (transit7 - transit8)/tau
    d/dt(transit9) <- (transit8 - transit9)/tau
    d/dt(transit10) <- (transit9 - transit10)/tau
    frac1 <- max(1 - abs(ncmt - 1), 0)
    frac2 <- max(1 - abs(ncmt - 2), 0)
    frac3 <- max(1 - abs(ncmt - 3), 0)
    frac4 <- max(1 - abs(ncmt - 4), 0)
    frac5 <- max(1 - abs(ncmt - 5), 0)
    frac6 <- max(1 - abs(ncmt - 6), 0)
    frac7 <- max(1 - abs(ncmt - 7), 0)
    frac8 <- max(1 - abs(ncmt - 8), 0)
    frac9 <- max(1 - abs(ncmt - 9), 0)
    frac10 <- max(1 - abs(ncmt - 10), 0)
    drug_effect <- frac1*transit1 + frac2*transit2 + frac3*transit3 + frac4*transit4 + frac5*transit5 + frac6*transit6 + frac7*transit7 + frac8*transit8 + frac9*transit9 + frac10*transit10
    d/dt(tumor) <- kng*tumor - drug_effect*tumor
    tumor ~ prop(propErr) + add(addErr)
  })
}

nlmixr(oncology_sdm_lobo_2002)
#>  ── rxode2-based free-form 11-cmt ODE model ───────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>      lkng      ltau     lec50      kmax      ncmt   propErr    addErr 
#> -3.912023  3.529297  3.912023  0.010000  4.000000  0.300000 50.000000 
#> 
#> States ($state or $stateDf): 
#>    Compartment Number Compartment Name
#> 1                   1         transit1
#> 2                   2         transit2
#> 3                   3         transit3
#> 4                   4         transit4
#> 5                   5         transit5
#> 6                   6         transit6
#> 7                   7         transit7
#> 8                   8         transit8
#> 9                   9         transit9
#> 10                 10        transit10
#> 11                 11            tumor
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     description <- "Signal transduction model for delayed concentration effects on cancer cell growth"
#>     reference <- "Lobo ED, Balthasar JP. Pharmacodynamic modeling of chemotherapeutic effects: Application of a transit compartment model to characterize methotrexate effects in vitro. AAPS J. 2002;4(4):212-222. doi:10.1208/ps040442"
#>     ini({
#>         lkng <- -3.91202300542815
#>         label("Cell net growth rate (growth minus death) (1/hr)")
#>         ltau <- c(0, 3.52929738428947, 6.21460809842219)
#>         label("Mean transit time of each transit compartment (hr)")
#>         lec50 <- c(0, 3.91202300542815, 7.60090245954208)
#>         label("Drug concentration reducing the cell growth by 50% (ug/mL)")
#>         kmax <- 0.01
#>         label("Maximum drug-related reduction in cell growth (1/hr)")
#>         ncmt <- c(1, 4, 10)
#>         label("Number of transit compartments")
#>         propErr <- c(0, 0.3)
#>         label("Proportional residual error (fraction)")
#>         addErr <- c(0, 50, 1000)
#>         label("Additive residual error (tumor volume units)")
#>     })
#>     model({
#>         kng <- exp(lkng)
#>         tau <- exp(ltau)
#>         taulast <- tau
#>         ec50 <- exp(lec50)
#>         edrug <- kmax * cp/(ec50 + cp)
#>         tumor(0) <- tumor0
#>         d/dt(transit1) <- (edrug - transit1)/tau
#>         d/dt(transit2) <- (transit1 - transit2)/tau
#>         d/dt(transit3) <- (transit2 - transit3)/tau
#>         d/dt(transit4) <- (transit3 - transit4)/tau
#>         d/dt(transit5) <- (transit4 - transit5)/tau
#>         d/dt(transit6) <- (transit5 - transit6)/tau
#>         d/dt(transit7) <- (transit6 - transit7)/tau
#>         d/dt(transit8) <- (transit7 - transit8)/tau
#>         d/dt(transit9) <- (transit8 - transit9)/tau
#>         d/dt(transit10) <- (transit9 - transit10)/tau
#>         frac1 <- max(1 - abs(ncmt - 1), 0)
#>         frac2 <- max(1 - abs(ncmt - 2), 0)
#>         frac3 <- max(1 - abs(ncmt - 3), 0)
#>         frac4 <- max(1 - abs(ncmt - 4), 0)
#>         frac5 <- max(1 - abs(ncmt - 5), 0)
#>         frac6 <- max(1 - abs(ncmt - 6), 0)
#>         frac7 <- max(1 - abs(ncmt - 7), 0)
#>         frac8 <- max(1 - abs(ncmt - 8), 0)
#>         frac9 <- max(1 - abs(ncmt - 9), 0)
#>         frac10 <- max(1 - abs(ncmt - 10), 0)
#>         drug_effect <- frac1 * transit1 + frac2 * transit2 + 
#>             frac3 * transit3 + frac4 * transit4 + frac5 * transit5 + 
#>             frac6 * transit6 + frac7 * transit7 + frac8 * transit8 + 
#>             frac9 * transit9 + frac10 * transit10
#>         d/dt(tumor) <- kng * tumor - drug_effect * tumor
#>         tumor ~ prop(propErr) + add(addErr)
#>     })
#> }

d_sim <- data.frame(id=1, time=0:24, tumor0=100, cp=1:25, dv=101:125)
# Simulation works
nlmixr(oncology_sdm_lobo_2002, data = d_sim, est="rxSolve")
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>      lkng      ltau     lec50      kmax      ncmt   propErr    addErr 
#> -3.912023  3.529297  3.912023  0.010000  4.000000  0.300000 50.000000 
#> ── Initial Conditions (x$inits): ──
#>  transit1  transit2  transit3  transit4  transit5  transit6  transit7  transit8 
#>         0         0         0         0         0         0         0         0 
#>  transit9 transit10     tumor 
#>         0         0         0 
#> ── First part of data (object): ──
#> # A tibble: 25 × 32
#>    time   kng   tau taulast  ec50    edrug frac1 frac2 frac3 frac4 frac5 frac6
#>   <dbl> <dbl> <dbl>   <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1     0  0.02  34.1    34.1    50 0.000196     0     0     0     1     0     0
#> 2     1  0.02  34.1    34.1    50 0.000385     0     0     0     1     0     0
#> 3     2  0.02  34.1    34.1    50 0.000566     0     0     0     1     0     0
#> 4     3  0.02  34.1    34.1    50 0.000741     0     0     0     1     0     0
#> 5     4  0.02  34.1    34.1    50 0.000909     0     0     0     1     0     0
#> 6     5  0.02  34.1    34.1    50 0.00107      0     0     0     1     0     0
#> # … with 19 more rows, and 20 more variables: frac7 <dbl>, frac8 <dbl>,
#> #   frac9 <dbl>, frac10 <dbl>, drug_effect <dbl>, ipredSim <dbl>, sim <dbl>,
#> #   transit1 <dbl>, transit2 <dbl>, transit3 <dbl>, transit4 <dbl>,
#> #   transit5 <dbl>, transit6 <dbl>, transit7 <dbl>, transit8 <dbl>,
#> #   transit9 <dbl>, transit10 <dbl>, tumor <dbl>, cp <dbl>, tumor0 <dbl>
# Estimation fails
nlmixr(oncology_sdm_lobo_2002, data = d_sim, est="focei")
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> Error in drug_effect * tumor : non-numeric argument to binary operator
#> Error: non-numeric argument to binary operator

^{Created on 2022-10-18 with reprex v2.0.2}

[billdenney]

I tracked it this far: https://github.com/nlmixr2/nlmixr2est/blob/dd5409f7fe467e6c9b4fbf17166171cc57da759d/R/nlmixr2Est.R#L175

[billdenney]

Further digging takes me here:

https://github.com/nlmixr2/nlmixr2est/blob/dd5409f7fe467e6c9b4fbf17166171cc57da759d/R/focei.R#L437

When I look there, I see that all of the frac variables (e.g. frac1) have the class of "try-error", and the error text is "Error in max(1 - abs(ncmt - 9), 0) : invalid 'type' (S4) of argument\n".

Are the max() and abs() functions not allowed? If they aren't, can a clearer error message be given about unsupported functions?

[billdenney]

The model progressed (with many complaints about integration) by switching to if-like coding:

oncology_sdm_lobo_2002 <- function() {
  description <- "Signal transduction model for delayed concentration effects on cancer cell growth"
  reference <- "Lobo ED, Balthasar JP. Pharmacodynamic modeling of chemotherapeutic effects: Application of a transit compartment model to characterize methotrexate effects in vitro. AAPS J. 2002;4(4):212-222. doi:10.1208/ps040442"
  # Values for lkng, ltau, lec50, and kmax are for methotrexate from Lobo 2002,
  # Table 2.  propErr and addErr are added as reasonable values though not from
  # Lobo 2002 where no value is apparent in the paper.
  ini({
    lkng <- log(0.02) ; label("Cell net growth rate (growth minus death) (1/hr)")
    ltau <- log(c(1, 34.1, 500)) ; label("Mean transit time of each transit compartment (hr)")
    lec50 <- log(c(1, 50, 2000)) ; label("Drug concentration reducing the cell growth by 50% (ug/mL)")
    kmax <- 0.01 ; label("Maximum drug-related reduction in cell growth (1/hr)")
    ncmt <- c(1, 4, 10) ; label("Number of transit compartments")

    propErr <- c(0, 0.3) ; label("Proportional residual error (fraction)")
    addErr <- c(0, 50, 1000) ; label("Additive residual error (tumor volume units)")
  })
  model({
    # cp is the drug concentration
    kng <- exp(lkng)
    tau <- exp(ltau)
    taulast <- tau
    ec50 <- exp(lec50)

    edrug <- kmax*cp/(ec50 + cp)

    tumor(initial) <- tumor0
    d/dt(transit1) <- (edrug - transit1)/tau
    d/dt(transit2) <- (transit1 - transit2)/tau
    d/dt(transit3) <- (transit2 - transit3)/tau
    d/dt(transit4) <- (transit3 - transit4)/tau
    d/dt(transit5) <- (transit4 - transit5)/tau
    d/dt(transit6) <- (transit5 - transit6)/tau
    d/dt(transit7) <- (transit6 - transit7)/tau
    d/dt(transit8) <- (transit7 - transit8)/tau
    d/dt(transit9) <- (transit8 - transit9)/tau
    d/dt(transit10) <- (transit9 - transit10)/tau
    frac1 <- (1 - abs(ncmt - 1))*(ncmt < 2)
    frac2 <- (1 - abs(ncmt - 2))*(ncmt > 1 & ncmt < 3)
    frac3 <- (1 - abs(ncmt - 3))*(ncmt > 2 & ncmt < 4)
    frac4 <- (1 - abs(ncmt - 4))*(ncmt > 3 & ncmt < 5)
    frac5 <- (1 - abs(ncmt - 5))*(ncmt > 4 & ncmt < 6)
    frac6 <- (1 - abs(ncmt - 6))*(ncmt > 5 & ncmt < 7)
    frac7 <- (1 - abs(ncmt - 7))*(ncmt > 6 & ncmt < 8)
    frac8 <- (1 - abs(ncmt - 8))*(ncmt > 7 & ncmt < 9)
    frac9 <- (1 - abs(ncmt - 9))*(ncmt > 8 & ncmt < 10)
    frac10 <- (1 - abs(ncmt - 10))*(ncmt > 9)
    drug_effect <- frac1*transit1 + frac2*transit2 + frac3*transit3 + frac4*transit4 + frac5*transit5 + frac6*transit6 + frac7*transit7 + frac8*transit8 + frac9*transit9 + frac10*transit10
    d/dt(tumor) <- kng*tumor - drug_effect*tumor
    tumor ~ prop(propErr) + add(addErr)
  })
}

[mattfidler]

In theory they are allowed, though they are a bit finicky when it comes to derivatives.

Perhaps they should be disallowed.

[mattfidler]

Change to logical operations, example:

https://stackoverflow.com/questions/30738923/max-implemented-with-basic-operators

[mattfidler]

It will work on focei with and without etas

library(nlmixr2est)
#> Loading required package: nlmixr2data

oncology_sdm_lobo_2002 <- function() {
  description <- "Signal transduction model for delayed concentration effects on cancer cell growth"
  reference <- "Lobo ED, Balthasar JP. Pharmacodynamic modeling of chemotherapeutic effects: Application of a transit compartment model to characterize methotrexate effects in vitro. AAPS J. 2002;4(4):212-222. doi:10.1208/ps040442"
  # Values for lkng, ltau, lec50, and kmax are for methotrexate from Lobo 2002,
  # Table 2.  propErr and addErr are added as reasonable values though not from
  # Lobo 2002 where no value is apparent in the paper.
  ini({
    lkng <- log(0.02) ; label("Cell net growth rate (growth minus death) (1/hr)")
    ltau <- log(c(1, 34.1, 500)) ; label("Mean transit time of each transit compartment (hr)")
    lec50 <- log(c(1, 50, 2000)) ; label("Drug concentration reducing the cell growth by 50% (ug/mL)")
    kmax <- 0.01 ; label("Maximum drug-related reduction in cell growth (1/hr)")
    ncmt <- c(1, 4, 10) ; label("Number of transit compartments")

    propErr <- c(0, 0.3) ; label("Proportional residual error (fraction)")
    eta.kng ~ 1
    addErr <- c(0, 50, 1000) ; label("Additive residual error (tumor volume units)")
  })
  model({
    # cp is the drug concentration
    kng <- exp(lkng + eta.kng)
    tau <- exp(ltau)
    taulast <- tau
    ec50 <- exp(lec50)

    edrug <- kmax*cp/(ec50 + cp)

    tumor(0) <- tumor0
    d/dt(transit1) <- (edrug - transit1)/tau
    d/dt(transit2) <- (transit1 - transit2)/tau
    d/dt(transit3) <- (transit2 - transit3)/tau
    d/dt(transit4) <- (transit3 - transit4)/tau
    d/dt(transit5) <- (transit4 - transit5)/tau
    d/dt(transit6) <- (transit5 - transit6)/tau
    d/dt(transit7) <- (transit6 - transit7)/tau
    d/dt(transit8) <- (transit7 - transit8)/tau
    d/dt(transit9) <- (transit8 - transit9)/tau
    d/dt(transit10) <- (transit9 - transit10)/tau
    frac1 <- max(1 - abs(ncmt - 1), 0)
    frac2 <- max(1 - abs(ncmt - 2), 0)
    frac3 <- max(1 - abs(ncmt - 3), 0)
    frac4 <- max(1 - abs(ncmt - 4), 0)
    frac5 <- max(1 - abs(ncmt - 5), 0)
    frac6 <- max(1 - abs(ncmt - 6), 0)
    frac7 <- max(1 - abs(ncmt - 7), 0)
    frac8 <- max(1 - abs(ncmt - 8), 0)
    frac9 <- max(1 - abs(ncmt - 9), 0)
    frac10 <- max(1 - abs(ncmt - 10), 0)
    drug_effect <- frac1*transit1 + frac2*transit2 + frac3*transit3 + frac4*transit4 + frac5*transit5 + frac6*transit6 + frac7*transit7 + frac8*transit8 + frac9*transit9 + frac10*transit10
    d/dt(tumor) <- kng*tumor - drug_effect*tumor
    tumor ~ prop(propErr) + add(addErr)
  })
}

oncology_sdm_lobo_2002 <- oncology_sdm_lobo_2002()
# Simulation works

oncology_sdm_lobo_2002$foceiModel
#> $inner
#> rxode2 NA model named rx_5ca927d7894c0489889587be8f130e27 model (✔ ready). 
#> $state: transit1, transit2, transit3, transit4, transit5, transit6, transit7, transit8, transit9, transit10, tumor, rx__sens_transit1_BY_ETA_1___, rx__sens_transit2_BY_ETA_1___, rx__sens_transit3_BY_ETA_1___, rx__sens_transit4_BY_ETA_1___, rx__sens_transit5_BY_ETA_1___, rx__sens_transit6_BY_ETA_1___, rx__sens_transit7_BY_ETA_1___, rx__sens_transit8_BY_ETA_1___, rx__sens_transit9_BY_ETA_1___, rx__sens_transit10_BY_ETA_1___, rx__sens_tumor_BY_ETA_1___
#> $params: THETA[1], THETA[2], THETA[3], THETA[4], THETA[5], THETA[6], THETA[7], ETA[1], cp, tumor0
#> $lhs: rx_pred_, rx__sens_rx_pred__BY_ETA_1___, rx_r_, rx__sens_rx_r__BY_ETA_1___
#> 
#> $predOnly
#> rxode2 NA model named rx_8dd22393399eee19b9a5ff1c0d59d283 model (✔ ready). 
#> $state: transit1, transit2, transit3, transit4, transit5, transit6, transit7, transit8, transit9, transit10, tumor
#> $params: THETA[1], THETA[2], THETA[3], THETA[4], THETA[5], THETA[6], THETA[7], ETA[1], cp, tumor0
#> $lhs: rx_pred_, rx_r_, lkng, ltau, lec50, kmax, ncmt, propErr, addErr, eta.kng, kng, tau, taulast, ec50, edrug, frac1, frac2, frac3, frac4, frac5, frac6, frac7, frac8, frac9, frac10, drug_effect, tad, dosenum
#> 
#> $extra.pars
#> NULL
#> 
#> $outer
#> NULL
#> 
#> $predNoLhs
#> rxode2 NA model named rx_2a6406a67beed57e4a4dc2088c2adf61 model (✔ ready). 
#> $state: transit1, transit2, transit3, transit4, transit5, transit6, transit7, transit8, transit9, transit10, tumor
#> $params: THETA[1], THETA[2], THETA[3], THETA[4], THETA[5], THETA[6], THETA[7], ETA[1], cp, tumor0
#> $lhs: rx_pred_, rx_r_
#> 
#> $theta
#> NULL
#> 
#> $pred.minus.dv
#> [1] TRUE
#> 
#> $log.thetas
#> integer(0)
#> 
#> $log.etas
#> integer(0)
#> 
#> $extraProps
#> list()
#> 
#> $eventTheta
#> [1] 0 0 0 0 0 0 0
#> 
#> $eventEta
#> [1] 0
#> 
#> attr(,"class")
#> [1] "foceiModelList"

^{Created on 2024-07-03 with reprex v2.1.0}

[mattfidler]

library(nlmixr2est)
#> Loading required package: nlmixr2data

oncology_sdm_lobo_2002 <- function() {
  description <- "Signal transduction model for delayed concentration effects on cancer cell growth"
  reference <- "Lobo ED, Balthasar JP. Pharmacodynamic modeling of chemotherapeutic effects: Application of a transit compartment model to characterize methotrexate effects in vitro. AAPS J. 2002;4(4):212-222. doi:10.1208/ps040442"
  # Values for lkng, ltau, lec50, and kmax are for methotrexate from Lobo 2002,
  # Table 2.  propErr and addErr are added as reasonable values though not from
  # Lobo 2002 where no value is apparent in the paper.
  ini({
    lkng <- log(0.02) ; label("Cell net growth rate (growth minus death) (1/hr)")
    ltau <- log(c(1, 34.1, 500)) ; label("Mean transit time of each transit compartment (hr)")
    lec50 <- log(c(1, 50, 2000)) ; label("Drug concentration reducing the cell growth by 50% (ug/mL)")
    kmax <- 0.01 ; label("Maximum drug-related reduction in cell growth (1/hr)")
    ncmt <- c(1, 4, 10) ; label("Number of transit compartments")

    propErr <- c(0, 0.3) ; label("Proportional residual error (fraction)")
    addErr <- c(0, 50, 1000) ; label("Additive residual error (tumor volume units)")
  })
  model({
    # cp is the drug concentration
    kng <- exp(lkng)
    tau <- exp(ltau)
    taulast <- tau
    ec50 <- exp(lec50)

    edrug <- kmax*cp/(ec50 + cp)

    tumor(0) <- tumor0
    d/dt(transit1) <- (edrug - transit1)/tau
    d/dt(transit2) <- (transit1 - transit2)/tau
    d/dt(transit3) <- (transit2 - transit3)/tau
    d/dt(transit4) <- (transit3 - transit4)/tau
    d/dt(transit5) <- (transit4 - transit5)/tau
    d/dt(transit6) <- (transit5 - transit6)/tau
    d/dt(transit7) <- (transit6 - transit7)/tau
    d/dt(transit8) <- (transit7 - transit8)/tau
    d/dt(transit9) <- (transit8 - transit9)/tau
    d/dt(transit10) <- (transit9 - transit10)/tau
    frac1 <- max(1 - abs(ncmt - 1), 0)
    frac2 <- max(1 - abs(ncmt - 2), 0)
    frac3 <- max(1 - abs(ncmt - 3), 0)
    frac4 <- max(1 - abs(ncmt - 4), 0)
    frac5 <- max(1 - abs(ncmt - 5), 0)
    frac6 <- max(1 - abs(ncmt - 6), 0)
    frac7 <- max(1 - abs(ncmt - 7), 0)
    frac8 <- max(1 - abs(ncmt - 8), 0)
    frac9 <- max(1 - abs(ncmt - 9), 0)
    frac10 <- max(1 - abs(ncmt - 10), 0)
    drug_effect <- frac1*transit1 + frac2*transit2 + frac3*transit3 + frac4*transit4 + frac5*transit5 + frac6*transit6 + frac7*transit7 + frac8*transit8 + frac9*transit9 + frac10*transit10
    d/dt(tumor) <- kng*tumor - drug_effect*tumor
    tumor ~ prop(propErr) + add(addErr)
  })
}

oncology_sdm_lobo_2002 <- oncology_sdm_lobo_2002()

d_sim <- data.frame(id=1, time=0:24, tumor0=100, cp=1:25, dv=101:125)
fit <- nlmixr(oncology_sdm_lobo_2002, data = d_sim, est="focei", control=foceiControl(print=0))
#> rxode2 2.1.3.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 1760
#> → compress parHistData in nlmixr2 object, save 181952

print(fit)
#> ── nlmixr² Population Only (outer: nlminb) ──
#> 
#>         OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> Pop 6.783208 66.73013 75.26227      -26.36507        3617.961        7.042943
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance table compress    other
#> elapsed 0.015019 0.022317   0.022318 0.025    0.011 1.166346
#> 
#> ── ($parFixed or $parFixedDf): ──
#> 
#>                                                          Parameter     Est.
#> lkng              Cell net growth rate (growth minus death) (1/hr)    -4.59
#> ltau            Mean transit time of each transit compartment (hr)     3.35
#> lec50   Drug concentration reducing the cell growth by 50% (ug/mL)      3.8
#> kmax          Maximum drug-related reduction in cell growth (1/hr)     5.17
#> ncmt                                Number of transit compartments     4.02
#> propErr                     Proportional residual error (fraction) 0.000792
#> addErr                Additive residual error (tumor volume units)     0.87
#>              SE   %RSE Back-transformed(95%CI) BSV(SD) Shrink(SD)%
#> lkng    0.00307 0.0671 0.0102 (0.0101, 0.0103)                    
#> ltau     0.0354   1.06       28.5 (26.6, 30.6)                    
#> lec50     0.104   2.74         44.9 (36.6, 55)                    
#> kmax      0.161   3.12       5.17 (4.85, 5.48)                    
#> ncmt     0.0373  0.929       4.02 (3.95, 4.09)                    
#> propErr                               0.000792                    
#> addErr                                    0.87                    
#>  
#>   Covariance Type ($covMethod): |r|
#>   Information about run found ($runInfo):
#>    • gradient problems with initial estimate and covariance; see $scaleInfo 
#>    • R matrix non-positive definite but corrected by R = sqrtm(R%*%R) 
#>    • last objective function was not at minimum, possible problems in optimization 
#>   Censoring ($censInformation): No censoring
#>   Minimization message ($message):  
#>     false convergence (8) 
#>   In an ODE system, false convergence may mean "useless" evaluations were performed.
#>   See https://tinyurl.com/yyrrwkce
#>   It could also mean the convergence is poor, check results before accepting fit
#>   You may also try a good derivative free optimization:
#>     nlmixr2(...,control=list(outerOpt="bobyqa"))
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 25 × 33
#>   ID     TIME    DV IPRED  IRES IWRES transit1  transit2    transit3    transit4
#>   <fct> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>     <dbl>       <dbl>       <dbl>
#> 1 1         0   101  100  1      1.15  0       0         0               0      
#> 2 1         1   102  101. 0.974  1.12  0.00388 0.0000685 0.000000835     7.97e-9
#> 3 1         2   103  102. 0.938  1.07  0.0113  0.000335  0.00000726      1.26e-7
#> # ℹ 22 more rows
#> # ℹ 23 more variables: transit5 <dbl>, transit6 <dbl>, transit7 <dbl>,
#> #   transit8 <dbl>, transit9 <dbl>, transit10 <dbl>, tumor <dbl>, kng <dbl>,
#> #   tau <dbl>, taulast <dbl>, ec50 <dbl>, edrug <dbl>, frac1 <dbl>,
#> #   frac2 <dbl>, frac3 <dbl>, frac4 <dbl>, frac5 <dbl>, frac6 <dbl>,
#> #   frac7 <dbl>, frac8 <dbl>, frac9 <dbl>, frac10 <dbl>, drug_effect <dbl>

^{Created on 2024-07-03 with reprex v2.1.0}

[mattfidler]

Needs some tests.

[mattfidler]

This should work, perhaps someone should review the transformations.

The max and min are recursive functions based on https://stackoverflow.com/questions/30738923/max-implemented-with-basic-operators

This gives rxGt and rxLt functions since they are handled in rxode2 dsl/symengine setup.

For abs(a) you would have:

abs(a) = 2*(a)*(a > 0) - a