Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Using monolixModel and nonmemModel attributes #113

Closed
john-harrold opened this issue Aug 2, 2024 · 20 comments
Closed

Using monolixModel and nonmemModel attributes #113

john-harrold opened this issue Aug 2, 2024 · 20 comments

Comments

@john-harrold
Copy link

This is related to this issue:

#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
Copy link
Member

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

  • Evaluating a function creates a non-function object with class
    rxUi. This is what is needed and used for piping, so any piping
    will return this object too.

  • You can always convert back to a function with as.function()
    though you would have to convert it back to access rxUi properties
    and functions like $nonmemModel

@mattfidler
Copy link
Member

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:

  • create a dummy datasets that contains two subjects (could be easily
    created with et()) and then adding the covariates found in
    my_model()$allCovs.

  • Use the nlmixr2 command to create model, datasets, and output 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
Copy link
Author

Awesome thank you.

@john-harrold
Copy link
Author

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 john-harrold reopened this Aug 3, 2024
@john-harrold
Copy link
Author

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

@mattfidler
Copy link
Member

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
Copy link
Member

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

@mattfidler
Copy link
Member

mattfidler commented Aug 3, 2024 •
edited
Loading

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
Copy link
Author

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
Copy link
Member

mattfidler commented Aug 4, 2024 via email

@john-harrold
Copy link
Author

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
Copy link
Member

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

@john-harrold
Copy link
Author

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
Copy link
Member

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

@john-harrold
Copy link
Author

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
Copy link
Author

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
Copy link
Member

mattfidler commented Sep 28, 2024 via email

@mattfidler
Copy link
Member

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

@john-harrold
Copy link
Author

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

@mattfidler
Copy link
Member

At the time you could but now it is gone.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@mattfidler @john-harrold