symengine-based RxODE release

mattfidler commented 4 years ago

symengine just went to CRAN, do this to say goodbye to python dependencies.

[x] Double check authorship.
[x] ~~Create backward-compatible simulation using mvnfast~~
[x] Deprecate old functions and arguments with .Deprecated().
[x] Update & rebuild pkgdown files.
[x] Check RxODE menu links.

[ ] Check spelling in the package with spelling::spell_check_package(".").

devtools::load_all();devtools::document();spelling::spell_check_package("~/src/RxODE")
devtools::load_all();spelling::spell_check_package("~/src/RxODE")

[x] Check with `rch/r-debug:

RDvalgrind -d "valgrind --tool=memcheck --leak-check=yes --show-reachable=yes --num-callers=20" -e "library(RxODE);rxTest(F)"

[x] ~~Check with RDsan~~ Currently RcppEigen mvnfast and ggforce won't compile
[x] Check with RDcsan

[x] Stress test function pointers/garbage collection

devtools::install("/home/src/RxODE");library(RxODE);sapply(1:100, function(...){rxTest("scale")})

[x] Check nlmixr works with this version of RxODE

mattfidler commented 4 years ago

See https://github.com/symengine/symengine.R/issues/60

Will arrive at https://cran.r-project.org/package=symengine eventually

mattfidler commented 4 years ago

Yea. It has arrived in CRAN

mattfidler commented 4 years ago

It has some issues currently....

rikardn commented 4 years ago

I could help testing this. Even if it doesn't fully work is it possible to test?

mattfidler commented 4 years ago

That is generous. If you can install symengine, then you can test if you wish; The instructions are below:

Install the appropriate C/C++ toolchain to compile symengine/RxODE/nlmixr from source:

First setup R with the instructions for use with the c++ stan tool-chain

Install `symengine` and `units`

Ubuntu Linux packages: In addition to some other (forgotten) packages you will need to have:

libudunits2-dev
cmake
libgmp-dev
libmpfr-dev
libmpc-dev
devscripts

Mac OSX packages using homebrew; This may not be necessary after CRAN provides mac binaries

cmake
gmp
mpfr
libmpc

Windows setup This simply works because the symenigne folks attached binaries to their build; All you need to do is install the cran package.

Installing development `RxODE` and `nlmixr`

Once complete with the prerequisites you can install using devtools:

## To use the latest github symengine you may wish to install from github
devtools::install_github("symengine/symengine.R") # 
devtools::install_github("nlmixrdevelopment/RxODE") # prunes Branching ie `if/else`
devtools::install_github("nlmixrdevelopment/nlmixr") # adds focei etas on dur and other things

What is incomplete and needs to be fixed before release

My current todo before releasing this is:

[x] Get CRAN requested updates to RxODE and nlmixr completed
[x] Fully test linCmt() models in nlmixr
[x] Update & test nlmixr graphics, including adding a stat for censoring support using ggproto
[x] Finish up the levels in RxODE (which is 50-90% complete)
[x] Finish up tests and test changes to the dynmodel code from our summer intern
[x] Fix any bugs that came up in testing

Then everything should be ready.

I have been stalled on this testing because of CRAN requests to update RxODE/nlmixr and other packages.

rikardn commented 4 years ago

Thanks for helping me getting started!

Installation

The first issue I encountered was the dependency of R version 3.6 for symengine. Since I had 3.5.2 I needed to upgrade. As 3.6 is only about 11 months old this dependency could be a problem for some users.

I successfully installed the development versions of RxODE and nlmixr as per the instructions. This was very painless.

Testing the pheno example

Now I tested the pheno example. Here is the final part of the output:

Calculating covariance matrix
Error in rxSolve_(object, .ctl, .nms, .xtra, params, events, inits, setupOnly = .setupOnly) : 
  Not initialized
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
[====|====|====|====|====|====|====|====|====|====] 0:00:00 

calculate sensitivities
[====|====|====|====|====|====|====|====|====|====] 0:00:00 

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
compiling EBE model...done
Calculating residuals/tables
done.
Compiling NPDE model...done
[====|====|====|====|====|====|====|====|====|====] 0:00:00 

Warning messages:
1: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE,  :
  Covariance matrix non-positive definite, corrected by sqrtm(fim %*% fim)
2: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE,  :
  Linearization of FIM could not be used to calculate covariance.
3: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE,  :
  NAs introduced by coercion

The warnings and the error in rxSolve_ wasn't shown using the CRAN version of nlmixr and RxODE for the same modelfit.

Also the fit object cannot be printed:

 Error: 'rxIsNotebook' is not an exported object from 'namespace:RxODE' 
5.  stop(gettextf("'%s' is not an exported object from 'namespace:%s'", 
    name, getNamespaceName(ns)), call. = FALSE, domain = NA) 
4.  getExportedValue(pkg, name) 
3.  RxODE::rxIsNotebook 
2.  print.nlmixrFitCore(x) 
1.  (function (x, ...) 
UseMethod("print"))(x)

I could reach components of the fit object though:

> fit$parFixedDf
         Estimate         SE      %RSE Back-transformed    CI Lower    CI Upper BSV(CV%) Shrink(SD)%
tcl     -4.994730 0.06686136  1.338638      0.006773548 0.005941604 0.007721981 52.11585    3.112099
tv       0.342029 0.05315817 15.542007      1.407801130 1.268507361 1.562390635 41.73459    1.462043
add.err  2.770992         NA        NA      2.770992195          NA          NA       NA          NA

Testing nlmixr from python

The original itch I was scratching was that nlmixr could not run through rpy2 in python beacuse of a weird cyclic python calls R calls python situation in the end crashing the python session. This is a known issue with rpy2 and reticulate. I tested again with this new version:

from rpy2.robjects import r, pandas2ri
pandas2ri.activate()

params = r('source("pheno.R", keep.source=TRUE); fit$parFixedDf')
print(params)

With the old nlmixr this crashed, but now it works! params is now a pandas DataFrame. So by removing the python dependency of nlmixr it can now somewhat paradoxically be called from python.

R Session info

My session info in case needed:

> devtools::session_info()
Session info --------------------------------------------------------------------------------------------------------------
 setting  value                       
 version  R version 3.6.2 (2019-12-12)
 system   x86_64, linux-gnu           
 ui       RStudio (1.1.463)           
 language (EN)                        
 collate  en_US.UTF-8                 
 tz       Europe/Stockholm            
 date     2020-03-13                  

Packages ------------------------------------------------------------------------------------------------------------------
 package       * version     date       source                                   
 assertthat      0.2.1       2019-03-21 cran (@0.2.1)                            
 backports       1.1.5       2019-10-02 cran (@1.1.5)                            
 base          * 3.6.2       2019-12-13 local                                    
 brew            1.0-6       2011-04-13 CRAN (R 3.4.4)                           
 checkmate       2.0.0       2020-02-06 cran (@2.0.0)                            
 cli             1.0.0       2017-11-05 CRAN (R 3.4.4)                           
 colorspace      1.4-1       2019-03-18 CRAN (R 3.6.2)                           
 compiler        3.6.2       2019-12-13 local                                    
 crayon          1.3.4       2017-09-16 CRAN (R 3.4.4)                           
 curl            4.3         2019-12-02 CRAN (R 3.6.2)                           
 data.table      1.12.8      2019-12-09 cran (@1.12.8)                           
 datasets      * 3.6.2       2019-12-13 local                                    
 devtools        1.13.6      2018-06-27 CRAN (R 3.4.4)                           
 digest          0.6.25      2020-02-23 CRAN (R 3.6.2)                           
 dparser         0.1.8       2017-11-13 cran (@0.1.8)                            
 dplyr           0.8.5       2020-03-07 cran (@0.8.5)                            
 evaluate        0.14        2019-05-28 cran (@0.14)                             
 generics        0.0.2       2018-11-29 cran (@0.0.2)                            
 ggplot2         3.3.0       2020-03-05 cran (@3.3.0)                            
 glue            1.3.2       2020-03-12 cran (@1.3.2)                            
 graphics      * 3.6.2       2019-12-13 local                                    
 grDevices     * 3.6.2       2019-12-13 local                                    
 grid            3.6.2       2019-12-13 local                                    
 gtable          0.3.0       2019-03-25 cran (@0.3.0)                            
 highr           0.7         2018-06-09 CRAN (R 3.4.4)                           
 httr            1.3.1       2017-08-20 CRAN (R 3.4.4)                           
 huxtable        4.7.1       2020-01-08 cran (@4.7.1)                            
 knitr           1.20        2018-02-20 CRAN (R 3.4.4)                           
 lattice         0.20-38     2018-11-04 CRAN (R 3.6.0)                           
 lazyeval        0.2.2       2019-03-15 CRAN (R 3.6.2)                           
 lbfgs           1.2.1       2014-08-31 cran (@1.2.1)                            
 lifecycle       0.2.0       2020-03-06 cran (@0.2.0)                            
 lotri           0.2.1       2020-02-12 cran (@0.2.1)                            
 magrittr        1.5         2014-11-22 CRAN (R 3.4.4)                           
 memoise         1.1.0       2017-04-21 CRAN (R 3.4.4)                           
 methods       * 3.6.2       2019-12-13 local                                    
 munsell         0.5.0       2018-06-12 CRAN (R 3.4.4)                           
 n1qn1           6.0.1-6     2020-03-02 cran (@6.0.1-6)                          
 nlme            3.1-143     2019-12-10 CRAN (R 3.6.2)                           
 nlmixr        * 1.1.1-5     2020-03-13 Github (nlmixrdevelopment/nlmixr@825aa5c)
 parallel        3.6.2       2019-12-13 local                                    
 pillar          1.4.3       2019-12-20 cran (@1.4.3)                            
 pkgconfig       2.0.3       2019-09-22 cran (@2.0.3)                            
 PreciseSums     0.3         2018-04-12 cran (@0.3)                              
 purrr           0.3.3       2019-10-18 cran (@0.3.3)                            
 R6              2.2.2       2017-06-17 CRAN (R 3.4.4)                           
 Rcpp            1.0.3       2019-11-08 CRAN (R 3.6.2)                           
 RcppArmadillo   0.9.850.1.0 2020-02-09 cran (@0.9.850)                          
 rex             1.1.2       2017-10-19 CRAN (R 3.4.4)                           
 rlang           0.4.5       2020-03-01 cran (@0.4.5)                            
 rstudioapi      0.7         2017-09-07 CRAN (R 3.4.4)                           
 RxODE         * 0.9.2-0     2020-03-13 Github (nlmixrdevelopment/RxODE@6fb8688) 
 scales          1.1.0       2019-11-18 CRAN (R 3.6.2)                           
 stats         * 3.6.2       2019-12-13 local                                    
 symengine       0.1.0       2020-03-06 CRAN (R 3.6.2)                           
 sys             3.3         2019-08-21 cran (@3.3)                              
 tibble          2.1.3       2019-06-06 cran (@2.1.3)                            
 tidyselect      1.0.0       2020-01-27 cran (@1.0.0)                            
 tools           3.6.2       2019-12-13 local                                    
 units           0.6-5       2019-10-08 cran (@0.6-5)                            
 utils         * 3.6.2       2019-12-13 local                                    
 vpc             1.1.0       2018-08-27 cran (@1.1.0)                            
 withr           2.1.2       2018-03-15 CRAN (R 3.4.4)                           
 yaml            2.2.1       2020-02-01 CRAN (R 3.6.2)

mattfidler commented 4 years ago

Thanks for the information @rikardn

The printout problem is a consequence of dropping rstudioapi yesterday as a dependency in RxODE. I did not update nlmixr to drop the rstudioapi yet.

In the new branch there is some differences in handling models which caused the second issue in saem. It seems to be intermittent, but needs to be taken care of before a release.

I am glad that it now works with python by dropping the python dependency though!

mattfidler commented 4 years ago

Printing is fixed

mattfidler commented 4 years ago

If you get a chance, I would be interested if the python will import the whole object as a pandas data frame:

from rpy2.robjects import r, pandas2ri
pandas2ri.activate()

params = r('source("pheno.R", keep.source=TRUE); fit')
print(params)

rikardn commented 4 years ago

So, I reinstalled RxODE and nlmixr and then first run the pheno example in R. Printing the fit now displays another error:

> fit
── nlmixr SAEM(ODE); OBJF by FOCEi approximation fit ────────────────────────────────────────────────────────────────────────

          OBJF      AIC      BIC Log-likelihood Condition Number
FOCEi 688.9656 985.8366 1004.097      -486.9183         18.32301

── Time (sec fit$time): ─────────────────────────────────────────────────────────────────────────────────────────────────────

          saem    setup table cwres covariance  npde    other
elapsed 28.496 1.364507 0.021 1.394      0.004 0.421 0.645493

── Population Parameters (fit$parFixed or fit$parFixedDf): ──────────────────────────────────────────────────────────────────

                         Parameter  Est.     SE %RSE    Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
tcl     typical value of clearance -4.99 0.0669 1.34 0.00677 (0.00594, 0.00772)     52.1      3.11% 
tv         typical value of volume 0.342 0.0532 15.5          1.41 (1.27, 1.56)     41.7      1.46% 
add.err       residual variability  2.77                                   2.77                     

  Covariance Type (fit$covMethod): |fim|
  Correlations in between subject variability (BSV) matrix:
 Hide Traceback

 Rerun with Debug
 Error: $ operator is invalid for atomic vectors 
4..getR(x) 
3..getCorPrint(x$omega) 
2.print.nlmixrFitCore(x) 
1.(function (x, ...) 
UseMethod("print"))(x)

Then I tried from python to see what happens with the fit object. Unfortunately it simply gives this DataFrame:

     ID   TIME    DV  EVID  ...         A1      EPRED       ERES      NPDE
1     1    2.0  17.3     0  ...  24.772881  18.423742  -1.123742 -0.103835
2     1  112.5  31.0     0  ...  40.068479  26.688259   4.311741  0.806433
3     2    2.0   9.7     0  ...  14.865701   9.358943   0.341057 -0.279513
4     2   63.5  24.6     0  ...  27.514981  16.593295   8.006705  1.388891
5     2  135.5  33.0     0  ...  38.903866  25.647451   7.352549  0.119393
..   ..    ...   ...   ...  ...        ...        ...        ...       ...
151  58   47.5  27.9     0  ...  31.580148  19.340219   8.559781  0.763733
152  58  131.8  31.0     0  ...  40.394087  20.710667  10.289333 -1.321555
153  59    1.8  22.6     0  ...  22.620130  13.052445   9.547555  0.779487
154  59   73.8  34.3     0  ...  32.220030  21.302482  12.997518  0.909069
155  59  146.8  40.2     0  ...  39.038945  22.562727  17.637273  0.119085

[155 rows x 23 columns]

Next I tried without the automatic pandas conversion of the R object. Below is the printout. I cannot see any entries to get to the diffrent sub-objects. But it is really no a big problem since it is always possible to do r('fit$anything') so I am happy.

R object with classes: ('nlmixrSaem', 'nlmixrFitData', 'nlmixrFitCore', 'tbl_df', 'tbl', 'data.frame') mapped to:
[IntSexpVe..., FloatSexp..., FloatSexp..., IntSexpVe..., ..., FloatSexp..., FloatSexp..., FloatSexp..., FloatSexp...]
  ID: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
[1, 1, 2, 2, ..., 58, 59, 59, 59]
  TIME: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[2.000000, 112.500000, 2.000000, 63.500000, ..., 131.800000, 1.800000, 73.800000, 146.800000]
  DV: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[17.300000, 31.000000, 9.700000, 24.600000, ..., 31.000000, 22.600000, 34.300000, 40.200000]
  EVID: <class 'rpy2.robjects.vectors.IntVector'>
  R object with classes: ('RTYPES.INTSXP',) mapped to:
[0, 0, 0, 0, ..., 0, 0, 0, 0]
...
  RES: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[24.772881, 40.068479, 14.865701, 27.514981, ..., 40.394087, 22.620130, 32.220030, 39.038945]
  WRES: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[18.423742, 26.688259, 9.358943, 16.593295, ..., 20.710667, 13.052445, 21.302482, 22.562727]
  IPRED: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[-1.123742, 4.311741, 0.341057, 8.006705, ..., 10.289333, 9.547555, 12.997518, 17.637273]
  IRES: <class 'rpy2.robjects.vectors.FloatVector'>
  R object with classes: ('numeric',) mapped to:
[-0.103835, 0.806433, -0.279513, 1.388891, ..., -1.321555, 0.779487, 0.909069, 0.119085]

mattfidler commented 4 years ago

Thanks @rikardn

For the python piece, I agree: accessing them from R is fine; I was more worried that you wouldn't be able to read in the main dataset since it had these extra pieces attached.

For me, I cannot reproduce your printing error:

library(nlmixr)
one.compartment <- function() {
    ini({
        tka <- 0.45 # Log Ka
        tcl <- 1 # Log Cl
        tv <- 3.45    # Log V
        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)
        d/dt(depot) = -ka * depot
        d/dt(center) = ka * depot - cl / v * center
        cp = center / v
        cp ~ add(add.sd)
    })
}
fit2 <- nlmixr(one.compartment, theo_sd, est="saem", control=list(print=0))
#> Compiling RxODE equations...done.
#> Building SAEM model...done
#> Calculating covariance matrix
#> Error in rxSolve_(object, .ctl, .nms, .xtra, params, events, inits, setupOnly = .setupOnly) : 
#>   Not initialized
#> 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
#> Calculating residuals/tables
#> done.
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : Covariance matrix non-positive definite, corrected by sqrtm(fim %*%
#> fim)
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : Linearization of FIM could not be used to calculate covariance.
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : NAs introduced by coercion
print(fit2)
#> ── nlmixr SAEM(ODE); OBJF not calculated fit ───────────────────────────────────
#> 
#> 
#>  Gaussian/Laplacian Likelihoods: AIC() or $objf etc. 
#>  FOCEi CWRES & Likelihoods: addCwres()
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#> 
#>           saem    setup table covariance    other
#> elapsed 14.391 1.245571 0.011      0.004 0.464429
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#> 
#> 
#>        Parameter  Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka       Log Ka  0.45  0.189 41.9       1.57 (1.08, 2.27)     70.3    -0.513% 
#> tcl       Log Cl  1.02 0.0833 8.19       2.76 (2.35, 3.25)     28.0      4.23% 
#> tv         Log V  3.45 0.0396 1.15         31.5 (29.1, 34)     12.6      9.93% 
#> add.sd           0.698                               0.698                     
#>  
#>   Covariance Type ($covMethod): |fim|
#>   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
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> 
#> # A tibble: 132 x 18
#>   ID     TIME    DV  EVID  PRED    RES IPRED   IRES  IWRES eta.ka eta.cl   eta.v
#>   <fct> <dbl> <dbl> <int> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
#> 1 1     0      0.74     0  0     0.74   0     0.74   1.06  0.0982 -0.480 -0.0819
#> 2 1     0.25   2.84     0  3.26 -0.418  3.84 -1.00  -1.44  0.0982 -0.480 -0.0819
#> 3 1     0.570  6.57     0  5.84  0.734  6.78 -0.213 -0.305 0.0982 -0.480 -0.0819
#> # … with 129 more rows, and 6 more variables: ka <dbl>, cl <dbl>, v <dbl>,
#> #   cp <dbl>, depot <dbl>, center <dbl>

^{Created on 2020-03-13 by the reprex package (v0.3.0)}

Session info

``` r devtools::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 3.6.0 (2019-04-26) #> os Gentoo/Linux #> system x86_64, linux-gnu #> ui X11 #> language (EN) #> collate en_US.utf8 #> ctype en_US.utf8 #> tz US/Central #> date 2020-03-13 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib #> assertthat 0.2.1 2019-03-21 [1] #> backports 1.1.5 2019-10-02 [1] #> brew 1.0-6 2011-04-13 [1] #> callr 3.4.2 2020-02-12 [1] #> checkmate 2.0.0 2020-02-06 [1] #> cli 2.0.2 2020-02-28 [1] #> colorspace 1.4-1 2019-03-18 [1] #> crayon 1.3.4 2017-09-16 [1] #> data.table 1.12.8 2019-12-09 [1] #> desc 1.2.0 2018-05-01 [1] #> devtools 2.2.2 2020-02-17 [1] #> digest 0.6.25 2020-02-23 [1] #> dparser 0.1.8 2017-11-13 [1] #> dplyr 0.8.5 2020-03-07 [1] #> ellipsis 0.3.0 2019-09-20 [1] #> evaluate 0.14 2019-05-28 [1] #> fansi 0.4.1 2020-01-08 [1] #> fs 1.3.2 2020-03-05 [1] #> generics 0.0.2 2018-11-29 [1] #> ggplot2 3.3.0 2020-03-05 [1] #> glue 1.3.2 2020-03-12 [1] #> gtable 0.3.0 2019-03-25 [1] #> highr 0.8 2019-03-20 [1] #> htmltools 0.4.0 2019-10-04 [1] #> huxtable 4.7.1 2020-01-08 [1] #> knitr 1.28 2020-02-06 [1] #> lattice 0.20-38 2018-11-04 [2] #> lazyeval 0.2.2 2019-03-15 [1] #> lifecycle 0.2.0 2020-03-06 [1] #> lotri 0.2.2 2020-03-13 [1] #> magrittr 1.5 2014-11-22 [1] #> Matrix 1.2-17 2019-03-22 [2] #> memoise 1.1.0 2017-04-21 [1] #> munsell 0.5.0 2018-06-12 [1] #> n1qn1 6.0.1-6 2020-03-02 [1] #> nlme 3.1-139 2019-04-09 [2] #> nlmixr * 1.1.1-5 2020-03-13 [1] #> pillar 1.4.3 2019-12-20 [1] #> pkgbuild 1.0.6 2019-10-09 [1] #> pkgconfig 2.0.3 2019-09-22 [1] #> pkgload 1.0.2 2018-10-29 [1] #> PreciseSums 0.3 2018-04-12 [1] #> prettyunits 1.1.1 2020-01-24 [1] #> processx 3.4.2 2020-02-09 [1] #> ps 1.3.2 2020-02-13 [1] #> purrr 0.3.3 2019-10-18 [1] #> R6 2.4.1 2019-11-12 [1] #> Rcpp 1.0.3 2019-11-08 [1] #> RcppArmadillo 0.9.850.1.0 2020-02-09 [1] #> remotes 2.1.1 2020-02-15 [1] #> rex 1.1.2 2017-10-19 [1] #> rlang 0.4.5 2020-03-01 [1] #> rmarkdown 2.1 2020-01-20 [1] #> rprojroot 1.3-2 2018-01-03 [1] #> RxODE * 0.9.2-0 2020-03-13 [1] #> scales 1.1.0 2019-11-18 [1] #> sessioninfo 1.1.1 2018-11-05 [1] #> stringi 1.4.6 2020-02-17 [1] #> stringr 1.4.0 2019-02-10 [1] #> symengine 0.1.0 2020-03-06 [1] #> sys 3.3 2019-08-21 [1] #> testthat 2.3.2 2020-03-02 [1] #> tibble 2.1.3 2019-06-06 [1] #> tidyselect 1.0.0 2020-01-27 [1] #> units 0.6-5 2019-10-08 [1] #> usethis 1.5.1 2019-07-04 [1] #> utf8 1.1.4 2018-05-24 [1] #> vctrs 0.2.4 2020-03-10 [1] #> vpc 1.1.0 2018-08-27 [1] #> withr 2.1.2 2018-03-15 [1] #> xfun 0.12 2020-01-13 [1] #> yaml 2.2.1 2020-02-01 [1] #> source #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> local #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> Github (nlmixrdevelopment/nlmixr@3396737) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> Github (nlmixrdevelopment/RxODE@0bd2cc1) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> CRAN (R 3.6.0) #> #> [1] /home/matt/R/x86_64-pc-linux-gnu-library/3.6 #> [2] /usr/lib64/R/library ```

Can you provide more information so I can fix the printing of the output?

rikardn commented 4 years ago

I ran this:

library(nlmixr)

pheno <- function() {
    ini({
        tcl <- log(0.008) # typical value of clearance
        tv <-  log(0.6)   # typical value of volume
        ## var(eta.cl)
        eta.cl + eta.v ~ c(1,
                           0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
        # interindividual variability on clearance and volume
        add.err <- 0.1    # residual variability
    })
    model({
        cl <- exp(tcl + eta.cl) # individual value of clearance
        v <- exp(tv + eta.v)    # individual value of volume
        ke <- cl / v            # elimination rate constant
        d/dt(A1) = - ke * A1    # model differential equation
        cp = A1 / v             # concentration in plasma
        cp ~ add(add.err)       # define error model
    })
}

fit <- nlmixr(pheno, pheno_sd, "saem", table=list(cwres=TRUE, npde=TRUE))

session info

``` Session info setting value version R version 3.6.2 (2019-12-12) system x86_64, linux-gnu ui RStudio (1.1.463) language (EN) collate en_US.UTF-8 tz Europe/Stockholm date 2020-03-13 Packages ------------- package * version assertthat 0.2.1 backports 1.1.5 base * 3.6.2 brew 1.0-6 checkmate 2.0.0 cli 2.0.2 colorspace 1.4-1 compiler 3.6.2 crayon 1.3.4 data.table 1.12.8 datasets * 3.6.2 devtools 1.13.6 digest 0.6.25 dparser 0.1.8 dplyr 0.8.5 fansi 0.4.1 generics 0.0.2 ggplot2 3.3.0 glue 1.3.2 graphics * 3.6.2 grDevices * 3.6.2 grid 3.6.2 gtable 0.3.0 huxtable 4.7.1 knitr 1.28 lattice lazyeval 0.2.2 lifecycle 0.2.0 lotri 0.2.1 magrittr 1.5 memoise 1.1.0 methods * 3.6.2 munsell 0.5.0 n1qn1 nlme nlmixr * 1.1.1-5 parallel 3.6.2 pillar 1.4.3 pkgconfig 2.0.3 PreciseSums 0.3 purrr 0.3.3 R6 2.4.1 Rcpp 1.0.3 RcppArmadillo rex 1.1.2 rlang 0.4.5 rstudioapi 0.11 RxODE * 0.9.2-0 scales 1.1.0 stats * 3.6.2 symengine 0.1.0 sys 3.3 tibble 2.1.3 tidyselect 1.0.0 tools 3.6.2 units 0.6-5 utils * 3.6.2 vpc 1.1.0 withr 2.1.2 xfun 0.12 yaml 2.2.1 ``` --------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------ date source 2019-03-21 cran (@0.2.1) 2019-10-02 cran (@1.1.5) 2019-12-13 local 2011-04-13 CRAN (R 3.4.4) 2020-02-06 cran (@2.0.0) 2020-02-28 cran (@2.0.2) 2019-03-18 CRAN (R 3.6.2) 2019-12-13 local 2017-09-16 CRAN (R 3.4.4) 2019-12-09 cran (@1.12.8) 2019-12-13 local 2018-06-27 CRAN (R 3.4.4) 2020-02-23 CRAN (R 3.6.2) 2017-11-13 cran (@0.1.8) 2020-03-07 cran (@0.8.5) 2020-01-08 cran (@0.4.1) 2018-11-29 cran (@0.0.2) 2020-03-05 cran (@3.3.0) 2020-03-12 cran (@1.3.2) 2019-12-13 local 2019-12-13 local 2019-12-13 local 2019-03-25 cran (@0.3.0) 2020-01-08 cran (@4.7.1) 2020-02-06 cran (@1.28) 0.20-38 2018-11-04 CRAN (R 3.6.0) 2019-03-15 CRAN (R 3.6.2) 2020-03-06 cran (@0.2.0) 2020-02-12 cran (@0.2.1) 2014-11-22 CRAN (R 3.4.4) 2017-04-21 CRAN (R 3.4.4) 2019-12-13 local 2018-06-12 CRAN (R 3.4.4) 6.0.1-6 2020-03-02 cran (@6.0.1-6) 3.1-143 2019-12-10 CRAN (R 3.6.2) 2020-03-13 Github (nlmixrdevelopment/nlmixr@3396737) 2019-12-13 local 2019-12-20 cran (@1.4.3) 2019-09-22 cran (@2.0.3) 2018-04-12 cran (@0.3) 2019-10-18 cran (@0.3.3) 2019-11-12 cran (@2.4.1) 2019-11-08 CRAN (R 3.6.2) 0.9.850.1.0 2020-02-09 cran (@0.9.850) 2017-10-19 CRAN (R 3.4.4) 2020-03-01 cran (@0.4.5) 2020-02-07 cran (@0.11) 2020-03-13 Github (nlmixrdevelopment/RxODE@0bd2cc1) 2019-11-18 CRAN (R 3.6.2) 2019-12-13 local 2020-03-06 CRAN (R 3.6.2) 2019-08-21 cran (@3.3) 2019-06-06 cran (@2.1.3) 2020-01-27 cran (@1.0.0) 2019-12-13 local 2019-10-08 cran (@0.6-5) 2019-12-13 local 2018-08-27 cran (@1.1.0) 2018-03-15 CRAN (R 3.4.4) 2020-01-13 cran (@0.12) 2020-02-01 CRAN (R 3.6.2)

mattfidler commented 4 years ago

Thanks @rikardn

Updated the nlmixr repository with the printing fixes.

rikardn commented 4 years ago

Great! I think this conclues my initial round of testing. I will continue using these versions of nlmixr/RxODE and open new issues if I find any problems.

mattfidler commented 4 years ago

Thanks @rikardn

mattfidler commented 4 years ago

Last fix will deal with the saem problem:

#Error in rxSolve_(object, .ctl, .nms, .xtra, params, events, inits, setupOnly = .setupOnly) :
#  Not initialized

It will also seems to deal with the random crashes when restarting saem.

billdenney commented 4 years ago

I was able to follow the installation instructions, and the installation was successful (with the addition of adding Makevars.win as described in #209).

I then tested the same pheno example that @rikardn used, and I got a different printing error:

library(nlmixr)

pheno <- function() {
  ini({
    tcl <- log(0.008) # typical value of clearance
    tv <-  log(0.6)   # typical value of volume
    ## var(eta.cl)
    eta.cl + eta.v ~ c(1,
                       0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
    # interindividual variability on clearance and volume
    add.err <- 0.1    # residual variability
  })
  model({
    cl <- exp(tcl + eta.cl) # individual value of clearance
    v <- exp(tv + eta.v)    # individual value of volume
    ke <- cl / v            # elimination rate constant
    d/dt(A1) = - ke * A1    # model differential equation
    cp = A1 / v             # concentration in plasma
    cp ~ add(add.err)       # define error model
  })
}

fit <- nlmixr(pheno, pheno_sd, "saem", table=list(cwres=TRUE, npde=TRUE))
#> Compiling RxODE equations...done.
#> 1:    -5.0270    0.3073    0.9500    0.9500   49.4073
# [Bill Denney removed 498 lines of minimization]
#> 500:   -4.9947   0.3420   0.2403   0.1606   7.6784
#> Calculating covariance matrix
#> 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 sensitivities
#> calculate d(f)/d(eta)
#> calculate d(R^2)/d(eta)
#> finding duplicate expressions in inner model...
#> optimizing duplicate expressions in inner model...
#> finding duplicate expressions in EBE model...
#> optimizing duplicate expressions in EBE model...
#> compiling inner model...done
#> compiling EBE model...done
#> Calculating residuals/tables
#> done.
#> Compiling NPDE model...done
fit
#> -- nlmixr SAEM(ODE); OBJF by FOCEi approximation fit ---------------------------------------------------------------------------------------------------------
#> 
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition Number
#> FOCEi 688.9656 985.8366 1004.097      -486.9183         6.240286
#> -- Time (sec fit$time): --------------------------------------------------------------------------------------------------------------------------------------
#> 
#> 
#>          saem setup table cwres covariance npde other
#> elapsed 11.75 2.199  0.04  2.24       0.03 0.45 0.601
#> -- Population Parameters (fit$parFixed or fit$parFixedDf): ---------------------------------------------------------------------------------------------------
#> 
#> 
#>                          Parameter  Est.     SE %RSE    Back-transformed(95%CI)
#> tcl     typical value of clearance -4.99 0.0749  1.5 0.00677 (0.00585, 0.00784)
#> tv         typical value of volume 0.342 0.0546   16          1.41 (1.26, 1.57)
#> add.err       residual variability  2.77                                   2.77
#>         BSV(CV%) Shrink(SD)%
#> tcl         52.1      3.11% 
#> tv          41.7      1.46% 
#> add.err                     
#>  
#>   Covariance Type (fit$covMethod): linFim
#>   Correlations in between subject variability (BSV) matrix:
#> Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

^{Created on 2020-06-09 by the reprex package (v0.3.0)}

mattfidler commented 4 years ago

Are you using, R 4.0? This looks like a R version issue.

On Tue, Jun 9, 2020, 8:39 AM Bill Denney notifications@github.com wrote:

I was able to follow the installation instructions, and the installation was successful (with the addition of adding Makevars.win as described in

209 https://github.com/nlmixrdevelopment/RxODE/issues/209).

I then tested the same pheno example that @rikardn https://github.com/rikardn used, and I got a different printing error:

library(nlmixr) pheno <- function() { ini({ tcl <- log(0.008) # typical value of clearance tv <- log(0.6) # typical value of volume

var(eta.cl)
eta.cl + eta.v ~ c(1,
                   0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
# interindividual variability on clearance and volume
add.err <- 0.1    # residual variability
}) model({ cl <- exp(tcl + eta.cl) # individual value of clearance v <- exp(tv + eta.v) # individual value of volume ke <- cl / v # elimination rate constant d/dt(A1) = - ke * A1 # model differential equation cp = A1 / v # concentration in plasma cp ~ add(add.err) # define error model }) } fit <- nlmixr(pheno, pheno_sd, "saem", table=list(cwres=TRUE, npde=TRUE))#> Compiling RxODE equations...done.#> 1: -5.0270 0.3073 0.9500 0.9500 49.4073# [Bill Denney removed 498 lines of minimization]#> 500: -4.9947 0.3420 0.2403 0.1606 7.6784#> Calculating covariance matrix#> 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 sensitivities#> calculate d(f)/d(eta)#> calculate d(R^2)/d(eta)#> finding duplicate expressions in inner model...#> optimizing duplicate expressions in inner model...#> finding duplicate expressions in EBE model...#> optimizing duplicate expressions in EBE model...#> compiling inner model...done#> compiling EBE model...done#> Calculating residuals/tables#> done.#> Compiling NPDE model...donefit#> -- nlmixr SAEM(ODE); OBJF by FOCEi approximation fit ---------------------------------------------------------------------------------------------------------#> #> #> OBJF AIC BIC Log-likelihood Condition Number#> FOCEi 688.9656 985.8366 1004.097 -486.9183 6.240286#> -- Time (sec fit$time): --------------------------------------------------------------------------------------------------------------------------------------#> #> #> saem setup table cwres covariance npde other#> elapsed 11.75 2.199 0.04 2.24 0.03 0.45 0.601#> -- Population Parameters (fit$parFixed or fit$parFixedDf): ---------------------------------------------------------------------------------------------------#> #> #> Parameter Est. SE %RSE Back-transformed(95%CI)#> tcl typical value of clearance -4.99 0.0749 1.5 0.00677 (0.00585, 0.00784)#> tv typical value of volume 0.342 0.0546 16 1.41 (1.26, 1.57)#> add.err residual variability 2.77 2.77#> BSV(CV%) Shrink(SD)%#> tcl 52.1 3.11% #> tv 41.7 1.46% #> add.err #> #> Covariance Type (fit$covMethod): linFim#> Correlations in between subject variability (BSV) matrix:#> Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

Created on 2020-06-09 by the reprex package https://reprex.tidyverse.org (v0.3.0)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/RxODE/issues/170#issuecomment-641303657, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWX32TLDMNG6NCY5FVTRVY3Q7ANCNFSM4LBRP56A .

billdenney commented 4 years ago

R 3.6.3

billdenney commented 4 years ago

I just updated everything. I now have R 4.0.1 and updated versions of libraries. I'm still getting the same error. The problem appears to be in print.nlmixrFitCore(). x$cor is NULL here:

https://github.com/nlmixrdevelopment/nlmixr/blob/17b688417334e6bc750efb7ac45a976922d2b1c0/R/print.R#L265

Edit: And, here is the traceback:

> traceback()
8: array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), 
       NULL) else NULL)
7: as.matrix.default(x)
6: as.matrix(x)
5: lower.tri(.rs)
4: .getR(x)
3: .getCorPrint(x$cor)
2: print.nlmixrFitCore(x)
1: (function (x, ...) 
   UseMethod("print"))(x)

mattfidler commented 4 years ago

Thanks. I will look at it ona Friday.

On Wed, Jun 10, 2020, 10:08 AM Bill Denney notifications@github.com wrote:

I just updated everything. I now have R 4.0.1 and updated versions of libraries. I'm still getting the same error. The problem appears to be in print.nlmixrFitCore(). x$cor is NULL here:

https://github.com/nlmixrdevelopment/nlmixr/blob/17b688417334e6bc750efb7ac45a976922d2b1c0/R/print.R#L265

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/RxODE/issues/170#issuecomment-642072235, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWWZ4NDDMAFBNCFB5UDRV6OVTANCNFSM4LBRP56A .

mattfidler commented 4 years ago

Out of curiosity, does tidy(fit) work?

On Wed, Jun 10, 2020, 11:44 AM Matthew Fidler matthew.fidler@gmail.com wrote:

Thanks. I will look at it ona Friday.

On Wed, Jun 10, 2020, 10:08 AM Bill Denney notifications@github.com wrote:

I just updated everything. I now have R 4.0.1 and updated versions of libraries. I'm still getting the same error. The problem appears to be in print.nlmixrFitCore(). x$cor is NULL here:

https://github.com/nlmixrdevelopment/nlmixr/blob/17b688417334e6bc750efb7ac45a976922d2b1c0/R/print.R#L265

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/RxODE/issues/170#issuecomment-642072235, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWWZ4NDDMAFBNCFB5UDRV6OVTANCNFSM4LBRP56A .

billdenney commented 4 years ago

tidy(fit) works (with a deprecation warning):

tidy(fit)
# A tibble: 6 x 7
  effect   group         term              estimate std.error statistic  p.value
  <chr>    <chr>         <chr>                <dbl>     <dbl>     <dbl>    <dbl>
1 fixed    NA            tcl                 -4.99     0.0749    -66.7   1.00e+0
2 fixed    NA            tv                   0.342    0.0546      6.26  1.94e-9
3 ran_pars ID            sd__eta.cl           0.490   NA          NA    NA      
4 ran_pars ID            sd__eta.v            0.401   NA          NA    NA      
5 ran_pars ID            cor__eta.v.eta.cl    0.948   NA          NA    NA      
6 ran_pars Residual(add) add.err              2.77    NA          NA    NA      
Warning message:
`as.tbl()` is deprecated as of dplyr 1.0.0.
Please use `tibble::as_tibble()` instead.

mattfidler commented 4 years ago

Ok. Fixed the issue.

mattfidler commented 3 years ago

For now I will close this issue. It seems a bit old.

nlmixrdevelopment / RxODE