usage of EVID

[csetraynor]

Hi, I could not fully understand how should EVID be used when using the pmx_solve_group function family. I see from the example (twocpt_population/twocpt_population.stan) that when EVID=1 and time=0 and the output from that time is an output to x. What would happen with other EVID, for example an EVID to reset a compartment? And is it possible to use the same EVID options that in mrgsolve or is limited to NONMEM? The problem at hand that I have requires to reset some (not necessary all) compartments and add some "amt" to a few compartments, I was usign EVID=8 in mrgsolve but I could not find that solution in NONMEM manual, so step by step is there EVID=8 option available in pmx_solve_group?

Just to clarify I can "implement" EVID=8 type solution by hand using pmx_integrate_ode. But it would be nice and neat (and fast) to have an implementation of EVID=8 in pmx_solve_group so it would be easy to use also with that family of functions.

Just also linking the doc of mrgsolve for this specific usage case: https://mrgsolve.github.io/user_guide/input-data-sets.html

[yizhang-yiz]

Very good point. This shouldn't be hard to implement.

[csetraynor]

Hi yes, would be actually a possibility that integer data was evid? I see that dummy_i is not used and dummy_r is by default rate? Although it may cause more issues than solving any as although it is not clear to me that usage of evid would allow for using statements like:

 if (evid==2 && t==time ){
x[1] = 0;
} 

[yizhang-yiz]

Although it may cause more issues than solving any as although it is not clear to me that usage of evid would allow for using statements like:

 if (evid==2 && t==time ){
x[1] = 0;
} 

This is unlikely a solution, because

• input EVID is not necessarily the final events series fed to solvers because there may be additional events generated internally, e.g. when ADDL events are present.
• a functional argument aka higher-order function in Stan only accept const reference type args, so x passed to ODE function cannot be modified within the function definition.

The way I'd do it is to follow current pattern to check event type. For example, the following code(for MPI solvers) checks event types and acts accordingly https://github.com/metrumresearchgroup/Torsten/blob/de93edcf8cad2b030c7be054cd697af8e34ea674/cmdstan/stan/lib/stan_math/stan/math/torsten/event_solver.hpp#L177-L197 For current issue one could add another branch like else if (events.evid(i) == 8) {...}.

[csetraynor]

Hi Sorry for not replying I have been on the road for the past few days. But I am still looking into this. I tried to add this simple line each time events.evid is used in events_solver as you recommended. https://github.com/csetraynor/Torsten/blob/7a73edd1880a262071eb38edfc82ff6def620a92/cmdstan/stan/lib/stan_math/stan/math/torsten/event_solver.hpp#L160-L162

I also tried a more fancy approach in this branch where I also tweaked the events manager and history files. https://github.com/csetraynor/Torsten/blob/bfead2c5c325940a79614a97ea5e483692320c5f/cmdstan/stan/lib/stan_math/stan/math/torsten/event_solver.hpp#L160 https://github.com/csetraynor/Torsten/blob/bfead2c5c325940a79614a97ea5e483692320c5f/cmdstan/stan/lib/stan_math/stan/math/torsten/events_manager.hpp#L93 https://github.com/csetraynor/Torsten/blob/bfead2c5c325940a79614a97ea5e483692320c5f/cmdstan/stan/lib/stan_math/stan/math/torsten/event_history.hpp#L220

I am running the model using this temporary repo where I also posted the simulated data sample and stan model I am using for testing. https://github.com/csetraynor/kolmo/tree/master/DerivedData https://github.com/csetraynor/kolmo/blob/master/Models/kolmo3.stan

Although it does not cause compilation errors, I am running into problems when sampling the model.

Rank 1 failed to solve id 9: CVode(mem, ts[i], y, &t1, CV_NORMAL) failed with error flag -4

The error is actually probably caused because an error in my modifications. I will keep looking into it but if you find that I am commiting any major mistakes in the modification of the files let me know please.

------ EDIT ------- Turns out I had a very easy to spot mistake in the fancier branch named replace now is fixed and the model samples without warnings but I still need to confirm that the inference is correct i.e. recovers parameters from the simulation. I will need to run at least a sample of 250 individuals for that (I tested parameter recovery with the analytical solution first and 200 individuals is roughly the minimum number to recover the parameters from the inverse sampling method used). Will update on this. The new file looks like: https://github.com/csetraynor/Torsten/blob/cc272c716f9d64d91bb8b729588b0e45ff3ac24c/cmdstan/stan/lib/stan_math/stan/math/torsten/event_solver.hpp#L160

[csetraynor]

An update on this issue I will be submitting my first PR in the next days.

[csetraynor]

21 No problem. This is a very simple snippet, twoCpt.r to quickly test that this is supposed to do what is meant to.

data <- read_rdump("./twoCpt.data.r")
data

data$nEvent <- 53
data$cmt <- c(data$cmt[1:4], 3, data$cmt[5:52])
data$evid <- c(data$evid[1:4], 8, data$evid[5:52])  
data$addl  <- c(data$addl,0)
data$ss  <- c(data$ss, 0)
data$amt <-  c(data$amt[1:4], 1000, data$amt[5:52])  
data$time <- c(data$time[1:4], data$time[4], data$time[5:52])  
data$rate <-  c(data$rate, 0)
data$ii <-  c(data$ii, 0) 

Results will not make sense but just using print option an 1 iteration 1 warmup shows that is doing what is supposed to.

  for(n in 1:nEvent){
    print("time = ", time[n]);
    print("EVID = ", evid[n]);
    print("DV =", mass[ ,n]);
  }  

I will also look into posting or adding an example case where this option is useful.

[yizhang-yiz]

Solved by https://github.com/metrumresearchgroup/torsten_math/issues/1.