Quickest way to run same simulation for different ETA parameters?

[rfaelens]

I am using RxODE to perform posthoc estimation of individual parameters. To achieve this, we combine RxODE with optim() to find the maximum likelihood estimate.

It seems using rxSolve() results in quite a bit of extra calculation. I try to read https://github.com/nlmixrdevelopment/RxODE/blob/9378b903020ecfeb28a172bab3f0163afa14e9f1/src/rxData.cpp, but it is quite difficult to follow.

What is the fastest way to perform multiple simulations, only changing parameter values? Parallel evaluation will not help here, as the next ETA to evaluate depends on the current result.

My current way of working (pseudo-code):

m1 <- RxODE(myModel)
et <- eventTable() %>% add.dosing() %>% add.sampling()
observed <- c(18, 12)
fun <- function(parameters) {
  covs <- iovToTimevaryingCovariate(parameters, et$get.EventTable()$time)
  pred <- rxSolve(m1, cbind(et, covs), parameters)
  sum( (pred$CONC - observed)^2 )
}
optim(fun, par=c(ETA1=0, ETA2=0))

[mattfidler]

Hi @rfaelens,

The quickest way is to use nlmixr to optimize; nlmixr(model, data, "posthoc")

The next quickest way is to skip the extra stuff in rxSolve by

rxSolve(m, cbind(et,covs), parameters, returnType="data.frame")

I think I need to add this to the speeding up RxODE vignette.

[mattfidler]

The reason why nlmixr is faster is:

1. It keeps the memory setup for solving
2. It uses n1qn1 and restarts the quadrature and ETAs on multi-individual solves
3. It doesn't even take the time to produce a data frame, but does the objective function calculation in C++

EDIT:

It is also faster because

4. It calculates the exact gradient using sympy/ODE forward senstivities

[mattfidler]

To reduce the time further with RxODE only keep the CONC in your model, that is:

m <- RxODE({
parameter1 ~ ...
d/dt(one) ~...
CONC=
})

[mattfidler]

I was thinking of allowing setup of a data-frame once and the multiple resolving of the setup data with different parameters, but I haven't implemented it yet.

[rfaelens]

That does not make a lot of sense; it seems my use of RxODE is a fringe case anyway. The same issue happens with plenty of other solvers (e.g. deSolve LSODA implementation): the setup step is often not reused for the next evaluation. I will definitely try your suggestions. Another implementation is definitely to use nlmixr; but from a research point of view, we prefer to keep control over the solver methods.

On Sat, Aug 17, 2019 at 3:54 PM Matthew Fidler notifications@github.com wrote:

I was thinking of allowing setup of a data-frame once and the multiple resolving of the setup data with different parameters, but I haven't implemented it yet.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/RxODE/issues/108?email_source=notifications&email_token=AAKQ5AW6OOBYTPV5SVKTAM3QE77JTA5CNFSM4IMOYMOKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4QLZGA#issuecomment-522239128, or mute the thread https://github.com/notifications/unsubscribe-auth/AAKQ5AWCZOESD3SY4WLBEFTQE77JTANCNFSM4IMOYMOA .

[mattfidler]

I'm unsure what you mean by the "keep control over the solver methods"; If you use nlmixr you can choose the solver method and options.

If you mean you wish to compare what deSolve mrgsolve and RxODE produce, that is a different question.

I plan to implement the setup skipping. The RxODE setup is more expensive than say deSolve or likely mrgsolve. This is because of the required data setup to allow parallel solving.

[mattfidler]

I did set this up @rfaelens (I know it has been almost a year), but it was difficult to test/maintain/use, so I took it out.

After more thought I think a better solution is to create a "compiled" event table that requires a model to setup. Hence many setup steps can be skipped and hence a faster run.

[mattfidler]

The etTrans(data,model) can be used to speed up the process slightly

library(RxODE)
library(ggplot2)
library(microbenchmark)

nsubg <- 200 # subjects per dose

rx <- RxODE({
    CL =  log(4)
    V = log(70)
    KA = log(1)
    CL = exp(CL + eta.CL)
    V = exp(V + eta.V)
    KA = exp(KA + eta.KA)
    d/dt(abs)    = -KA*abs;
    d/dt(centr)  =  KA*abs-(CL/V)*centr;
    C2=centr/V;
})

omega <- lotri(eta.CL ~ 0.09,
               eta.V ~ 0.09,
               eta.KA ~ 0.09)

doses <- c(10, 30, 60, 120)

ev <- do.call("rbind",
        lapply(seq_along(doses), function(i){
            et() %>%
                et(amt=doses[i]) %>% # Add single dose
                et(0) %>% # Add 0 observation
                ## Generate 4 samples in 24 hour period
                et(lapply(1:4, function(...){c(0, 24)})) %>%
                et(id=seq(1, nsubg) + (i - 1) * nsubg) %>%
                ## Convert to data frame to skip sorting the data
                ## When binding the data together
                as.data.frame 
        }))

ev2 <- etTrans(ev, rx)

## To better compare, use the same output, that is data.table
mb <- microbenchmark(dt=rxSolve(rx, ev, omega=omega, returnType="data.table"),
                               cmp=rxSolve(rx, ev2, omega=omega, returnType="data.table"))

print(mb)
#> Unit: milliseconds
#>  expr      min        lq     mean   median       uq     max neval
#>    dt 57.57990 229.31958 216.9635 257.6452 264.9177 296.987   100
#>   cmp 22.23381  54.13186 175.1774 224.4673 234.8725 273.955   100

autoplot(mb)
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.

^{Created on 2020-07-10 by the reprex package (v0.3.0)}