nlmixrdevelopment / nlmixr

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

Feature request: More specific parameter near boundary error message #544

Closed billdenney closed 2 years ago

billdenney commented 2 years ago

I just ran a model where a parameter ended near its boundary. This is an obvious model problem, but a minor update that I think would help is to give the name of the parameter in the warning message to help more rapidly identify the issue.

The specific error message is "parameter estimate near boundary; covariance not calculated". It appears to come from here:

https://github.com/nlmixrdevelopment/nlmixr/blob/e91c61a127eeb565e3841c93e1134ba3b879e71e/src/inner.cpp#L5308-L5311

(I had hoped to be able to make a PR for this, but since it's C++ code, I'm not the best person to code it.)

mattfidler commented 2 years ago

Hi @billdenney

Its been awhile, thank you for your patience in this issue.

This now works. Take the example below:

library(nlmixr)

one.cmt <- function() {
  ini({
    ## You may label each parameter with a comment
    tka <- 0.45 # Log Ka
    tcl <- c(1.1, 1.2, 3) # Log Cl
    ## This works with interactive models
    ## You may also label the preceding line with label("label text")
    tv <- 3.45; label("log V")
    ## the label("Label name") works with all models
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.sd)
  })
}

fit <- nlmixr(one.cmt, theo_sd, list(print=0), est="focei")
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> → compiling events FD model...
#> ✔ done
#> RxODE 1.1.1 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, : initial ETAs were nudged; (can control by foceiControl(etaNudge=.,
#> etaNudge2=))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : ETAs were reset to zero during optimization; (Can control by
#> foceiControl(resetEtaP=.))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : parameter estimate near boundary; covariance not calculated:
#>    "tcl" 
#>  use 'getVarCov' to calculate anyway
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate; see $scaleInfo
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : dur() specified in model, but no modeled duration (rate=-2) in the
#> dataset
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : rate() specified in model, but no modeled rates (rate=-1) in the
#> dataset

print(fit)
#> ── nlmixr FOCEi (outer: nlminb) fit ────────────────────────────────────────────
#> 
#>           OBJF      AIC     BIC Log-likelihood
#> FOCEi 117.8486 374.4484 394.628      -180.2242
#> 
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#>            setup optimize covariance table    other
#> elapsed 3.336028 0.000101   0.000102 0.013 1.018769
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#> 
#>        Parameter  Est. Back-transformed BSV(CV%) Shrink(SD)%
#> tka       Log Ka 0.452             1.57     70.7      1.94% 
#> tcl       Log Cl   1.1                3     28.5      6.98% 
#> tv         log V  3.45             31.6     13.6      11.1% 
#> add.sd           0.693            0.693                     
#>  
#>   Covariance Type ($covMethod): Boundary issue; Get SEs with `getVarCov()`: "tcl" 
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Minimization message ($message):  
#>     relative convergence (4) 
#> 
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 132 × 20
#>   ID     TIME    DV  PRED    RES   WRES IPRED   IRES  IWRES    CPRED   CRES
#>   <fct> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>    <dbl>  <dbl>
#> 1 1      0     0.74  0     0.74   1.07   0     0.74   1.07  1.62e-16  0.74 
#> 2 1      0.25  2.84  3.25 -0.406 -0.217  3.84 -1.00  -1.45  3.20e+ 0 -0.361
#> 3 1      0.57  6.57  5.80  0.765  0.307  6.78 -0.213 -0.308 5.75e+ 0  0.817
#> # … with 129 more rows, and 9 more variables: CWRES <dbl>, eta.ka <dbl>,
#> #   eta.cl <dbl>, eta.v <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> #   dosenum <dbl>

Created on 2021-08-28 by the reprex package (v2.0.0)

Here the warning states

parameter estimate near boundary; covariance not calculated:
    "tcl" 
 use 'getVarCov' to calculate anyway

And in the output the following is listed:

 Covariance Type ($covMethod): Boundary issue; Get SEs with `getVarCov()`: "tcl" 
mattfidler commented 2 years ago

Mistyped commit in log 49f2091e

billdenney commented 2 years ago

Thank you!