nlmixr2 / babelmixr2

The goal of babelmixr2 is to interface nlmixr2 with other pharmacometric tools
https://nlmixr2.github.io/babelmixr2/
8 stars 1 forks source link

Using monolixModel and nonmemModel attributes #113

Closed john-harrold closed 2 months ago

john-harrold commented 3 months ago

This is related to this issue:

https://github.com/nlmixr2/babelmixr2/issues/107

Consider this example here.

library(rxode2)
library(babelmixr2)

my_model <- function() {
  ini({
    # Typical Value of System Parameters 
    TV_Vc           = c(.Machine$double.eps, 1.0, Inf)
    TV_Vt           = c(.Machine$double.eps, 1.0, Inf)
    TV_CL           = c(.Machine$double.eps, 1.0, Inf)
    TV_Q            = c(.Machine$double.eps, 1.0, Inf)
    TV_ka           = c(.Machine$double.eps, 1.0, Inf)
    TV_fb           = c(.Machine$double.eps, 1.0, Inf)

    # Error model parameters
    prop_err     =  c(.Machine$double.eps, 0.1, .Machine$double.xmax)
    add_err      =  c(.Machine$double.eps, 0.1, .Machine$double.xmax)

  })
  model({ 
    # System Parameters 
    Vc          = TV_Vc
    Vt          = TV_Vt
    CL          = TV_CL
    Q           = TV_Q
    ka          = TV_ka
    fb          = TV_fb

    # Placeholders for infusion rates.
    # These should be defined in the dataset.
    # time scale: 1 units: (hours) 
    # mass scale: 1 units: (mg/hour) 
    Dinf      = 0.0

    # Defining ODEs
    d/dt(At)    = (-ka*At)
    d/dt(Cc)    = (ka*At*fb/Vc - CL/Vc*Cc -Q*(Cc-Cp)/Vc  + Dinf/Vc)
    d/dt(Cp)    = (+Q*(Cc-Cp)/Vt)

    # Outputs and error models
    Cc_mg_ml    = Cc
    Cc_mg_ml ~ prop(prop_err) + add(add_err)

  })
}

One option I've tried is the following:

mlx1 = my_model()$monolixModel
nm1  = my_model()$nonmemModel

The monolix model works fine but the NONMEM model fails with:

Error: 'mat' must be a matrix

I think this is because I have no IIV in the model, so I add a IIV to one of the parameters:

# If I add IIV to a parameters 
my_model = my_model |> 
  model({Vc = exp(TV_Vc + eta_Vc)})

But I'm unable to access my_model as a function any longer:

mlx1 = my_model()$monolixModel
Error in my_model() : could not find function "my_model"

Now I try to access these as list elements:

nm2  = my_model$nonmemModel
mlx2 = my_model$monolixModel

Which seems to work, however the monolixModel element has a length of 3.

Can you tell me what I'm doing wrong?

mattfidler commented 3 months ago

But I'm unable to access my_model as a function any longer:

That is because that it isn't a function. If you get back the functional form you can do the following:

library(rxode2)
#> rxode2 2.1.3.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(babelmixr2)
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data

my_model <- function() {
  ini({
    # Typical Value of System Parameters
    TV_Vc           = c(.Machine$double.eps, 1.0, Inf)
    TV_Vt           = c(.Machine$double.eps, 1.0, Inf)
    TV_CL           = c(.Machine$double.eps, 1.0, Inf)
    TV_Q            = c(.Machine$double.eps, 1.0, Inf)
    TV_ka           = c(.Machine$double.eps, 1.0, Inf)
    TV_fb           = c(.Machine$double.eps, 1.0, Inf)

    # Error model parameters
    prop_err     =  c(.Machine$double.eps, 0.1, .Machine$double.xmax)
    add_err      =  c(.Machine$double.eps, 0.1, .Machine$double.xmax)

  })
  model({
    # System Parameters
    Vc          = TV_Vc
    Vt          = TV_Vt
    CL          = TV_CL
    Q           = TV_Q
    ka          = TV_ka
    fb          = TV_fb

    # Placeholders for infusion rates.
    # These should be defined in the dataset.
    # time scale: 1 units: (hours)
    # mass scale: 1 units: (mg/hour)
    Dinf      = 0.0

    # Defining ODEs
    d/dt(At)    = (-ka*At)
    d/dt(Cc)    = (ka*At*fb/Vc - CL/Vc*Cc -Q*(Cc-Cp)/Vc  + Dinf/Vc)
    d/dt(Cp)    = (+Q*(Cc-Cp)/Vt)

    # Outputs and error models
    Cc_mg_ml    = Cc
    Cc_mg_ml ~ prop(prop_err) + add(add_err)

  })
}

my_model = my_model |>
    model({Vc = exp(TV_Vc + eta_Vc)})
#> ℹ add between subject variability `eta_Vc` and set estimate to 1

print(is.function(my_model))
#> [1] FALSE

my_model = my_model |>
    as.function()

print(is.function(my_model))
#> [1] TRUE

Created on 2024-08-03 with reprex v2.1.1

mattfidler commented 3 months ago

As far as "what is wrong", again this interface, ie the $nonmemModel and $monolixModel are not meant to be used directly.

I would still would do the following instead to create NONMEM control streams and monolix mlxtran files:

For example, the following code will create the NONMEM control stream without running NONMEM

nlmixr2(model, dummyData, "nonmem", nonmemControl(modelName="exportModel", runCommand=NA))

A similar command would work with Monolix translation

nlmixr2(model, dummyData, "monolix", monolixControl(modelName="exportModel", runCommand=NA))

This should combine any of the non-compliant output you see into a model you can read back.

It will also perform the checks on the model that you ran into that requires etas (for instance) so you will have better experience building an interface and debugging outputs.

If you ever want to couple an actual rxode2 compatible dataset, then this would also be easier.

john-harrold commented 3 months ago

Awesome thank you.

john-harrold commented 3 months ago

Would it make sense to create functions like as_monolix and as_nonmem that can take in a rxode2 or nlmixr2 object and just do this under the hood?

john-harrold commented 3 months ago

I can do it and add it to bablemixr2 if you're cool with it.

mattfidler commented 3 months ago

The problem is that the models change based on the data, so I don't think there really should be these functions. Especially with monolix as the data specification is in their model file.

If you really want to, that is fine, but we don't use snake_case here either. (We would use base R convention for all the as. functions and use dots.

mattfidler commented 3 months ago

Mu referencing will change based on the data you use and there for both the NONMEM and Monolix models will change.

mattfidler commented 3 months ago

If you write them, because of my strong reservations, perhaps either some warning is always produced or something similar.

Like "these models are produced under the assumption of non-time-varying covariates and should not be used when you have time-varying covariates"

john-harrold commented 3 months ago

No worries. I'll make the function and put it into ruminate since I want to use this functionality there as well.

Are there any quick tests that can be done to determine if an object is an rxode2 or nlmixr2 object. Something like is.rxode2() or is.nlmixr2()?

mattfidler commented 3 months ago

rxode2::assertRxUi(model)

Will assert the input is compatible with rxode2 funtion or can be converted into it.

It will return a rxode2 ui object if it is a funtion string etc that cam be converted.

On Sun, Aug 4, 2024, 11:12 AM John Harrold @.***> wrote:

No worries. I'll make the function and put it into ruminate since I want to use this functionality there as well.

Are there any quick tests that can be done to determine if an object is an rxode2 or nlmixr2 object. Something like is.rxode2() or is.nlmixr2()?

— Reply to this email directly, view it on GitHub https://github.com/nlmixr2/babelmixr2/issues/113#issuecomment-2267593521, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWQXPXBGS53MO3SU4PTZPZHH3AVCNFSM6AAAAABL5LAFH2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENRXGU4TGNJSGE . You are receiving this because you commented.Message ID: @.***>

john-harrold commented 3 months ago

The instructions above for NONMEM works but when I try it with Monolix I get the error below. Here is the example I'm using:

library(rxode2)
library(nlmixr2)
library(babelmixr2)

rx_fun = function () 
{
    ini({
        TV_Vc <- c(2.22044604925031e-16, 1)
        TV_Vt <- c(2.22044604925031e-16, 1)
        TV_CL <- c(2.22044604925031e-16, 1)
        TV_Q <- c(2.22044604925031e-16, 1)
        TV_ka <- c(2.22044604925031e-16, 1)
        TV_fb <- c(2.22044604925031e-16, 1)
        prop_err <- c(2.22044604925031e-16, 0.1)
        add_err <- c(2.22044604925031e-16, 0.1)
        eta_Vc ~ 1
    })
    model({
        Vc <- exp(TV_Vc + eta_Vc)*BOB
        #Vc <- TV_Vc
        Vt = TV_Vt
        CL = TV_CL
        Q = TV_Q
        ka = TV_ka
        fb = TV_fb
        Dinf = 0
        d/dt(At) = (-ka * At)
        d/dt(Cc) = (ka * At * fb/Vc - CL/Vc * Cc - Q * (Cc - 
            Cp)/Vc + Dinf/Vc)
        d/dt(Cp) = (+Q * (Cc - Cp)/Vt)
        Cc_mg_ml = Cc
        Cc_mg_ml ~ add(add_err) + prop(prop_err)
    })
}

rx_obj = rxode2(rx_fun)

dataset = et(id=c(1,2), time=c(0,1), evid=0)
dataset[["dv"]] = NA
for(tmp_cov in rx_obj$allCovs){
  dataset[[tmp_cov]] = 1
}

res = nlmixr2(rx_obj, dataset, "monolix", babelmixr2::monolixControl(modelName="exportModel", runCommand=NA))

The error I get is:

Error : invalid value for monolixControl(runCommand=)
Error: invalid value for monolixControl(runCommand=)

I'm using babelmixr2 0.1.12.

mattfidler commented 3 months ago

This should be fixed by #114. A small fix in nlmixr2est will make this run even better.

john-harrold commented 2 months ago

In case anyone is interested I created a function that uses babelmixr2 to convert rxode2 models into NONMEM and Monolix.

https://ruminate.ubiquity.tools/reference/rx2other.html

This vignette shows examples:

https://r.ubiquity.tools/articles/Howto.html#getting-your-model-in-nonmem-format

mattfidler commented 2 months ago

Great John. You can also use this to create PopED scripts :smile:

john-harrold commented 2 months ago

I'm happy to do it if you can just tell me what I need to do with a dataset and rxode2 object :)

john-harrold commented 2 months ago

Hey Matt

Can you tell me what I need to change here to generate the input files for PopEd?

rx_fun = function () 
{
    ini({
        TV_Vc <- c(2.22044604925031e-16, 1)
        TV_Vt <- c(2.22044604925031e-16, 1)
        TV_CL <- c(2.22044604925031e-16, 1)
        TV_Q <- c(2.22044604925031e-16, 1)
        TV_ka <- c(2.22044604925031e-16, 1)
        TV_fb <- c(2.22044604925031e-16, 1)
        prop_err <- c(2.22044604925031e-16, 0.1)
        add_err <- c(2.22044604925031e-16, 0.1)
        eta_Vc ~ 1
    })
    model({
        Vc <- exp(TV_Vc + eta_Vc)*BOB
        #Vc <- TV_Vc
        Vt = TV_Vt
        CL = TV_CL
        Q = TV_Q
        ka = TV_ka
        fb = TV_fb
        Dinf = 0
        d/dt(At) = (-ka * At)
        d/dt(Cc) = (ka * At * fb/Vc - CL/Vc * Cc - Q * (Cc - 
            Cp)/Vc + Dinf/Vc)
        d/dt(Cp) = (+Q * (Cc - Cp)/Vt)
        Cc_mg_ml = Cc
        Cc_mg_ml ~ add(add_err) + prop(prop_err)
    })
}

rx_obj = rxode2(rx_fun)
#> Error in rxode2(rx_fun): could not find function "rxode2"

dataset = et(id=c(1,2), time=c(0,1), evid=0)
#> Error in et(id = c(1, 2), time = c(0, 1), evid = 0): could not find function "et"
dataset[["dv"]] = NA
#> Error: object 'dataset' not found
for(tmp_cov in rx_obj$allCovs){
  dataset[[tmp_cov]] = 1
}
#> Error: object 'rx_obj' not found

res = nlmixr2(rx_obj, dataset, "poped", babelmixr2::popedControl(modelName="exportModel", runCommand=NA))
#> Error in nlmixr2(rx_obj, dataset, "poped", babelmixr2::popedControl(modelName = "exportModel", : could not find function "nlmixr2"

Created on 2024-09-28 with reprex v2.1.1

mattfidler commented 2 months ago

Sorry John,

As of now it simply runs the poped analysis... you definitely need the data.

On Sat, Sep 28, 2024, 11:08 AM John Harrold @.***> wrote:

Hey Matt

Can you tell me what I need to change here to generate the input files for PopEd?

rx_fun = function () { ini({ TV_Vc <- c(2.22044604925031e-16, 1) TV_Vt <- c(2.22044604925031e-16, 1) TV_CL <- c(2.22044604925031e-16, 1) TV_Q <- c(2.22044604925031e-16, 1) TV_ka <- c(2.22044604925031e-16, 1) TV_fb <- c(2.22044604925031e-16, 1) prop_err <- c(2.22044604925031e-16, 0.1) add_err <- c(2.22044604925031e-16, 0.1) eta_Vc ~ 1 }) model({ Vc <- exp(TV_Vc + eta_Vc)*BOB

Vc <- TV_Vc

    Vt = TV_Vt
    CL = TV_CL
    Q = TV_Q
    ka = TV_ka
    fb = TV_fb
    Dinf = 0
    d/dt(At) = (-ka * At)
    d/dt(Cc) = (ka * At * fb/Vc - CL/Vc * Cc - Q * (Cc -
        Cp)/Vc + Dinf/Vc)
    d/dt(Cp) = (+Q * (Cc - Cp)/Vt)
    Cc_mg_ml = Cc
    Cc_mg_ml ~ add(add_err) + prop(prop_err)
})

} rx_obj = rxode2(rx_fun)#> Error in rxode2(rx_fun): could not find function "rxode2" dataset = et(id=c(1,2), time=c(0,1), evid=0)#> Error in et(id = c(1, 2), time = c(0, 1), evid = 0): could not find function "et"dataset[["dv"]] = NA#> Error: object 'dataset' not foundfor(tmp_cov in rx_obj$allCovs){ dataset[[tmp_cov]] = 1 }#> Error: object 'rx_obj' not found res = nlmixr2(rx_obj, dataset, "poped", babelmixr2::popedControl(modelName="exportModel", runCommand=NA))#> Error in nlmixr2(rx_obj, dataset, "poped", babelmixr2::popedControl(modelName = "exportModel", : could not find function "nlmixr2"

Created on 2024-09-28 with reprex v2.1.1 https://reprex.tidyverse.org

— Reply to this email directly, view it on GitHub https://github.com/nlmixr2/babelmixr2/issues/113#issuecomment-2380721671, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWSJGGNU2RQJQCCE4D3ZY3H7PAVCNFSM6AAAAABL5LAFH2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGOBQG4ZDCNRXGE . You are receiving this because you commented.Message ID: @.***>

mattfidler commented 2 months ago

ie, there is no script generated, it creates a PopED database and loads the rxode2 model into the memory.

john-harrold commented 2 months ago

Wait you said I could do this. Now you're just toying with me :).

mattfidler commented 2 months ago

At the time you could but now it is gone.