nlmixrdevelopment / RxODE

RxODE is an R package that facilitates easy simulations in R
https://nlmixrdevelopment.github.io/RxODE/
GNU General Public License v3.0
54 stars 14 forks source link

Notation of current time point value in an equation #62

Closed motyocska closed 5 years ago

motyocska commented 5 years ago

Matt,

how would one signify the current time point in the solver if the value of the time point itself is used in an equation within the model in and ODE model using nlmixr? I am looking for the FORTRAN equivalent of T (or t), or with an earlier example of yours from RxODE you may recall the time varying co variate issue with a model where you had an equation as this:

Kin =   Kin0 +amp *cos(2*pi*(ctime-Tz)/24)

how would you refer to ctime as the current time point in the solver? I tried simply enter "time" but will not take it as seems to ask for it to be declared?

appreciate your help

Andras

mattfidler commented 5 years ago

Hi @motyocska,

You have to have ctime in the solving dataset or specify it directly in the covariate In the example you referred to, ctime is in a separate dataset as a time-varying covariate.

library(RxODE)
mod3 <- RxODE({
    KA=2.94E-01;
    CL=1.86E+01; 
    V2=4.02E+01; 
    Q=1.05E+01;
    V3=2.97E+02; 
    Kin0=1;
    Kout=1;
    EC50=200;
    ## The linCmt() picks up the variables from above
    C2   = linCmt();
    Tz= 8
    amp=0.1
    eff(0) = 1  ## This specifies that the effect compartment starts at 1.
    ## Kin changes based on time of day (like cortosol)
    Kin =   Kin0 +amp *cos(2*pi*(ctime-Tz)/24)
    d/dt(eff) =  Kin - Kout*(1-C2/(EC50+C2))*eff;
})

ev <- eventTable(amount.units="mg", time.units="hours") %>%
    add.dosing(dose=10000, nbr.doses=1, dosing.to=2) %>%
    add.sampling(seq(0,48,length.out=100));

 ## Create data frame of  8 am dosing for the first dose
cov.df  <- data.frame(ctime =(seq(0,48,length.out=100)+8) %% 24);
Now there is a covariate present, the system can be solved using the cov option

## Solve with Next Observation Carried Backward (NONMEM's interpolation method)
r1 <- solve(mod3, ev, covs=cov.df,covsInterpolation="nocb")

which gives

r1
▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ 
── Parameters (r1$params): ───────────────────────────────────────────────────── 
        KA         CL         V2          Q         V3       Kin0       Kout 
  0.294000  18.600000  40.200000  10.500000 297.000000   1.000000   1.000000 
      EC50         Tz        amp         pi 
200.000000   8.000000   0.100000   3.141593 
── Covariates (r1$covs): ─────────────────────────────────────────────────────── 
# A tibble: 100 x 1
  ctime
  <dbl>
1  8   
2  8.48
3  8.97
4  9.45
5  9.94
6 10.4 
# ... with 94 more rows
── Initial Conditions (r1$inits): ────────────────────────────────────────────── 
eff 
  1 
── First part of data (object): ──────────────────────────────────────────────── 
# A tibble: 100 x 4
   time    C2   Kin   eff
  <dbl> <dbl> <dbl> <dbl>
1 0       0    1.10  1   
2 0.485  27.8  1.10  1.07
3 0.970  43.7  1.09  1.14
4 1.45   51.8  1.09  1.21
5 1.94   54.8  1.08  1.26
6 2.42   54.6  1.07  1.29
# ... with 94 more rows
▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂

Another option is to have a dataset with both pieces in it


dat <- ev$get.EventTable() ## Get the event table dataset
dat$ctime <- c(8,seq(0,48,length.out=100)+8 %% 24)
## Solve with Last Observation Carried Forward (Monolix/nlmixr's interpolation method)
r1 <- solve(mod3, dat, covsInterpolation="locf")

Which gives the same solution:

> r1
▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ 
── Parameters (r1$params): ───────────────────────────────────────────────────── 
        KA         CL         V2          Q         V3       Kin0       Kout 
  0.294000  18.600000  40.200000  10.500000 297.000000   1.000000   1.000000 
      EC50         Tz        amp         pi 
200.000000   8.000000   0.100000   3.141593 
── Initial Conditions (r1$inits): ────────────────────────────────────────────── 
eff 
  1 
── First part of data (object): ──────────────────────────────────────────────── 
# A tibble: 100 x 4
   time    C2   Kin   eff
  <dbl> <dbl> <dbl> <dbl>
1 0       0    1.1   1   
2 0.485  27.8  1.10  1.07
3 0.970  43.7  1.10  1.15
4 1.45   51.8  1.09  1.22
5 1.94   54.8  1.09  1.27
6 2.42   54.6  1.08  1.31
# ... with 94 more rows
▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
mattfidler commented 5 years ago

I'm assuming this works for you @motyocska; If not let me know.

motyocska commented 5 years ago

Matt, it does solve my initial problem, yes...I did try t also I just did not have it in the area where the diffeqs were set up in the model code, but once included there it all worked... Thanks for the example , I actually was trying to implement exactly what you were showing with the cosin equation... I do continue to run into some trouble though with my real data set I am trying to fit though and thought would ask to see if you have thoughts about it although depending on the status of nlmixr development there may be limitation at this time. Long story short I am trying to fit a likely poor data set ( poor I mean as I have no solid confirmation on the amount of control the investigator had on recording timing of events accurately)  and perhaps some model miss-specifications as there are no standards in this area I can find, but I would think irrespective of that the system should (for the most part) run through and generate some sort of results... after running saem I am getting "near expected' results, but I need to implement time varying covariates, at least try to evaluate, as the physiologic/real life process would require, but based on a previous post I saw saem does not have that capability yet? I also tried to run it with focei (the version without implementation of time varying covariates) but getting an error message on my 3rd output equation during the initial calculation of ETA based prediction error and derivatives. Some of this is certainly over my head, especially from the theoretical point of view with regards to the calculations and how they are done in nlmixr and what may be causing the crash, but I do get reasonable output with SADAPT suggesting perhaps the run should work on other systems? nlmixr is certainly appealing as saem did a very good job on run #1 versus took me a while to get to some likely parameter space that should be searched with SADAPT. The second issue I seem to run into is memory with focei, seems to be eating up a lot, almost never able to run on my 4GB machine... Putting all questions about model structure and appropriateness a side (unless that is the real cause of issues) thought I would ask if you would be willing to give it a spin on your machines to see if your experience is the same? I could share model and data file oof the github email list if that works and if solution can be identified then am happy to post an abbreviated version in the github issues list as may be useful for others, also... thanks Andras  On Wednesday, December 5, 2018, 12:52:15 PM EST, Matthew Fidler notifications@github.com wrote:

I'm assuming this works for you @motyocska; If not let me know.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

mattfidler commented 5 years ago

Hi Andreas,

I would be willing to try to debug the issues; Send me an email if you do not wish to share the data/code on a public forum.

On a side note, FOCEi does support time-varying covariates. Also FOCEi loads the solving space and data into memory, which could cause your memory issues.

mattfidler commented 5 years ago

My current development environment is 8GB instead of 4GB, perhaps that is the memory limitation...?

motyocska commented 5 years ago

Matt, yes that certainly could be the problem with the differences in memory... this is the r script I am using and the data file is attached, the model structure can certainly be simplified to have less parameters but please don't mind that (all the K13 and ups to one constant, etc...) but please don't mind that for now unless really have to: library(nlmixr)library(ggplot2)

data<-read.csv("rdirectory\data1csfcosyn.csv")

one.compartment.saem <- function() {  ini({     tVD<-c(log(4),log(12),log(50))    tamp<-c(log(0.01),log(0.5),log(1.2))    tampKREL<-c(log(0.01),log(0.5),log(1.2))    tampKRELB<-c(log(0.01),log(0.5),log(1.2))    tV30<-c(log(0.01),log(0.025),log(0.05))    tCL<-c(log(0.01),log(4),log(8))    tK12<-log(0.2)    tK21<-log(0.3)    tK13<-c(log(0.0000001),log(0.000006),log(0.00001))    tK31<-c(log(0.0000001),log(0.000006),log(0.00001))    tK14<-c(log(0.0000001),log(0.000006),log(0.00001))    tK41<-c(log(0.0000001),log(0.000006),log(0.00001))    tK15<-c(log(0.0000001),log(0.000006),log(0.00001))    tK51<-c(log(0.0000001),log(0.000006),log(0.00001))    tK16<-c(log(0.0000001),log(0.000006),log(0.00001))    tK61<-c(log(0.0000001),log(0.000006),log(0.00001))    tKRELECFF<-c(log(0.0001),log(0.004),log(0.03))    tKREL0<-c(log(0.0001),log(0.1),log(0.9))    tKRELB0<-c(log(0.0001),log(0.01),log(0.1))    tKSHUNT<-c(log(0.0001),log(0.03),log(0.1))    tV6<-c(log(0.0001),log(0.02),log(0.12))     eta.VD~1    eta.amp~1    eta.ampKREL~1    eta.ampKRELB~1    eta.V30~1    eta.CL~1    eta.K12~1    eta.K21~1    eta.K13~1    eta.K31~1    eta.K14~1    eta.K41~1    eta.K15~1    eta.K51~1    eta.K16~1    eta.K61~1    eta.KRELECFF~1    eta.KREL0~1    eta.KRELB0~1    eta.KSHUNT~1    eta.V6~1     add.err1 <- 0.1    add.err2 <- 3    add.err3 <- 0.1   prop.err1<-0.2   prop.err2<-0.3   prop.err3<-0.2   })  model({     VD<-exp(tVD + eta.VD)    amp<-exp(tamp + eta.amp)    ampKREL<-exp(tampKREL + eta.ampKREL)    ampKRELB<-exp(tampKRELB + eta.ampKRELB)    V30<-exp(tV30 + eta.V30)    CL<-exp(tCL + eta.CL)    K12<-exp(tK12 + eta.K12)    K21<-exp(tK21 + eta.K21)    K13<-exp(tK13 + eta.K13)    K31<-exp(tK31 + eta.K31)    K14<-exp(tK14 + eta.K14)    K41<-exp(tK41 + eta.K41)    K15<-exp(tK15 + eta.K15)    K51<-exp(tK51+ eta.K51)    K16<-exp(tK16 + eta.K16)    K61<-exp(tK61 + eta.K61)    KRELECFF<-exp(tKRELECFF + eta.KRELECFF)    KREL0<-exp(tKREL0 + eta.KREL0)    KRELB0<-exp(tKRELB0 + eta.KRELB0)    KSHUNT<-exp(tKSHUNT + eta.KSHUNT)    V6<-exp(tV6 + eta.V6)        AB=cos(2pi(t-0.0001190476)/0.0002380952)    d/dt(x1) = -((CL/VD)+K12+K13+K14+K15+K16)x1+K21x2+K31x3+K41x4+K51x5+K61x6+(KREL0 exp(ampKREL AB))x7    d/dt(x2) = K12x1-K21x2    d/dt(x3) = K13x1-K31x2+(KRELB0 exp(ampKRELB AB))x4-KRELECFFx3    d/dt(x4) = K14 x1 - K41 x4 - (KRELB0 exp(ampKRELB AB)) x4+KRELECFFx3- (KREL0 exp(ampKREL AB)) x4+(KRELB0 exp(ampKRELB AB))x5- KSHUNTCLAMP x4    d/dt(x5) = K15 x1 - K51 x5 + (KRELB0 exp(ampKRELB AB))x6-(KRELB0 exp(ampKRELB AB)) x5+ (KREL0 exp(ampKREL AB)) x4- (KREL0 exp(ampKREL AB)) x5    d/dt(x6) = K16 x1 - K61 x6 + (KRELB0 exp(ampKRELB AB))x7-(KRELB0 exp(ampKRELB AB)) x6+ (KREL0 exp(ampKREL AB)) x5- (KREL0 exp(ampKREL AB))x6    d/dt(x7) = -(KREL0 exp(ampKREL AB)) x7+ (KREL0 exp(ampKREL AB)) x6- (KRELB0 exp(ampKRELB AB))x7    cp = x1 / VD     ccsf = x4/(V30 exp(amp AB))    ctech = x7/V6       cp ~add(add.err1) + prop(prop.err1)|x1   ccsf ~add(add.err2) + prop(prop.err2)|x4   ctech ~add(add.err3) + prop(prop.err3)|x7    #  cp ~ add(add.err1)|x1 # ccsf~ add(add.err2)|x4#  ctech ~add(add.err3)|x7   })} nlmixr(one.compartment.saem)

fit <- nlmixr(one.compartment.saem, data,est="focei",              #table=tableControl(cwres=TRUE,npde=TRUE),              foceiControl(sigdig = 2, maxInnerIterations = 1000,                           maxOuterIterations = 5000, method = c("lsoda"),  atol = 1e-4, rtol = 1e-4,                           maxstepsOde = 5000000, hmin = 0L, hmax = 0.01, hini = 0.01,                           covsInterpolation = c("constant"),                            derivSwitchTol = NULL))

appreciate any thoughts you have... impact of the time varying covariate CLAMP may also be negligble but have it in just to test the example. the idea with that is the last obs gets carried forward, ie if at time 0 CLAMP=0 and at time 1 CLAMP=1 then between the times 0 and 1 the value should be 0, then at time 1 it switches to 1, hence I selected the "constant" for covsInterpolation  thanks, Andras  On Thursday, December 6, 2018, 5:09:15 PM EST, Matthew Fidler notifications@github.com wrote:

My current development environment is 8GB instead of 4GB, perhaps that is the memory limitation...?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.