Feature Request: Allow rxSolve to take in an un-solved nlmixr function definition

[billdenney]

Related to #548 in that I'm still trying to make that simulation...

It would be helpful if rxSolve could allow the object to be an nlmixr function definition (i.e. a function with an ini() and model()).

An example of the desired functionality is below. If that is out of scope, a clearer error message would help rather than "rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits, : Not compatible with STRSXP: [type=NULL]."

library(nlmixr)

linear_twocmt_iiv_growth <- function() {
  ini({
    tcl <- fixed(log(0.01)) ; label("Clearance (L/kg/hr)")
    tq <- fixed(log(0.05)) ; label("Intercompartmental clearance (L/kg/hr)")
    tv <- fixed(log(1.5)) ; label("Volume central (L/kg)")
    tv1 <- fixed(log(2)) ; label("Volume peripheral (L/kg)")
    
    tktr <- log(c(0.00001, 0.03, 3)) ; label("Transit compartment rate (1/hr)")
    slope <- c(-2, -0.02, 2) ; label("Slope on growth rate ((ng/mL)/hr)")
    ltumorgrowth <- log(c(0.0001, 0.01, 1)) ; label("Baseline tumor exponential growth rate (1/hr)")

    iiv_tumorgrowth ~ 0.1 ; label("Inter-individual variability in baseline tumor exponential growth rate")

    prop_sd <- 0.2 ; label("Proportional error for tumor volume (fraction)")
    add_sd <- 100 ; label("Additive error for tumor volume (mm^3)")
  })
  drake::no_deps(model({
    cl <- exp(tcl)
    q <- exp(tq)
    v <- exp(tv)
    v1 <- exp(tv1)
    kcp <- q/v
    kpc <- q/v1
    ktr <- exp(tktr)
    TUMOR(0) <- TUMORBL
    
    d/dt(CENTRAL) = -cl/v * CENTRAL - kcp*CENTRAL + kpc*PERIPH1
    d/dt(PERIPH1) = kcp*CENTRAL - kpc*PERIPH1
    # unit conversion
    cp <- CENTRAL/v*1000
    d/dt(TRANSIT) = ktr*(cp - TRANSIT)
    tumorgrowth <- exp(ltumorgrowth + iiv_tumorgrowth)
    ktumor <- tumorgrowth + slope*TRANSIT
    d/dt(TUMOR) = ktumor*TUMOR
    TUMOR ~ add(add_sd) + prop(prop_sd)
  }))
}

sim_dose <-
  data.frame(
    ID=1:50,
    time=0,
    amt=rep(seq(0, 20, by=5), each=10),
    evid=1,
    cmt="CENTRAL",
    mdv=1,
    TUMORBL=100
  )

sim_obs <-
  expand.grid(
    ID=1:50,
    time=(0:60)*24,
    amt=0,
    evid=0,
    cmt="TUMOR",
    mdv=0,
    TUMORBL=100
  )

sim_data <- rbind(sim_dose, sim_obs)
sim_data <- sim_data[order(sim_data$ID, sim_data$time), ]

simulation <- rxSolve(linear_twocmt_iiv_growth, events=sim_data)
#> Error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits, : Not compatible with STRSXP: [type=NULL].

^{Created on 2021-07-26 by the reprex package (v2.0.0)}

[mattfidler]

This is what I thought you were requesting in nlmixrdevelopment/RxODE#442 Is this a duplicate of the issue or am I missing something.

[billdenney]

The only marginal/subtle difference is the use of rxSolve() in this issue and nlmxirSim() in the other issue. My guess is that they are duplicates in the end as long as nlmixrSim() gets the same functionality. (Feel free to close; I'm leaving it open for the moment in case you think that the difference between nlmixrSim() and rxSolve() is an actual difference.)

[mattfidler]

The only possible difference is nlmixrSim() reads the information about the model to simulate with uncertainty. That is the number of subjects modeled and the number of observations modeled.

[billdenney]

That sounds possibly like a difference. Up to you if it is a separate issue or not.

My underlying preference is that from a user's perspective, nlmixrSim() and rxSolve() should work almost identically including with unsolved models and with UI-style functions. Or, if they don't work the same (I think that would be because you should be using rxSolve() for unsolved models and UI-style functions and nlmxirSim() for fit models), then the error message should clearly indicate: "use rxSolve() for unsolved models and nlmixr UI functions" or something like that. (But, it seems like it would be easier for everyone to just do that work for the user rather than raise an error.)

[mattfidler]

While I haven't fixed this yet, I have added the ability to save information before a RxODE solve is started. This may make the reprex easier to generate. With the github version of RxODE and nlmixr, you would simply use:

options(RxODE.debug = TRUE)

The solve information will be saved in "last-rxode.qs" which can be read by qs::qread()

This file is a list containing:

• Normalized model text
• RxODE control options
• names internal variable
• extra internal informatioin
• Parameters supplied to RxODE
• Events supplied to RxODE
• Initial conditions
• Boolean to see if this is setup for a in memory solve like with nlmixr

This way you can get the model you are trying to reproduce without the bother of trying to figure out what nlmixr is doing to the model before it gets to where you are at.

I think this backdoor is useful for debugging

[billdenney]

That sounds helpful. I think that this information would be helpful in a "contributing" guideline or some other way to describe how to submit bugs and reprexes.

[mattfidler]

I agree. Documentation is lagging. Maybe even a CRAN type of issue

[dzianismr]

I also think that a function which can simulate an un-solved nlmixr model will be helpful.

[mattfidler]

It will be coming in rxode2 and nlmixr2. I have proof of concept code there