Final estimates violates the boundaries (SAEM)

[gbygbygby]

Hi!

I am building a sequential zero and first order absorption model, using SAEM, here is the code:

run1 <- function(){
  ini({
    tv <- c(0, 130)
    tcl <- c(0, 7)
    tmat <- c(0, 2)
    td1 <- c(0, 0.3, 1)

    eta.v ~ 0.9
    eta.cl ~ 0.9
    eta.mat ~ 0.9

    add_sigma <- 0.9
  })
  model({
    mu1 <- log(tcl)
    mu2 <- log(tv)
    mu3 <- log(tmat)
    
    cl <- exp(mu1 + eta.cl) 
    v <- exp(mu2 + eta.v)
    mat <- exp(mu3 + eta.mat)

    d1 <- mat * (1 - td1)
    ka <- 1 / (mat - d1)
    
    d/dt(depot) = - ka * depot
    d/dt(centr) = ka * depot - (cl/v) * centr
    dur(depot) = d1

    depot(0) = 0
    centr(0) = 0
    
    cp = log(1000 * centr / v + 0.00001)
    
    cp ~ add(add_sigma)
  })
}
model1 <- nlmixr2(run1, dt, "saem")
print(model1)

But it gives me some very strange final estimations:

is there any options I can set to avoid violating the boundaries？

Thanks in advance！

[mattfidler]

Note, the algorithm saem does not honor any boundaries, so you need to use the transformation to allow the model to estimate the values.

Also, as written, none of your expressions are mu-referenced according to nlmixr2 (you can see what is mu-referenced by nlmixr2(model)). To get the best estimates in mu referenced algorithms like saem, you need to maximize the number of mu-referenced expressions.

So my suggestions for this model are:

• Drop all boundaries in the ini() section, and use parameter transformations in the model
• Use the back-transformations to get the estimated values in saem
• I would also suggest using lognormal distribution instead of the log add distribution that you have in your model. The objective functions will be comparable to other error models on the non-transformed scale, and it will simulate on the normal scale as well.

Here is the model that you had (and I use nlmixr2() to examine it). You will see the difference in the mu referencing between the two models when it is printed out.

library(nlmixr2)
#> Loading required package: nlmixr2data

run1 <- function(){
  ini({
    tv <- c(0, 130)
    tcl <- c(0, 7)
    tmat <- c(0, 2)
    td1 <- c(0, 0.3, 1)

    eta.v ~ 0.9
    eta.cl ~ 0.9
    eta.mat ~ 0.9

    add_sigma <- 0.9
  })
  model({
    mu1 <- log(tcl)
    mu2 <- log(tv)
    mu3 <- log(tmat)

    cl <- exp(mu1 + eta.cl)
    v <- exp(mu2 + eta.v)
    mat <- exp(mu3 + eta.mat)

    d1 <- mat * (1 - td1)
    ka <- 1 / (mat - d1)

    d/dt(depot) = - ka * depot
    d/dt(centr) = ka * depot - (cl/v) * centr
    dur(depot) = d1

    depot(0) = 0
    centr(0) = 0

    cp = log(1000 * centr / v + 0.00001)

    cp ~ add(add_sigma)
  })
}

nlmixr2(run1)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>        tv       tcl      tmat       td1 add_sigma 
#>     130.0       7.0       2.0       0.3       0.9 
#> 
#> Omega ($omega): 
#>         eta.v eta.cl eta.mat
#> eta.v     0.9    0.0     0.0
#> eta.cl    0.0    0.9     0.0
#> eta.mat   0.0    0.0     0.9
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2            centr
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         tv <- c(0, 130)
#>         tcl <- c(0, 7)
#>         tmat <- c(0, 2)
#>         td1 <- c(0, 0.3, 1)
#>         add_sigma <- c(0, 0.9)
#>         eta.v ~ 0.9
#>         eta.cl ~ 0.9
#>         eta.mat ~ 0.9
#>     })
#>     model({
#>         mu1 <- log(tcl)
#>         mu2 <- log(tv)
#>         mu3 <- log(tmat)
#>         cl <- exp(mu1 + eta.cl)
#>         v <- exp(mu2 + eta.v)
#>         mat <- exp(mu3 + eta.mat)
#>         d1 <- mat * (1 - td1)
#>         ka <- 1/(mat - d1)
#>         d/dt(depot) = -ka * depot
#>         d/dt(centr) = ka * depot - (cl/v) * centr
#>         dur(depot) = d1
#>         depot(0) = 0
#>         centr(0) = 0
#>         cp = log(1000 * centr/v + 1e-05)
#>         cp ~ add(add_sigma)
#>     })
#> }

run2 <- function(){
  ini({
    tv <- log(130)
    tcl <- log(7) # Lower boundary is implied by exp()
    tmat <- log(2)
    td1 <- logit(0.3)

    eta.v ~ 0.9
    eta.cl ~ 0.9
    eta.mat ~ 0.9

    lnorm.sd <- 0.9
  })
  model({
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    mat <- exp(tmat + eta.mat)

    d1 <- mat * (1 - expit(td1))
    ka <- 1 / (mat - d1)

    d/dt(depot) = - ka * depot
    d/dt(centr) = ka * depot - (cl/v) * centr
    dur(depot) = d1

    depot(0) = 0
    centr(0) = 0

    cp = 1000 * centr / v + 0.00001

    cp ~ lnorm(lnorm.sd)
  })
}

nlmixr2(run2)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>         tv        tcl       tmat        td1   lnorm.sd 
#>  4.8675345  1.9459101  0.6931472 -0.8472979  0.9000000 
#> 
#> Omega ($omega): 
#>         eta.v eta.cl eta.mat
#> eta.v     0.9    0.0     0.0
#> eta.cl    0.0    0.9     0.0
#> eta.mat   0.0    0.0     0.9
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2            centr
#>  ── μ-referencing ($muRefTable): ──  
#>   theta     eta level
#> 1   tcl  eta.cl    id
#> 2    tv   eta.v    id
#> 3  tmat eta.mat    id
#> 
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         tv <- 4.86753445045558
#>         tcl <- 1.94591014905531
#>         tmat <- 0.693147180559945
#>         td1 <- -0.847297860387204
#>         lnorm.sd <- c(0, 0.9)
#>         eta.v ~ 0.9
#>         eta.cl ~ 0.9
#>         eta.mat ~ 0.9
#>     })
#>     model({
#>         cl <- exp(tcl + eta.cl)
#>         v <- exp(tv + eta.v)
#>         mat <- exp(tmat + eta.mat)
#>         d1 <- mat * (1 - expit(td1))
#>         ka <- 1/(mat - d1)
#>         d/dt(depot) = -ka * depot
#>         d/dt(centr) = ka * depot - (cl/v) * centr
#>         dur(depot) = d1
#>         depot(0) = 0
#>         centr(0) = 0
#>         cp = 1000 * centr/v + 1e-05
#>         cp ~ lnorm(lnorm.sd)
#>     })
#> }

^{Created on 2024-12-02 with reprex v2.1.1}