Problem with variable name `b1`

[tamas-ferenci]

Interestingly, nlmixr doesn't seem to like if I name a variable b1. (But b2 is OK!) As an example:

one.compartment.saem <- function() {
ini({
 tb1 <- .5   # Log b1
 tcl <- -3.2 # Log Cl
 tv <- -1    # Log V
 eta.b1 ~ 1
 eta.cl ~ 2
 eta.v ~ 1
 add.err <- 0.1
})
model({
 b1 <- exp(tb1 + eta.b1)
 cl <- exp(tcl + eta.cl)
 v <- exp(tv + eta.v)
 d/dt(depot) = -b1 * depot
 d/dt(center) = b1 * depot - cl / v * center
 cp = center / v
 cp ~ add(add.err)
})
}

I consulted the help for naming issues ("These parameters can be called by almost any R-compatible name. Please note that:"), but I haven't found anything that'd make b1 an illegal name...

[mattfidler]

To be more specific SAEM does not support b1 (FOCEi does)

library(nlmixr)
one.compartment.saem <- function() {
    ini({
        tb1 <- .5   # Log b1
        tcl <- -3.2 # Log Cl
        tv <- -1    # Log V
        eta.b1 ~ 1
        eta.cl ~ 2
        eta.v ~ 1
        add.err <- 0.1
    })
    model({
        b1 <- exp(tb1 + eta.b1)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv + eta.v)
        d/dt(depot) = -b1 * depot
        d/dt(center) = b1 * depot - cl / v * center
        cp = center / v
        cp ~ add(add.err)
    })
}

fit2 <- nlmixr(one.compartment.saem, theo_sd, est="focei", control=foceiControl(print=500))
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gradient
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |    #| Objective Fun |       tb1 |       tcl |        tv |   add.err |
#> |.....................|           |           |           |...........|
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> Calculating covariance matrix
#>                            done
#> Calculating residuals/tables
#> done.
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Gradient problems with initial estimate and covariance;
#> see $scaleInfo
fit2
#> ── nlmixr FOCEi (outer: nlminb) fit ───────────────────────────────────────
#>           OBJF      AIC      BIC Log-likelihood Condition Number
#> FOCEi 116.0091 372.6089 392.7885      -179.3044         102.2606
#> 
#> ── Time (sec; fit2$time): ─────────────────────────────────────────────────
#>            setup optimize covariance table
#> elapsed 0.558933 0.361135   0.361137 0.013
#> 
#> ── Population Parameters (fit2$parFixed or fit2$parFixedDf): ──────────────
#>         Parameter   Est.     SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tb1        Log b1  0.451   0.19 42.2       1.57 (1.08, 2.28)    72.0%
#> tcl        Log Cl  -3.22 0.0829 2.58 0.0401 (0.0341, 0.0472)    26.8%
#> tv          Log V -0.776 0.0424 5.47       0.46 (0.423, 0.5)    13.8%
#> add.err            0.685                               0.685         
#>         Shrink(SD)%
#> tb1         0.525% 
#> tcl          3.70% 
#> tv           10.6% 
#> add.err             
#> 
#>   Covariance Type (fit2$covMethod): r,s
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance (fit2$omega) or correlation (fit2$omegaR; diagonals=SDs)
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in fit2$shrink
#>   Minimization message (fit2$message): relative convergence (4)
#> 
#> ── Fit Data (object fit2 is a modified tibble): ───────────────────────────
#> # A tibble: 132 x 21
#>   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.08    0     0.74   1.08   0    0.74 
#> 2 1     0.25   2.84  2.80 0.0356 0.0213  3.86 -1.02  -1.49   2.64 0.205
#> 3 1     0.570  6.57  5.02 1.55   0.699   6.78 -0.208 -0.304  4.77 1.80 
#> # … with 129 more rows, and 10 more variables: CWRES <dbl>, eta.b1 <dbl>,
#> #   eta.cl <dbl>, eta.v <dbl>, b1 <dbl>, cl <dbl>, v <dbl>, cp <dbl>,
#> #   depot <dbl>, center <dbl>

fit <- nlmixr(one.compartment.saem, theo_sd, est="saem", control=saemControl(print=500))
#> Compiling RxODE equations...
#> done.
#> PKG_CXXFLAGS=-I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/nlmixr//include -I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders//include -I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/Rcpp//include -I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/RcppArmadillo//include -I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/RcppEigen//include -I/home/matt/R/x86_64-pc-linux-gnu-library/3.5/BH//include
#> PKG_LIBS= $(BLAS_LIBS) $(LAPACK_LIBS)
#> Error in dyn.load(saem.dll): unable to load shared object '/home/matt/.rxCache/saem1db116c0eee2.so':
#>   /home/matt/.rxCache/saem1db116c0eee2.so: cannot open shared object file: No such file or directory
fit
#> Error in eval(expr, envir, enclos): object 'fit' not found

^{Created on 2019-01-19 by the reprex package (v0.2.1)}

[mattfidler]

I have a fix, but I need it reviewed first.

[tamas-ferenci]

Thank you very much in advance! (Totally unrelated sidenote: my R output looks different for the very same input. But I absolutely don't know the reprex package, so it might be an intentional feature; sorry in that case for mentioning this.)

[mattfidler]

Sometimes parts of the calculation are cached, so you may see more output with the same input. Of course, colors and lines look different in different systems. My output looks like the following:

The reprex package makes sure the issue is reproducible and gives output similar to knitr, I believe

[tamas-ferenci]

Thanks for the explanation (re reprex)!

[mattfidler]

This should now be fixed. You may test if you wish.

[tamas-ferenci]

I checked it, works perfectly, thank you very much!