residual error not assigned to variable defined using lincmt

[crgfn]

In code below, variable Cc is not assigned proportional error as expected. However, the output "sim" variable does exhibit residual error.

rm(list=ls())
require(rxode2)

mymdl <- function(){
  ini({
    lcl <- 1.6
    lvc <- log(90)
    lka <- 0.1
    prop.err <- 0.2
    eta.cl ~ 0.1
    eta.vc ~ 0.1
    eta.ka ~ 0.1
  })
  model({
    cl <- exp(lcl + eta.cl)
    vc <- exp(lvc + eta.vc)
    ka <- exp(lka + eta.ka)
    Cc <- linCmt()
    Cc ~ prop(prop.err)
  })
}  

ev  <- et(amt=100, ii=24, addl=3, cmt="depot")
mysim <- rxSolve(mymdl, ev, nSub=1)
plot(mysim)

[mattfidler]

Hi @crgfn

In nlmixr2 style simulations the sim includes the proportional error, Cc is the linCmt() concentrations and ipredSim is the simulated individual predictions.

This is what is shown in your plots, so this is expected behavior.

The underlying rxode2 low-level model can show this a bit more clearly:

require(rxode2)
#> Loading required package: rxode2
#> Warning: package 'rxode2' was built under R version 4.4.2
#> rxode2 3.0.2.9000 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mymdl <- function(){
  ini({
    lcl <- 1.6
    lvc <- log(90)
    lka <- 0.1
    prop.err <- 0.2
    eta.cl ~ 0.1
    eta.vc ~ 0.1
    eta.ka ~ 0.1
  })
  model({
    cl <- exp(lcl + eta.cl)
    vc <- exp(lvc + eta.vc)
    ka <- exp(lka + eta.ka)
    Cc <- linCmt()
    Cc ~ prop(prop.err)
  })
}  

mymdl <- mymdl()

summary(mymdl$simulationModel)
#> rxode2 3.0.2.9000 model named rx_000ac6be307786ed4194123ebb6a5938_x6 model (✔ ready). 
#> DLL: C:\Users\fidlema3\AppData\Local\Temp\1\RTMPK0~1\rxode2\RX_000~1.RXD/rx_000ac6be307786ed4194123ebb6a5938_x64.dll
#> NULL
#> 
#> Calculated Variables:
#> [1] "cl"       "vc"       "ka"       "Cc"       "ipredSim" "sim"     
#> ── rxode2 Model Syntax ──
#> rxode2({
#>     param(lcl, lvc, lka, prop.err, eta.cl, eta.vc, eta.ka)
#>     cl = exp(lcl + eta.cl)
#>     vc = exp(lvc + eta.vc)
#>     ka = exp(lka + eta.ka)
#>     Cc = linCmt()
#>     rx_yj_ ~ 2
#>     rx_lambda_ ~ 1
#>     rx_low_ ~ 0
#>     rx_hi_ ~ 1
#>     rx_pred_f_ ~ Cc
#>     rx_pred_ ~ rx_pred_f_
#>     rx_r_ ~ (rx_pred_f_ * prop.err)^2
#>     ipredSim = rxTBSi(rx_pred_, rx_lambda_, rx_yj_, rx_low_, 
#>         rx_hi_)
#>     sim = rxTBSi(rx_pred_ + sqrt(rx_r_) * rxerr.Cc, rx_lambda_, 
#>         rx_yj_, rx_low_, rx_hi_)
#>     cmt(Cc)
#>     dvid(1)
#> })

^{Created on 2025-02-20 with reprex v2.1.1}

You can see that the sim has the proportional error.

[crgfn]

Thanks, that's very helpful. Two more clarifying questions:
1., Under what circumstances would Cc and ipredSim be different?
2. How can sim be accessed in the model code block?

[mattfidler]

1. Under what circumstances would Cc and ipredSim be different?

This would occur in situations where Cc may be not the only predicted variable.

rm(list=ls())
require(rxode2)
#> Loading required package: rxode2
#> rxode2 3.0.4 using 10 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

pk.turnover.emax3 <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err) | pca
  })
}

ev  <- et(amt=100, ii=24, addl=3, cmt=1) %>%
  et(seq(0, 100, length.out=60), cmt=5) %>%
  et(seq(0, 100, length.out=60), cmt=6) %>%
  et(id=1:3)

mysim <- rxSolve(pk.turnover.emax3, ev, nSub=1)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
plot(mysim, cp, effect, ipredSim, sim)

^{Created on 2025-02-21 with reprex v2.1.1}

[mattfidler]

Also ipredSim will always be the prediction (whereas Cc may not always be the prediction).

[mattfidler]

How can sim be accessed in the model code block?

You can force the model to an old style rxode2 model and then access it that way:

require(rxode2)
#> Loading required package: rxode2
#> rxode2 3.0.4 using 10 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mymdl <- function(){
  ini({
    lcl <- 1.6
    lvc <- log(90)
    lka <- 0.1
    prop.err <- 0.2
    eta.cl ~ 0.1
    eta.vc ~ 0.1
    eta.ka ~ 0.1
  })
  model({
    cl <- exp(lcl + eta.cl)
    vc <- exp(lvc + eta.vc)
    ka <- exp(lka + eta.ka)
    Cc <- linCmt()
    Cc ~ prop(prop.err)
  })
}

mymdl <- mymdl()

simModel <- mymdl$simulationModel

# Convert it to a string
simModel <- rxode2::rxNorm(simModel)

simModel <- rxode2(simModel)
simModel <- as.function(simModel)

print(simModel)
#> function () 
#> {
#>     model({
#>         param(lcl, lvc, lka, prop.err, eta.cl, eta.vc, eta.ka)
#>         cl = exp(lcl + eta.cl)
#>         vc = exp(lvc + eta.vc)
#>         ka = exp(lka + eta.ka)
#>         Cc = linCmtA(rx__PTR__, t, 0, 1, 1, cl, vc, 0, 0, 0, 
#>             0, 0, 1, 0, 0, ka, 0, 1, 0, 0)
#>         rx_yj_ ~ 2
#>         rx_lambda_ ~ 1
#>         rx_low_ ~ 0
#>         rx_hi_ ~ 1
#>         rx_pred_f_ ~ Cc
#>         rx_pred_ ~ rx_pred_f_
#>         rx_r_ ~ (rx_pred_f_ * prop.err)^2
#>         ipredSim = rxTBSi(rx_pred_, rx_lambda_, rx_yj_, rx_low_, 
#>             rx_hi_)
#>         sim = rxTBSi(rx_pred_ + sqrt(rx_r_) * rxerr.Cc, rx_lambda_, 
#>             rx_yj_, rx_low_, rx_hi_)
#>         cmt(Cc)
#>         dvid(1)
#>     })
#> }
#> <environment: 0x569029dcb600>

# Now you can modify this with model piping
simModel %>%
  model({sim=2})
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#>  ── rxode2-based Pred model ───────────────────────────────────────────────────── 
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     model({
#>         param(lcl, lvc, lka, prop.err, eta.cl, eta.vc, eta.ka)
#>         cl = exp(lcl + eta.cl)
#>         vc = exp(lvc + eta.vc)
#>         ka = exp(lka + eta.ka)
#>         Cc = linCmtA(rx__PTR__, t, 0, 1, 1, cl, vc, 0, 0, 0, 
#>             0, 0, 1, 0, 0, ka, 0, 1, 0, 0)
#>         rx_yj_ ~ 2
#>         rx_lambda_ ~ 1
#>         rx_low_ ~ 0
#>         rx_hi_ ~ 1
#>         rx_pred_f_ ~ Cc
#>         rx_pred_ ~ rx_pred_f_
#>         rx_r_ ~ (rx_pred_f_ * prop.err)^2
#>         ipredSim = rxTBSi(rx_pred_, rx_lambda_, rx_yj_, rx_low_, 
#>             rx_hi_)
#>         sim <- 2
#>         cmt(Cc)
#>         dvid(1)
#>     })
#> }

^{Created on 2025-02-21 with reprex v2.1.1}