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 when piping omega blog #276

Closed emmavestesson closed 4 years ago

emmavestesson commented 4 years ago

Hi,

I get an error when I try to add an omega block to a model through piping. I've included an example below (I've removed some of the printed iterations to make it shorter). I'm using version 1.1.1-5 of nlmixr, version 3.6.1 of R and python 2.7. Please let me know if you need more information. 

library(tidyverse)
library(nlmixr)
library(broom)

one_cmt <- function(){
  #  Note we will log transform parameters
  ini({
    tvcl <- log(0.008) 
    tv <-  log(0.6)   
    eta.cl ~ 0.1      
    eta.v ~ 0.1       
    add.err <- 1      
  })
  model({
    cl <- exp(tvcl + eta.cl)
    v <- exp(tv +eta.v )   
    ke <- cl / v            
    d/dt(A1) = - ke * A1    
    cp = A1 / v             
    cp ~ add(add.err)       
  })
}

one_cmt_mod_add_err <- nlmixr(one_cmt, data=theo_sd, est='saem')
#> Compiling RxODE equations...done.
#> 1:   -4.8107e+00   2.7433e+00   9.5000e-02   9.5000e-02   2.4872e+02
#> 2:   -4.8179e+00   3.0754e+00   1.0546e-01   9.0250e-02   1.1300e+02
#> 3:    -4.8815    3.3915    0.1002    0.0857   45.8380
#> 4:    -4.8501    3.6551    0.0995    0.0815   22.0038
#> 5:    -4.8989    3.9173    0.0945    0.0774   11.2863
[Removed some printing to make example shorter]
#> 496:   -8.8320e+00   4.1428e+00   6.0776e-01   2.1962e-05   8.2007e+00
#> 497:   -8.8317e+00   4.1428e+00   6.0665e-01   2.1966e-05   8.2007e+00
#> 498:   -8.8326e+00   4.1428e+00   6.0638e-01   2.1961e-05   8.2006e+00
#> 499:   -8.8333e+00   4.1428e+00   6.0659e-01   2.1956e-05   8.2006e+00
#> 500:   -8.8335e+00   4.1428e+00   6.0599e-01   2.1987e-05   8.2007e+00
#> Calculating covariance matrix
#> Calculating residuals/tables
#> done.

one_cmt_mod_add_err_omega <- one_cmt_mod_add_err %>%
  ini(  eta.cl + eta.v ~ c(0.1, 
                           0.01, 0.1)) %>%
  nlmixr(data=theo_sd, est="saem")
#> Error in ini(., eta.cl + eta.v ~ c(0.1, 0.01, 0.1)): Arguments must be named

Created on 2020-02-15 by the reprex package (v0.3.0)

mattfidler commented 4 years ago

Hi @emmavestesson,

Thank you for reporting this; I'm unsure if I have implemented the eta.cl+eta.v ~ syntax yet. I will look into it soon.

mattfidler commented 4 years ago

I checked @emmavestesson and I did not implement adding correlations yet. I have started and will let you know when you can try it out.

mattfidler commented 4 years ago

Hi @emmavestesson ,

Thank you for reporting the issue; The covariance piping now works.

See the reprex with the github versions of RxODE and nlmixr below:

``` r 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=saemControl(print=0)) #> Compiling RxODE equations...done. #> Building SAEM model...done #> Calculating covariance matrix #> Load into sympy...Using sympy via SnakeCharmR #> done #> Freeing Python/SymPy memory...done #> ################################################################################ #> Optimizing expressions in Predictions/EBE model...done #> Compiling Predictions/EBE model...done #> Standardized prediction/ebe models produced. #> Calculating residuals/tables #> done. 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 23.3 11.203 0.13 0.02 1.297 #> #> -- Population Parameters ($parFixed or $parFixedDf): ----------------------- #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)% #> tka Log Ka 0.451 0.196 43.5 1.57 (1.07, 2.31) 71.9 0.411% #> tcl Log Cl 1.02 0.0836 8.22 2.77 (2.35, 3.26) 27.0 3.36% #> tv Log V 3.45 0.0469 1.36 31.5 (28.7, 34.5) 14.0 10.0% #> add.sd 0.692 0.692 #> #> #> Covariance Type ($covMethod): linFim #> 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 #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.07 0.105 -0.487 -0.0800 #> 2 1 0.25 2.84 0 3.26 -0.423 3.86 -1.02 -1.48 0.105 -0.487 -0.0800 #> 3 1 0.570 6.57 0 5.84 0.725 6.81 -0.235 -0.340 0.105 -0.487 -0.0800 #> # ... with 129 more rows, and 6 more variables: ka , cl , v , #> # cp , depot , center fit2 %>% update(eta.cl + eta.v ~ c(0.1, 0.01, 0.1)) %>% nlmixr(est="saem", control=saemControl(print=0)) -> fit3 #> Calculating covariance matrix #> Calculating residuals/tables #> done. print(fit3) #> -- 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 6.99 0.612 0.02 0.01 0.248 #> #> -- Population Parameters ($parFixed or $parFixedDf): ----------------------- #> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)% #> tka Log Ka 0.454 0.195 42.9 1.57 (1.07, 2.31) 71.1 -1.25% #> tcl Log Cl 1.03 0.073 7.07 2.81 (2.44, 3.24) 22.8 -1.45% #> tv Log V 3.45 0.0471 1.37 31.5 (28.8, 34.6) 14.3 -1.37% #> add.sd 0.697 0.697 #> #> #> Covariance Type ($covMethod): linFim #> Correlations in between subject variability (BSV) matrix: #> cor:eta.v,eta.cl #> 0.0318 #> #> #> 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 #> #> 1 1 0 0.74 0 0 0.74 0 0.74 1.06 -0.0773 -0.303 -0.188 #> 2 1 0.25 2.84 0 3.26 -0.422 3.70 -0.860 -1.23 -0.0773 -0.303 -0.188 #> 3 1 0.570 6.57 0 5.84 0.732 6.73 -0.165 -0.237 -0.0773 -0.303 -0.188 #> # ... with 129 more rows, and 6 more variables: ka , cl , v , #> # cp , depot , center devtools::session_info() #> - Session info --------------------------------------------------------------- #> setting value #> version R version 3.6.1 (2019-07-05) #> os Windows 10 x64 #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate English_United States.1252 #> ctype English_United States.1252 #> tz America/Chicago #> date 2020-02-19 #> #> - 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.3.2 2019-09-22 [1] #> cli 2.0.1 2020-01-08 [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.1 2019-09-24 [1] #> digest 0.6.24 2020-02-12 [1] #> dparser 0.1.8 2017-11-13 [1] #> dplyr 0.8.4 2020-01-31 [1] #> ellipsis 0.3.0 2019-09-20 [1] #> evaluate 0.14 2019-05-28 [1] #> fansi 0.4.1 2020-01-08 [1] #> farver 2.0.3 2020-01-16 [1] #> fs 1.3.1 2019-05-06 [1] #> generics 0.0.2 2018-11-29 [1] #> ggforce 0.3.1 2019-08-20 [1] #> ggplot2 3.2.1 2019-08-10 [1] #> glue 1.3.1 2019-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] #> jsonlite 1.6.1 2020-02-02 [1] #> knitr 1.28 2020-02-06 [1] #> lattice 0.20-38 2018-11-04 [1] #> lazyeval 0.2.2 2019-03-15 [1] #> lbfgs 1.2.1 2014-08-31 [1] #> lifecycle 0.1.0 2019-08-01 [1] #> lotri 0.2.1 2020-02-12 [1] #> magrittr 1.5 2014-11-22 [1] #> MASS 7.3-51.4 2019-03-31 [1] #> Matrix 1.2-17 2019-03-22 [1] #> memoise 1.1.0 2017-04-21 [1] #> munsell 0.5.0 2018-06-12 [1] #> mvnfast 0.2.5 2018-01-31 [1] #> n1qn1 6.0.1-3 2018-09-17 [1] #> nlme 3.1-140 2019-05-12 [1] #> nlmixr * 1.1.1-5 2020-02-19 [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] #> polyclip 1.10-0 2019-03-14 [1] #> PreciseSums 0.3 2018-04-12 [1] #> prettyunits 1.1.0 2020-01-09 [1] #> processx 3.4.1 2019-07-18 [1] #> ps 1.3.0 2018-12-21 [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.4 2020-01-28 [1] #> rmarkdown 2.1 2020-01-20 [1] #> rprojroot 1.3-2 2018-01-03 [1] #> RxODE * 0.9.1-9 2020-02-19 [1] #> scales 1.1.0 2019-11-18 [1] #> sessioninfo 1.1.1 2018-11-05 [1] #> SnakeCharmR 1.0.7 2019-11-14 [1] #> stringi 1.4.6 2020-02-17 [1] #> stringr 1.4.0 2019-02-10 [1] #> sys 3.3 2019-08-21 [1] #> testthat 2.3.0 2019-11-05 [1] #> tibble 2.1.3 2019-06-06 [1] #> tidyselect 1.0.0 2020-01-27 [1] #> tweenr 1.0.1 2018-12-14 [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.2 2020-01-24 [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.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.2) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> local #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.0) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> local #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> Github (nlmixrdevelopment/SnakeCharmR@b88ef0c) #> CRAN (R 3.6.2) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> CRAN (R 3.6.1) #> #> [1] c:/R/nlmixr_1.1.1-3/R/library ``` Created on 2020-02-19 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

Note I use update instead of ini, though both work. I also suppress the iteration printing with saemControl(print=0)

mattfidler commented 4 years ago

Let me know if it works for you; If it does not, feel free to re-open an issue. Thank you again for your report.