metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
129 stars 36 forks source link

[F] Reinstate $ENV to provide R objects when parsing the model #115

Closed kylebaron closed 7 years ago

kylebaron commented 8 years ago

For example:

    library(mrgsolve)

    code <- '
    $ENV
    par_conv <- 24

    $PARAM
    CL = 1.2*par_conv, VC = 200, KA = 1.89*par_conv
    '

    mod <- mcode("ENV", code, quiet=TRUE)

    param(mod)

    ## 
    ##  Model parameters (N=3):
    ##  name value . name value
    ##  CL   28.8  | VC   200  
    ##  KA   45.4  | .    .

I'd like to expand this:

vjd commented 8 years ago

Really cool option! I am trying to run through all the scenarios in my head where this can be useful -

  1. omega matrix
  2. sigma matrix

Where else can this be really useful other than keeping the R object call within the model specification?

kylebaron commented 8 years ago

Have this working now:

 code <- '
    $ENV
    mat <- cmat(1,0.8,2)

    $OMEGA object="mat"

    $SIGMA object="mat"

    '

    mod <- mcode("ENV2", code, quiet=TRUE)

    revar(mod)

    ## $omega
    ## $...
    ##         [,1]     [,2]
    ## 1:  1.000000 1.131371
    ## 2:  1.131371 2.000000
    ## 
    ## 
    ## $sigma
    ## $...
    ##         [,1]     [,2]
    ## 1:  1.000000 1.131371
    ## 2:  1.131371 2.000000

    exists("mat")

    ## [1] FALSE
kylebaron commented 8 years ago

@vjd Yes ... we've been doing some PBPK modeling with mrgsolve lately and some huge OMEGA matrices with off diagonals etc .... it's really doesn't work well with the NONMEM-like specification. But would be easy if you could just run a couple of lines of R.

This works ...

    code <- '
    $ENV
    foo <- list(Y = 5, Z = 1)
    yes <- TRUE
    conv.factor <- 2.2

    $PARAM >> object="foo"

    $PARAM >> annotated=yes
    Cl : 1.2*conv.factor : Clearance
    Vc : 200/conv.factor : Volume

    '

    mod <- mcode("ENV3", code, quiet=TRUE)

    param(mod)

    ## 
    ##  Model parameters (N=4):
    ##  name value . name value
    ##  Cl   2.64  | Y    5    
    ##  Vc   90.9  | Z    1

    exists("conv.factor")

    ## [1] FALSE
kylebaron commented 8 years ago

Just an aside:

Keep in mind this only operates when the model is being parsed. It just helps you write more complicated but structured model code.

For a minute, I was also thinking of implementing something like this:

$PARAM CL = 1, VC = 2, KE = CL/VC

It's kind-of along the same lines and sometimes people ask for this. We could do it, but deciding not to. It gives the illusion that KE is inherently some function of CL and VC. But it's not. Maybe when the model is parsed, it is. But if you ever change CL or VC, KE stays the same. The way to do it is:

$PARAM CL = 1, VC = 2
$MAIN double KE = CL/VC;

The values in $PARAM are always just fixed numbers. Anything that is a function of parameters needs be derived in $MAIN.