EVID doesn't seem to work - Githubissues

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

EVID doesn't seem to work #121

Closed tamas-ferenci closed 5 years ago

tamas-ferenci commented 5 years ago

Take the basic example with theo_sd from the vignette:

fit <- nlmixr(one.compartment.saem, theo_sd, est="saem")
fit$depot
fit$center

Everything is fine, we can see that dosing goes to the first compartment (i.e. depot), just as specified in the EVIDs of theo_sd. However:

theo_sd2 <- theo_sd
theo_sd2[ theo_sd2$EVID==101, ]$EVID <- 102
fit <- nlmixr(one.compartment.saem, theo_sd2, est="saem")
fit$depot
fit$center

I.e. I have changed that dosing should go to the second compartment (center), yet it still goes to the first, as if nothing has changed! Despite the fact that EVIDs are now set to 102.

mattfidler commented 5 years ago

EVID should be set to 201 for dosing to the second compartment.

library(nlmixr)
one.compartment.saem <- function() {
    ini({
        tka <- .5   # Log Ka
        tcl <- -3.2 # Log Cl
        tv <- -1    # Log V
        eta.ka ~ 1
        eta.cl ~ 2
        eta.v ~ 1
        add.err <- 0.1
    })
    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.err)
    })
}

theo_sd2 <- theo_sd
theo_sd2[ theo_sd2$EVID==101, ]$EVID <- 201
fit <- nlmixr(one.compartment.saem, theo_sd2, 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)
#> done.
#> 1:     0.2373   -3.5659   -0.1803    0.9500    1.9000    0.9500   12.7088
#> 500:   -11.6487   -4.0256   -0.1778    1.1082    0.0113    0.0009    7.8409
#> Calculating residuals/tables
#> done.
fit$depot
#>   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [106] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
fit$center
#>   [1] 4.020000 3.998764 3.971743 3.925722 3.851565 3.707432 3.608258
#>   [8] 3.463659 3.318578 3.109613 2.398831 4.400000 4.374534 4.351083
#>  [15] 4.306405 4.222052 4.081031 3.949843 3.782754 3.625918 3.399457
#>  [22] 2.609606 4.530000 4.503924 4.474167 4.432265 4.338490 4.192564
#>  [29] 4.063722 3.894408 3.737000 3.493616 2.701893 4.400000 4.367084
#>  [36] 4.343721 4.300133 4.203432 4.081668 3.950733 3.784750 3.625817
#>  [43] 3.402737 2.592897 5.860000 5.822851 5.795755 5.737069 5.614330
#>  [50] 5.440901 5.268400 5.049646 4.831862 4.543781 3.497229 4.000000
#>  [57] 3.977068 3.950899 3.903225 3.830755 3.707156 3.595984 3.446023
#>  [64] 3.286932 3.091423 2.407144 4.950000 4.923778 4.897693 4.843872
#>  [71] 4.742030 4.597180 4.451110 4.267704 4.088453 3.831940 2.958924
#>  [78] 4.530000 4.505805 4.479817 4.435883 4.338136 4.200043 4.065503
#>  [85] 3.886633 3.730047 3.495646 2.702152 3.100000 3.080194 3.058550
#>  [92] 3.031222 2.969038 2.874770 2.784708 2.659655 2.568642 2.419494
#>  [99] 1.839435 5.500000 5.456927 5.410736 5.382063 5.265532 5.100328
#> [106] 4.940337 4.731730 4.506068 4.253028 3.323961 4.920000 4.893676
#> [113] 4.867490 4.817601 4.715304 4.554178 4.417512 4.230979 4.053262
#> [120] 3.793226 2.934624 5.300000 5.271787 5.243722 5.188034 5.078428
#> [127] 4.916252 4.756240 4.557384 4.370644 4.097752 3.164906

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

mattfidler commented 5 years ago

The event specification is located at https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-events.html

I hope to use NONMEM -style events in the future.

mattfidler commented 5 years ago

From the Page 2018 course, the easiest graphic to describe this is below:

mattfidler commented 5 years ago

If you need any more clarification let me know.

tamas-ferenci commented 5 years ago

Aaaaaaargh... thank you very much! (To write something constructive: it'd be perhaps nice if - in case it receives a syntactically invalid compartment number through EVID - the nlmixr would give at least a warning, instead of silently redirecting the dosing.)

mattfidler commented 5 years ago

Unfortunately, 102 is also legal. However, in this example if you dose to compartment 3 there is a warning generated:

library(nlmixr)
one.compartment.saem <- function() {
    ini({
        tka <- .5   # Log Ka
        tcl <- -3.2 # Log Cl
        tv <- -1    # Log V
        eta.ka ~ 1
        eta.cl ~ 2
        eta.v ~ 1
        add.err <- 0.1
    })
    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.err)
    })
}

theo_sd2 <- theo_sd
theo_sd2[ theo_sd2$EVID==101, ]$EVID <- 301
fit <- nlmixr(one.compartment.saem, theo_sd2, 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)
#> done.
#> 1:     0.2238   -3.4586   -1.0580    0.9500    1.9850    0.9500   32.7653
#> 500:    -6.7212    4.6264   -4.5298    8.5717    9.7720    4.7478   32.7653
#> Calculating residuals/tables
#> done.
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=12)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=11)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=10)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=9)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=8)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=7)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=6)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=5)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=4)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=3)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=2)
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Dose to Compartment 3 ignored (not in ODE; id=1)
fit$depot
#>   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [106] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
fit$center
#>   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [106] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

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

mattfidler commented 5 years ago

If you are trying to convert nonmem datasets, you can supply them and nlmixr attempts to convert them through the nmDataConvert function

mattfidler commented 5 years ago

Obligatory example using NONMEM-style data:

library(nlmixr)
one.compartment.saem <- function() {
    ini({
        tka <- .5   # Log Ka
        tcl <- -3.2 # Log Cl
        tv <- -1    # Log V
        eta.ka ~ 1
        eta.cl ~ 2
        eta.v ~ 1
        add.err <- 0.1
    })
    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.err)
    })
}

theo_sd2 <- theo_sd
theo_sd2[ theo_sd2$EVID==101, ]$CMT <- 1
theo_sd2[ theo_sd2$EVID==101, ]$EVID <- 1
fit <- nlmixr(one.compartment.saem, theo_sd2, 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)
#> done.
#> 1:    0.1503  -3.3314  -0.6447   0.9500   1.9000   0.9500   9.6167
#> 500:    0.4546  -3.2153  -0.7848   0.4413   0.0713   0.0180   0.4774
#> Calculating residuals/tables
#> done.
fit
#> ── nlmixr SAEM(ODE); OBJF not calculated fit ──────────────────────────────
#>       OBJF AIC BIC Log-likelihood Condition Number
#> SAEMg   NA  NA  NA             NA         21.72835
#> 
#> ── Time (sec; x$time): ────────────────────────────────────────────────────
#>    saem    setup optimize covariance table
#>  13.428 0.567771    2e-06      3e-06 0.004
#> 
#> ── Population Parameters (x$parFixed): ────────────────────────────────────
#>         Parameter   Est.     SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tka        Log Ka  0.455  0.198 43.4       1.58 (1.07, 2.32)    74.5%
#> tcl        Log Cl  -3.22 0.0814 2.53 0.0401 (0.0342, 0.0471)    27.2%
#> tv          Log V -0.785 0.0426 5.43     0.456 (0.42, 0.496)    13.5%
#> add.err            0.691                               0.691         
#>         Shrink(SD)%
#> tka        -0.520% 
#> tcl          3.76% 
#> tv           11.0% 
#> add.err             
#> 
#>   Covariance Type (x$covMethod): fim
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance (x$omega) or correlation (x$omegaR; diagonals=SDs)
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in x$shrink
#> 
#> ── Fit Data (object x is a modified tibble): ──────────────────────────────
#> # A tibble: 132 x 17
#>   ID     TIME    DV  PRED     RES IPRED   IRES  IWRES eta.ka eta.cl  eta.v
#>   <fct> <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
#> 1 1     0      0.74  0    0.74     0     0.74   1.07   0.117 -0.625 -0.212
#> 2 1     0.25   2.84  2.84 0.00441  3.86 -1.02  -1.48   0.117 -0.625 -0.212
#> 3 1     0.570  6.57  5.07 1.50     6.79 -0.219 -0.316  0.117 -0.625 -0.212
#> # … with 129 more rows, and 6 more variables: ka <dbl>, cl <dbl>, v <dbl>,
#> #   cp <dbl>, depot <dbl>, center <dbl>

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

mattfidler commented 5 years ago

You can use the CRAN version of nlmixr with the attached fit.Rdata (inside fit.zip so I can attach it to the issue) so you don't have to setup the python pieces if you wish. (It requires reticulate on mac/linux and a custom SnakeCharmR on windows systems; The python requires sympy)

tamas-ferenci commented 5 years ago

Thank you very much for the so detailed answer!