search

nlmixrdevelopment / nlmixr

nlmixr: an R package for population PKPD modeling
https://nlmixrdevelopment.github.io/nlmixr/
GNU General Public License v2.0
115 stars 45 forks source link

Error: linCmt are not supported in RxODE #543

Open billdenney opened 3 years ago

billdenney commented 3 years ago

Based on the discussion in #540, I thought the following was allowed. Can you please help me understand how to use linCmt on the right hand side?

library(nlmixr)

nlmixr_onecmt_fixed <- function() {
  ini({
    tcl <- log(c(0, .1, 10)) ; label("Clearance (L/kg/hr)")
    tv <- log(0.1); label("log V (L/kg)")
    add_sd <- 0.7
  })
  model({
    cl <- exp(tcl)
    v <- exp(tv)
    # unit conversion
    cp <- linCmt()*1000/v
    cp ~ add(add_sd)
  })
}

fit <- nlmixr(nlmixr_onecmt_fixed, data=theo_sd, control=list(print=0), est="focei")
#> > creating full model...
#> > pruning branches (`if`/`else`)...
#> v done
#> > loading into symengine environment...
#> Error: function 'linCmt' or its derivatives are not supported in RxODE

Created on 2021-07-12 by the reprex package (v2.0.0)

mattfidler commented 3 years ago

What version of RxODE/nlmixr are you using?

By the way linCmt() returns concentrations so a corrected model would be:

nlmixr_onecmt_fixed <- function() {
  ini({
    tcl <- log(c(0, .1, 10)) ; label("Clearance (L/kg/hr)")
    tv <- log(0.1); label("log V (L/kg)")
    add_sd <- 0.7
  })
  model({
    cl <- exp(tcl)
    v <- exp(tv)/1000
    # unit conversion
    cp <- linCmt()
    cp ~ add(add_sd)
  })
}
mattfidler commented 3 years ago

This may be a limitation of the parsing in nlmixr's UI. Currently is splits things into parameter block and ODE block. It has some limitations, especially with linCmt().

I am hoping to get some time to work on the UI again in RxODE for a more robust, parsing-based UI

mattfidler commented 3 years ago

I guess that #540 is true in nlmixr. The following still work;

library(nlmixr)

nlmixr_onecmt_fixed <- function() {
  ini({
    tcl <- log(c(0.001, .1, 10)) ; label("Clearance (L/kg/hr)")
    eta.cl <- 0.1
    tv <- log(3); label("log V (L/kg)")
    add_sd <- 0.1
  })
  model({
    cl <- exp(tcl + eta.cl)
    v <- exp(tv) / 1000
    # unit conversion
    linCmt() ~ add(add_sd)
  })
}

fit <- nlmixr(nlmixr_onecmt_fixed, data=theo_sd, control=list(print=0), est="focei")
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : parameter estimate near boundary; covariance not calculated
#>  use 'getVarCov' to calculate anyway
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate; see $scaleInfo

print(fit)
#> ── nlmixr Population Only (outer: nlminb) fit ──────────────────────────────────
#> 
#>         OBJF     AIC      BIC Log-likelihood
#> Pop 409.7362 660.336 671.8672       -326.168
#> 
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#>            setup optimize covariance table    other
#> elapsed 1.259577 0.000174   0.000175 0.008 0.337074
#> 
#> ──  ($parFixed or $parFixedDf): ────────────────────────────────────────────────
#> 
#>         Est. Back-transformed
#> tcl      2.3              2.3
#> eta.cl -18.6            -18.6
#> tv      11.1         6.41e+04
#> add_sd  2.88             2.88
#>  
#>   Covariance Type ($covMethod): Boundary issue; Get SEs with getVarCov
#>   Minimization message ($message):  
#>     relative convergence (4) 
#> 
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 x 10
#>   ID     TIME    DV IPRED  IRES  IWRES           cl     v   tad dosenum
#>   <fct> <dbl> <dbl> <dbl> <dbl>  <dbl>        <dbl> <dbl> <dbl>   <dbl>
#> 1 1      0     0.74  4.99 -4.25 -1.48  0.0000000866  64.1  0          1
#> 2 1      0.25  2.84  4.99 -2.15 -0.746 0.0000000866  64.1  0.25       1
#> 3 1      0.57  6.57  4.99  1.58  0.550 0.0000000866  64.1  0.57       1
#> # … with 129 more rows

nlmixr_onecmt_fixed <- function() {
  ini({
    tcl <- log(c(0.001, .1, 10)) ; label("Clearance (L/kg/hr)")
    tv <- log(3); label("log V (L/kg)")
    add_sd <- 0.1
  })
  model({
    cl <- exp(tcl)
    v <- exp(tv) / 1000
    # unit conversion
    linCmt() ~ add(add_sd)
  })
}

fit <- nlmixr(nlmixr_onecmt_fixed, data=theo_sd, control=list(print=0), est="focei")
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : bad solve during optimization
#> Warning in (function (uif, data, est = NULL, control = list(), ...,
#> sum.prod = FALSE, : Hessian reset during optimization; (can control by
#> foceiControl(resetHessianAndEta=.))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : parameter estimate near boundary; covariance not calculated
#>  use 'getVarCov' to calculate anyway
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate; see $scaleInfo

print(fit)
#> ── nlmixr Population Only (outer: nlminb) fit ──────────────────────────────────
#> 
#>         OBJF      AIC      BIC Log-likelihood
#> Pop 602.7245 851.3242 859.9726      -422.6621
#> 
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#>            setup optimize covariance table    other
#> elapsed 1.088069  7.6e-05    7.8e-05 0.006 0.471777
#> 
#> ──  ($parFixed or $parFixedDf): ────────────────────────────────────────────────
#> 
#>        Est. Back-transformed
#> tcl    2.25             9.45
#> tv     30.5         1.69e+13
#> add_sd 4.76             4.76
#>  
#>   Covariance Type ($covMethod): Boundary issue; Get SEs with getVarCov
#>   Minimization message ($message):  
#>     both X-convergence and relative convergence (5) 
#> 
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 x 10
#>   ID     TIME    DV        IPRED  IRES IWRES    cl            v   tad dosenum
#>   <fct> <dbl> <dbl>        <dbl> <dbl> <dbl> <dbl>        <dbl> <dbl>   <dbl>
#> 1 1      0     0.74 0.0000000190 0.740 0.155  9.45 16885658654.  0          1
#> 2 1      0.25  2.84 0.0000000190 2.84  0.596  9.45 16885658654.  0.25       1
#> 3 1      0.57  6.57 0.0000000190 6.57  1.38   9.45 16885658654.  0.57       1
#> # … with 129 more rows

Created on 2021-07-12 by the reprex package (v2.0.0)

billdenney commented 3 years ago

I'm using the CRAN versions right now: nlmixr version 2.0.4 and RxODE version 1.0.9.