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

error in eventTable #75

Closed palar359 closed 5 years ago

palar359 commented 5 years ago

Hi Matt, I have started using RxODE and faced an error message, which I don't know how to turn around:

"Events: Parameters: Error in rxSolve.default(mod, inits, Params, ev) : Need some event information (observation/dosing times) to solve. You can use either 'eventTable' or an RxODE compatible data.frame/matrix."

##################################################################
# model
ode <- "
  Gc(0) = Gss*Vg;
  Ge(0) = Gss;
  Gp(0) = Gss*Vp;
  Gpro(0) = Gss*(Clg + Clgi*Iss);
  Gcm = (Ge/Gss)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+absg50)
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec = Iss*Cli*Gcm*Iabsg;
  ktr = (n+1)/mtt

  d/dt(depot) = transit(n,mtt,bio)-ka*depot; 
  d/dt(Gc) = ka*depot + Gpro + (Q/Vp)*Gp - (Clg + Clgi*Ie + Q)*Gc/Vg;
  d/dt(Gp) = Q*(Gc/Vg - Gp/Vp);
  d/dt(Ge) = keog*Gc/Vg - keog*Ge;
  d/dt(I) = Isec - Cli/Vi*I;
  d/dt(Ie) = keoi*I/Vi - keoi*Ie
"

mod <- RxODE(model=ode, modName = "mod")

# initial conditions for ODEs
inits <- c(0,0,0,0,0,0,0)

# thetas
Params <- c(Gss=0.891, Vg=9.33, Vp=8.56, Q=0.442, 
            Clg=0.0287, Clgi=0.0059, Iss=0.93, Cli=1.22, 
            Vi=6.09, Sinc=0.001, IPRG=1.42, keog=0.0289, keoi=0.0213,  
            ka=0.5, mtt=34.9 , bio=0.811, n=1.27, absg50=14.8, Emax=1.47)
ev <- eventTable(amount.units = "g", time.units = "min")
ev$add.dosing(dose=75, nbr.doses=1, start.time = 0, dosing.to = 1)
ev$add.sampling(seq(0,24,1))

igi <- as.data.frame(rxSolve(mod, inits, Params, ev))

###########################################################

Any help is appreciated.

Regards Arindam

mattfidler commented 5 years ago

Hi Arindam

Thank you for the reproducible example.

You need to change the following:

Here is my R information when trying your model:

library(RxODE)
##################################################################
# model
ode <- "
  Gc(0) = Gss*Vg;
  Ge(0) = Gss;
  Gp(0) = Gss*Vp;
  Gpro(0) = Gss*(Clg + Clgi*Iss);
  Gcm = (Ge/Gss)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+absg50)
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec = Iss*Cli*Gcm*Iabsg;
  ktr = (n+1)/mtt

  d/dt(depot) = transit(n,mtt,bio)-ka*depot;
  d/dt(Gc) = ka*depot + Gpro + (Q/Vp)*Gp - (Clg + Clgi*Ie + Q)*Gc/Vg;
  d/dt(Gp) = Q*(Gc/Vg - Gp/Vp);
  d/dt(Ge) = keog*Gc/Vg - keog*Ge;
  d/dt(I) = Isec - Cli/Vi*I;
  d/dt(Ie) = keoi*I/Vi - keoi*Ie
"

mod <- RxODE(model=ode, modName = "mod")

# initial conditions for ODEs
inits <- c(0,0,0,0,0,0,0)

# thetas
Params <- c(Gss=0.891, Vg=9.33, Vp=8.56, Q=0.442,
            Clg=0.0287, Clgi=0.0059, Iss=0.93, Cli=1.22,
            Vi=6.09, Sinc=0.001, IPRG=1.42, keog=0.0289, keoi=0.0213,
            ka=0.5, mtt=34.9 , bio=0.811, n=1.27, absg50=14.8, Emax=1.47)
ev <- eventTable(amount.units = "g", time.units = "min")
ev$add.dosing(dose=75, nbr.doses=1, start.time = 0, dosing.to = 1)
ev$add.sampling(seq(0,24,1))

igi <- as.data.frame(rxSolve(mod, inits, Params, ev))
#> Events:
#> Parameters:
#> Error in rxSolve_(object, rxControl(..., events = events, params = params), : Need some event information (observation/dosing times) to solve.
#> You can use either 'eventTable' or an RxODE compatible data.frame/matrix.

###########################################################

## The issue here is the order of parameters
## You can fix it by specifying inits=inits

igi <- rxSolve(mod, inits=inits, Params, ev)
#> Error in rxSolve_(object, rxControl(..., events = events, params = params), : Length mismatch
#> req: c(depot, Gc, Gp, Ge, I, Ie)
#> vec: c(0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000)

## But in this case you need to match the number of initial conditions
## specified by your ODE.  You seem to be missing gPro

## OR you can drop the inits argument since they are all equal to zero
## Note you can solve this

igi <- rxSolve(mod, Params, ev)
#> Warning in rxSolve_(object, rxControl(..., events = events, params =
#> params), : Assumed transit compartment model since 'podo' is in the model.
#> Model:
#> 
#> Gc(0)=Gss*Vg;
#> Ge(0)=Gss;
#> Gp(0)=Gss*Vp;
#> Gpro(0)=Gss*(Clg+Clgi*Iss);
#> Gcm=(Ge/Gss)^IPRG;
#> Iabsg=1+(Emax*absg)/(absg+absg50);
#> Ie(0)=Iss;
#> I(0)=Iss*Vi;
#> Isec=Iss*Cli*Gcm*Iabsg;
#> ktr=(n+1)/mtt;
#> d/dt(depot)=transit(n,mtt,bio)-ka*depot;
#> d/dt(Gc)=ka*depot+Gpro+(Q/Vp)*Gp-(Clg+Clgi*Ie+Q)*Gc/Vg;
#> d/dt(Gp)=Q*(Gc/Vg-Gp/Vp);
#> d/dt(Ge)=keog*Gc/Vg-keog*Ge;
#> d/dt(I)=Isec-Cli/Vi*I;
#> d/dt(Ie)=keoi*I/Vi-keoi*Ie;
#> Error in rxSolve_(object, rxControl(..., events = events, params = params), : The following parameter(s) are required for solving: absg, Gpro

## This shows you are missing gPro

Created on 2019-03-29 by the reprex package (v0.2.1)

mattfidler commented 5 years ago

If you have other problems feel free to submit an issue. Always happy to help.

mattfidler commented 5 years ago

If you get the model running I would be interested in seeing it...

palar359 commented 5 years ago

Thanks Matt! I have figured out few other problems with the above code. I will fix them up first and try running the code again. Will get back to you if I face anything unusual and if the model runs successfully I will put it up here in the forum.

Best regards Arindam

palar359 commented 5 years ago

Hi Matt, I have reworked on the model and get it working but struggling to map omegas for each parameter. There'a an omega matrix which applies for Vg, Q and Vi but for other parameters (Clg, Cli, Clgi, mtt) the omegas have fixed values. My question is how to map those omegas in a single block including the matrix for Vg, Q and Vi?


# Jauslin's IGI (ogtt) model
ode <- "
# volumes in L
# Cl in L/min

  kcp = Q/Vg;
  kpc = Q/Vp;
  kg = Clg/Vg;
  kgi = Clgi/Vg;
  ki = Cli/Vi;
  ktr = n/mtt;
  lnfac = log(2.5066)+(n+0.5)*log(n)-n+log(1+1/(12*n));
  absg = exp(log(75000*(0.811/10))+log(ktr) + n*log(ktr*t+0.00001)-ktr*t-lnfac);
  #Gcm1 = ((Ge1/Gss)+0.001)^GPRG;
  Gcm2 = ((Ge2/Gss)+0.001)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+CA50);
  Isec = Iss*Cli*Gcm2*Iabsg;
  Gpro = Gss*(kg + (kgi*Iss))*Vg; 

  depot(0) = 0;
  trn(0) = 0;
  Gc(0) = Gss*Vg;   
  #Ge1(0) = Gss;
  Ge2(0) = Gss;
  Gp(0) = kcp*Gss*Vg/kpc;   
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec(0) = Iss*ki*Vi;

  d/dt(depot) = -absg;
  d/dt(trn) = absg - ka*trn;  
  d/dt(Gc) = ka*trn + Gss*(kg + (kgi*Iss))*Vg - (kg + (kgi*Ie))*Gc - (Gc*kcp - Gp*kpc); 
  d/dt(Gp) = Gc*kcp - Gp*kpc;
  d/dt(Ge2) = keog*(Gc/Vg) - keog*Ge2;
  d/dt(I) = (Iss*ki*Vi*(((Ge2/Gss)+0.001)^IPRG)*Iabsg) - ki*I;
  d/dt(Ie) = keoi*(I/Vi) - keoi*Ie;

  Cg = Gc/Vg   
  Cp = Gp/Vp
  Ci = I/Vi

"

mod <- RxODE(model=ode, modName = "mod")

# initial conditions for ODEs
#ini <- c(depot=0, Gc=0, Gp=0, Ge1=0, Ge2=0, Gpro=0, I=0, Ie=0, Isec=0, Iabsg=1)

# thetas
theta <- c(Vg=9.33, Vp=8.56, Q=0.442, 
           Clg=0.0287, Clgi=0.0059, 
           Iss=9.3, Gss=150,
           Cli=1.22, Vi=6.09, 
           IPRG=1.42, 
           keog=0.0289, keoi=0.0213,  
           mtt=34.9 , n=1.27, 
           Emax=1.47, ka=0.02865, CA50=14.8)

# omega <- matrix(c(0.0887, -0.192, 0.0855, 0.0, 0.73, -0.12, 0.0, 0.0, 0.165),3,3, 
                 dimnames=list(NULL,c("eta.Vg", "eta.Q", "eta.Vi")));

et <- eventTable(amount.units = "mg", time.units = "min");
et$add.dosing(dose=75000, nbr.doses=1, start.time = 0, dosing.to = 1);
et$add.sampling(0:360);

igi <- rxSolve(mod, theta, et)

Regards Arindam

palar359 commented 5 years ago

Of course, I will inculde etas with the typical values of each parameters which I didn't mention in the model presented.

mattfidler commented 5 years ago

We do support a transit compartment. See:

https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-transit-compartments.html

You could simplify the model some by using the transit comparment, and it would also be a little more accurate.

palar359 commented 5 years ago

Hi Matt, Yes, that can be used but the primary problem is creating an omega matrix with so many omegas. In NONMEM it's easy to do as a value can be specified against each parameter and a block can be made wherever there's a correlation established or found. In this problem I am grappling with an omega block (with correlated parameters) plus omegas from other parameters.

I have added the omegas in a very naive way as:

ode <- "
.
..
.
omega1 <- matrix(c(0.0887, -0.192, 0.0855, 0.0, 0.73, -0.12, 0.0, 0.0, 0.165),3,3, 
                 dimnames=list(NULL,c("eta.Vg", "eta.Q", "eta.Vi")));
omega2 <- matrix(0.352,dimnames=list(NULL,c("eta.Clg")));
omega3 <- matrix(0.207,dimnames=list(NULL,c("eta.Clgi")));
omega4 <- matrix(0.0852,dimnames=list(NULL,c("eta.Cli")));
.
.
igi_sim <- rxSolve(mod, theta, omega = c(omega1, omega2, omega3, omega4,
                   omega5, omega6, omega7, omega8, omega9, omega10, omega11,
                   omega12, omega13, omega14, omega15), et, nSub=10)

." but it throws an obvious error message as it is trying to find an overall matrix of omega, which I am struggling with:

"Error in rxSolve.default(mod, theta, omega = c(omega1, omega2, omega3, : Not a matrix."

mattfidler commented 5 years ago
library(RxODE)
set.seed(42)
## Jauslin's IGI (ogtt) model
ode <- "
# volumes in L
# Cl in L/min
  Vg = tVg*exp(eta.Vg);
  Q =tQ*exp(eta.Q);
  Vi = tVi*exp(eta.Vi);
  kcp = Q/Vg;
  kpc = Q/Vp;
  kg = Clg/Vg;
  kgi = Clgi/Vg;
  ki = Cli/Vi;
  ktr = n/mtt;
  lnfac = log(2.5066)+(n+0.5)*log(n)-n+log(1+1/(12*n));
  absg = exp(log(75000*(0.811/10))+log(ktr) + n*log(ktr*t+0.00001)-ktr*t-lnfac);
  #Gcm1 = ((Ge1/Gss)+0.001)^GPRG;
  Gcm2 = ((Ge2/Gss)+0.001)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+CA50);
  Isec = Iss*Cli*Gcm2*Iabsg;
  Gpro = Gss*(kg + (kgi*Iss))*Vg;

  depot(0) = 0;
  trn(0) = 0;
  Gc(0) = Gss*Vg;
  #Ge1(0) = Gss;
  Ge2(0) = Gss;
  Gp(0) = kcp*Gss*Vg/kpc;
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec(0) = Iss*ki*Vi;

  d/dt(depot) = -absg;
  d/dt(trn) = absg - ka*trn;
  d/dt(Gc) = ka*trn + Gss*(kg + (kgi*Iss))*Vg - (kg + (kgi*Ie))*Gc - (Gc*kcp - Gp*kpc);
  d/dt(Gp) = Gc*kcp - Gp*kpc;
  d/dt(Ge2) = keog*(Gc/Vg) - keog*Ge2;
  d/dt(I) = (Iss*ki*Vi*(((Ge2/Gss)+0.001)^IPRG)*Iabsg) - ki*I;
  d/dt(Ie) = keoi*(I/Vi) - keoi*Ie;

  Cg = Gc/Vg
  Cp = Gp/Vp
  Ci = I/Vi

"

mod <- RxODE(model=ode)

# initial conditions for ODEs
#ini <- c(depot=0, Gc=0, Gp=0, Ge1=0, Ge2=0, Gpro=0, I=0, Ie=0, Isec=0, Iabsg=1)

# thetas
theta <- c(tVg=9.33, Vp=8.56, tQ=0.442,
           Clg=0.0287, Clgi=0.0059,
           Iss=9.3, Gss=150,
           Cli=1.22, tVi=6.09,
           IPRG=1.42,
           keog=0.0289, keoi=0.0213,
           mtt=34.9 , n=1.27,
           Emax=1.47, ka=0.02865, CA50=14.8)

omega <- matrix(c(0.0887, -0.192, 0.0855, 0.0, 0.73, -0.12, 0.0, 0.0, 0.165),3,3,
                 dimnames=list(NULL,c("eta.Vg", "eta.Q", "eta.Vi")));

et <- eventTable(amount.units = "mg", time.units = "min");
et$add.dosing(dose=75000, nbr.doses=1, start.time = 0, dosing.to = 1);
et$add.sampling(0:360);

et <- et %>% et(id=1:7)

igi <- rxSolve(mod, theta, et, omega=omega)

plot(igi, "Cp")

Created on 2019-04-10 by the reprex package (v0.2.1)

mattfidler commented 5 years ago

I didn't see your comment. You have a larger matrix for omega, then?

palar359 commented 5 years ago

Yes, the omega matrix will contain omegas for all the parameters listed under "theta". Of them only Vg, Q and Vi are correlated but all of other parameters have respective omegas.

mattfidler commented 5 years ago

The key is to create a block matrix with the Matrix package and then convert to a regular matrix. Unfortunately it drops the names, so you will either name all the parameters at the end or to combine the names at the end (which is what I did). An example follows:

library(RxODE)
set.seed(42)
## Jauslin's IGI (ogtt) model
ode <- "
# volumes in L
# Cl in L/min
  Vg = tVg*exp(eta.Vg);
  Q =tQ*exp(eta.Q);
  Vi = tVi*exp(eta.Vi);
  Clg = tClg*exp(eta.Clg)
  Clgi = tClgi*exp(eta.Clgi)
  Cli = tCli*exp(eta.Cli)
  kcp = Q/Vg;
  kpc = Q/Vp;
  kg = Clg/Vg;
  kgi = Clgi/Vg;
  ki = Cli/Vi;
  ktr = n/mtt;
  lnfac = log(2.5066)+(n+0.5)*log(n)-n+log(1+1/(12*n));
  absg = exp(log(75000*(0.811/10))+log(ktr) + n*log(ktr*t+0.00001)-ktr*t-lnfac);
  #Gcm1 = ((Ge1/Gss)+0.001)^GPRG;
  Gcm2 = ((Ge2/Gss)+0.001)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+CA50);
  Isec = Iss*Cli*Gcm2*Iabsg;
  Gpro = Gss*(kg + (kgi*Iss))*Vg;

  depot(0) = 0;
  trn(0) = 0;
  Gc(0) = Gss*Vg;
  #Ge1(0) = Gss;
  Ge2(0) = Gss;
  Gp(0) = kcp*Gss*Vg/kpc;
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec(0) = Iss*ki*Vi;

  d/dt(depot) = -absg;
  d/dt(trn) = absg - ka*trn;
  d/dt(Gc) = ka*trn + Gss*(kg + (kgi*Iss))*Vg - (kg + (kgi*Ie))*Gc - (Gc*kcp - Gp*kpc);
  d/dt(Gp) = Gc*kcp - Gp*kpc;
  d/dt(Ge2) = keog*(Gc/Vg) - keog*Ge2;
  d/dt(I) = (Iss*ki*Vi*(((Ge2/Gss)+0.001)^IPRG)*Iabsg) - ki*I;
  d/dt(Ie) = keoi*(I/Vi) - keoi*Ie;

  Cg = Gc/Vg
  Cp = Gp/Vp
  Ci = I/Vi

"

mod <- RxODE(model=ode)

# initial conditions for ODEs
#ini <- c(depot=0, Gc=0, Gp=0, Ge1=0, Ge2=0, Gpro=0, I=0, Ie=0, Isec=0, Iabsg=1)

# thetas
theta <- c(tVg=9.33, Vp=8.56, tQ=0.442,
           tClg=0.0287, tClgi=0.0059,
           Iss=9.3, Gss=150,
           tCli=1.22, tVi=6.09,
           IPRG=1.42,
           keog=0.0289, keoi=0.0213,
           mtt=34.9 , n=1.27,
           Emax=1.47, ka=0.02865, CA50=14.8)

omega1 <- matrix(c(0.0887, -0.192, 0.0855, 0.0, 0.73, -0.12, 0.0, 0.0, 0.165),3,3,
                dimnames=list(NULL,c("eta.Vg", "eta.Q", "eta.Vi")));
omega2 <- matrix(0.352,dimnames=list(NULL,c("eta.Clg")));
omega3 <- matrix(0.207,dimnames=list(NULL,c("eta.Clgi")));
omega4 <- matrix(0.0852,dimnames=list(NULL,c("eta.Cli")));

library(Matrix);

omega <- as.matrix(bdiag(list(omega1, omega2, omega3, omega4)))

n  <- unlist(c(dimnames(omega1)[2], dimnames(omega2)[2], dimnames(omega3)[2], dimnames(omega4)[2]))

dimnames(omega) <- list(n, n)

et <- eventTable(amount.units = "mg", time.units = "min");
et$add.dosing(dose=75000, nbr.doses=1, start.time = 0, dosing.to = 1);
et$add.sampling(0:360);

et <- et %>% et(id=1:7)

igi <- rxSolve(mod, theta, et, omega=omega)

plot(igi, "Cp")

Created on 2019-04-10 by the reprex package (v0.2.1)

mattfidler commented 5 years ago

If you update RxODE, you can simply say rxSolve(...,omega=list(omega1,omega2,...))

library(RxODE)
set.seed(42)
## Jauslin's IGI (ogtt) model
ode <- "
# volumes in L
# Cl in L/min
  Vg = tVg*exp(eta.Vg);
  Q =tQ*exp(eta.Q);
  Vi = tVi*exp(eta.Vi);
  Clg = tClg*exp(eta.Clg)
  Clgi = tClgi*exp(eta.Clgi)
  Cli = tCli*exp(eta.Cli)
  kcp = Q/Vg;
  kpc = Q/Vp;
  kg = Clg/Vg;
  kgi = Clgi/Vg;
  ki = Cli/Vi;
  ktr = n/mtt;
  lnfac = log(2.5066)+(n+0.5)*log(n)-n+log(1+1/(12*n));
  absg = exp(log(75000*(0.811/10))+log(ktr) + n*log(ktr*t+0.00001)-ktr*t-lnfac);
  #Gcm1 = ((Ge1/Gss)+0.001)^GPRG;
  Gcm2 = ((Ge2/Gss)+0.001)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+CA50);
  Isec = Iss*Cli*Gcm2*Iabsg;
  Gpro = Gss*(kg + (kgi*Iss))*Vg;

  depot(0) = 0;
  trn(0) = 0;
  Gc(0) = Gss*Vg;
  #Ge1(0) = Gss;
  Ge2(0) = Gss;
  Gp(0) = kcp*Gss*Vg/kpc;
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec(0) = Iss*ki*Vi;

  d/dt(depot) = -absg;
  d/dt(trn) = absg - ka*trn;
  d/dt(Gc) = ka*trn + Gss*(kg + (kgi*Iss))*Vg - (kg + (kgi*Ie))*Gc - (Gc*kcp - Gp*kpc);
  d/dt(Gp) = Gc*kcp - Gp*kpc;
  d/dt(Ge2) = keog*(Gc/Vg) - keog*Ge2;
  d/dt(I) = (Iss*ki*Vi*(((Ge2/Gss)+0.001)^IPRG)*Iabsg) - ki*I;
  d/dt(Ie) = keoi*(I/Vi) - keoi*Ie;

  Cg = Gc/Vg
  Cp = Gp/Vp
  Ci = I/Vi

"

mod <- RxODE(model=ode)

# initial conditions for ODEs
#ini <- c(depot=0, Gc=0, Gp=0, Ge1=0, Ge2=0, Gpro=0, I=0, Ie=0, Isec=0, Iabsg=1)

# thetas
theta <- c(tVg=9.33, Vp=8.56, tQ=0.442,
           tClg=0.0287, tClgi=0.0059,
           Iss=9.3, Gss=150,
           tCli=1.22, tVi=6.09,
           IPRG=1.42,
           keog=0.0289, keoi=0.0213,
           mtt=34.9 , n=1.27,
           Emax=1.47, ka=0.02865, CA50=14.8)

omega1 <- matrix(c(0.0887, -0.192, 0.0855, 0.0, 0.73, -0.12, 0.0, 0.0, 0.165),3,3,
                dimnames=list(NULL,c("eta.Vg", "eta.Q", "eta.Vi")));
omega2 <- matrix(0.352,dimnames=list(NULL,c("eta.Clg")));
omega3 <- matrix(0.207,dimnames=list(NULL,c("eta.Clgi")));
omega4 <- matrix(0.0852,dimnames=list(NULL,c("eta.Cli")));

et <- eventTable(amount.units = "mg", time.units = "min");
et$add.dosing(dose=75000, nbr.doses=1, start.time = 0, dosing.to = 1);
et$add.sampling(0:360);

et <- et %>% et(id=1:7)

igi <- rxSolve(mod, theta, et, omega=list(omega1, omega2, omega3, omega4))

plot(igi, "Cp")

Created on 2019-04-11 by the reprex package (v0.2.1)

palar359 commented 5 years ago

Thank you Matt for the apparent solution. So far I haven't used the Matrix package (may be need to learn through documentation) so was not aware of that but yes I would try the "list" function.

Kind regards Arindam

palar359 commented 5 years ago

Hi Matt, Sorry!

palar359 commented 5 years ago

When I tried to run it in R, it gave me few unusual messages: with, et <- et %>% et(id=1:7) Error in et(., id = 1:7) : could not find function "et"

when, I blocked this line (et <- et %>% et(id=1:7) and included nSub it thrw another message: igi_sim <- rxSolve(mod, theta, et, omega = list(omega1, omega2, omega3, omega4, omega5, omega6, omega7, omega8, omega9, omega10, omega11, omega12, omega13, omega14, omega15), nSub=10) Error in rxSolve.default(mod, theta, et, omega = list(omega1, omega2, : Not compatible with requested type: [type=list; target=double].

Last but not the least, when I tried exactly running the solution you provided it gave back, Error in rxSolve.default(mod, theta, et, omega = list(omega1, omega2, : The following parameter(s) are required for solving: eta.Cli, eta.Clgi, eta.Clg, eta.Vi, eta.Q, eta.Vg

These are totally unexplained!! The same running model is now playing around.

mattfidler commented 5 years ago

Note that I also updated your

You have not updated RxODE from github if the et is not found

palar359 commented 5 years ago

Hi Matt, The model is still not executable! I have run your edited model in R but it still threw the same two errors:

Error in et(., id = 1:7) : could not find function "et"

Error in rxSolve.default(mod, params = theta, events = et, omega = list(omega1, : The following parameter(s) are required for solving: eta.Cli, eta.Clgi, eta.Clg, eta.Vi, eta.Q, eta.Vg

Can you also kindly let me know how should I update RxODE from GitHub? I have followed the installation from https://github.com/nlmixrdevelopment/nlmixr/releases but honestly don't know how to update RxODE.

Regards Arindam

mattfidler commented 5 years ago

My guess is there is two different .libPaths; You can check:

.libPaths()

This can sometimes cause issues when trying to install the releases.

Regardless, you should be able to do the following:

devtools::install_github("nlmixrdevelopment/RxODE")

This should update RxODE to the latest version

mattfidler commented 5 years ago

Once that is done you should be able to see the et function:

> library(RxODE)
> et
function (..., envir = parent.frame()) 
{
    .lst <- as.list(match.call()[-1])
    if (rxIs(.lst[1], "numeric") || rxIs(.lst[1], "integer") || 
        rxIs(.lst[1], "list") || rxIs(.lst[1], "rxEt")) {
        do.call(et.default, .lst, envir = envir)
    }
    else {
        UseMethod("et")
    }
}
<bytecode: 0x000000001a4ddd48>
<environment: namespace:RxODE>
> 
palar359 commented 5 years ago

After uopdate it returned a message: Error in i.p(...) : (converted from warning) installation of package ‘C:/Users/Arindam/AppData/Local/Temp/RtmpuQ1Wpt/file45a420161109/RxODE_0.9.0-5.tar.gz’ had non-zero exit status.

what is this about?

mattfidler commented 5 years ago

This is a new feature of remotes

From the documentation, I suppose you should use:

Sys.setenv("R_REMOTES_NO_ERRORS_FROM_WARNINGS"="true")
devtools::install_github("nlmixrdevelopment/RxODE")

If you want to use the installer route, the newest installer includes these updates:

https://github.com/nlmixrdevelopment/nlmixr/releases

Look at version 1.1.0-7 for a installer that should fix these issues.

mattfidler commented 5 years ago

I developed an easier way to specify matrices, more like nlmixr and similar to NONMEM.

While doing that I noticed your omega matrix is not symmetric, which is required for non cholesky decomposed matrices. It seems to not matter here, though it probably should throw an error.

library(RxODE)
set.seed(42)
## Jauslin's IGI (ogtt) model
ode <- "
# volumes in L
# Cl in L/min
  Vg = tVg*exp(eta.Vg);
  Q =tQ*exp(eta.Q);
  Vi = tVi*exp(eta.Vi);
  Clg = tClg*exp(eta.Clg)
  Clgi = tClgi*exp(eta.Clgi)
  Cli = tCli*exp(eta.Cli)
  kcp = Q/Vg;
  kpc = Q/Vp;
  kg = Clg/Vg;
  kgi = Clgi/Vg;
  ki = Cli/Vi;
  ktr = n/mtt;
  lnfac = log(2.5066)+(n+0.5)*log(n)-n+log(1+1/(12*n));
  absg = exp(log(75000*(0.811/10))+log(ktr) + n*log(ktr*t+0.00001)-ktr*t-lnfac);
  #Gcm1 = ((Ge1/Gss)+0.001)^GPRG;
  Gcm2 = ((Ge2/Gss)+0.001)^IPRG;
  Iabsg = 1 + (Emax*absg)/(absg+CA50);
  Isec = Iss*Cli*Gcm2*Iabsg;
  Gpro = Gss*(kg + (kgi*Iss))*Vg;

  depot(0) = 0;
  trn(0) = 0;
  Gc(0) = Gss*Vg;
  #Ge1(0) = Gss;
  Ge2(0) = Gss;
  Gp(0) = kcp*Gss*Vg/kpc;
  Ie(0) = Iss;
  I(0) = Iss*Vi;
  Isec(0) = Iss*ki*Vi;

  d/dt(depot) = -absg;
  d/dt(trn) = absg - ka*trn;
  d/dt(Gc) = ka*trn + Gss*(kg + (kgi*Iss))*Vg - (kg + (kgi*Ie))*Gc - (Gc*kcp - Gp*kpc);
  d/dt(Gp) = Gc*kcp - Gp*kpc;
  d/dt(Ge2) = keog*(Gc/Vg) - keog*Ge2;
  d/dt(I) = (Iss*ki*Vi*(((Ge2/Gss)+0.001)^IPRG)*Iabsg) - ki*I;
  d/dt(Ie) = keoi*(I/Vi) - keoi*Ie;

  Cg = Gc/Vg
  Cp = Gp/Vp
  Ci = I/Vi

"

mod <- RxODE(model=ode)

# initial conditions for ODEs
#ini <- c(depot=0, Gc=0, Gp=0, Ge1=0, Ge2=0, Gpro=0, I=0, Ie=0, Isec=0, Iabsg=1)

# thetas
theta <- c(tVg=9.33, Vp=8.56, tQ=0.442,
           tClg=0.0287, tClgi=0.0059,
           Iss=9.3, Gss=150,
           tCli=1.22, tVi=6.09,
           IPRG=1.42,
           keog=0.0289, keoi=0.0213,
           mtt=34.9 , n=1.27,
           Emax=1.47, ka=0.02865, CA50=14.8)

omega  <- rxMatrix({
    eta.Vg + eta.Q + eta.Vi ~
        c(0.0887,
          -0.1920, 0.73,
          0.0855, -0.12, 0.165);
    eta.Clg ~ 0.352
    eta.Clgi ~ 0.207
    eta.Cli ~ 0.0852
})

et <- eventTable(amount.units = "mg", time.units = "min");
et$add.dosing(dose=75000, nbr.doses=1, start.time = 0, dosing.to = 1);
et$add.sampling(0:360);

et <- et %>% et(id=1:7)

igi <- rxSolve(mod, theta, et, omega=omega)

plot(igi, "Cp")

Created on 2019-04-25 by the reprex package (v0.2.1)

palar359 commented 5 years ago

Hi Matt, Sorry to bug you!

When I ran it in my Rstudio environment it threw this error: "Error in rxMatrix({ : could not find function "rxMatrix"" I tried to install the latest v1.1.0-7, but it gave me this message: "NSIS error Installer integrity check has failed. Common causes include incomplete download and damaged media. Contact the installer's author to obtain a new copy."

This is not enabling me to install the latest installer. Any idea how to turn around?

Thanks Arindam

mattfidler commented 5 years ago

Hi Arindam,

I think you need to delete the installer and try again. Do you have a 64 bit system? I'm assuming you do since the last installer worked well.

Matt.

palar359 commented 5 years ago

Hi Matt, I have deleted the installer, downloaded the latest version, nlmixr - 1.1.0-9. Updated both R and Rstudio, kept the installer in the same path as that of R. when I ran .libPaths() it showed: [1] "C:/Users/Arindam/Documents/R/win-library/3.6" "c:/R/nlmixr_1.1.0-9/R/library"

and then ran your above code (latest one), but still I am getting this error message: Error in rxMatrix({ : could not find function "rxMatrix".

I have only 64 bit system

mattfidler commented 5 years ago

Hi @palar359

Use lotri instead of rxMatrix

palar359 commented 5 years ago

when I installed and loaded "lotri", it gave the following message:

library(lotri) Error in FUN(X[[i]], ...) : cannot open file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/R/lotri.rdb': No such file or directory In addition: Warning messages: 1: In read.dcf(file.path(p, "DESCRIPTION"), c("Package", "Version")) : cannot open compressed file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/DESCRIPTION', probable reason 'No such file or directory' 2: In FUN(X[[i]], ...) : restarting interrupted promise evaluation

when I tried lotri as: omega <- lotri({ eta.Vg + eta.Q + eta.Vi ~ c(0.0887, -0.1920, 0.73, 0.0855, -0.12, 0.165); eta.Clg ~ 0.352 eta.Clgi ~ 0.207 eta.Cli ~ 0.0852 })

I got this message:

Error: cannot open file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/R/lotri.rdb': No such file or directory In addition: Warning message: restarting interrupted promise evaluation

although installing the package never had a problem.

mattfidler commented 5 years ago

I would remove the package and reinstall

Restart R remove.package("lotri") install.package("lotri") library(lotri)

On Mon, Jun 3, 2019 at 10:52 AM palar359 notifications@github.com wrote:

when I installed and loaded "lotri", it gave the following message:

library(lotri) Error in FUN(X[[i]], ...) : cannot open file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/R/lotri.rdb': No such file or directory In addition: Warning messages: 1: In read.dcf(file.path(p, "DESCRIPTION"), c("Package", "Version")) : cannot open compressed file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/DESCRIPTION', probable reason 'No such file or directory' 2: In FUN(X[[i]], ...) : restarting interrupted promise evaluation

when I tried lotri as: omega <- lotri({ eta.Vg + eta.Q + eta.Vi ~ c(0.0887, -0.1920, 0.73, 0.0855, -0.12, 0.165); eta.Clg ~ 0.352 eta.Clgi ~ 0.207 eta.Cli ~ 0.0852 })

I got this message:

Error: cannot open file 'c:/R/nlmixr_1.1.0-9/R/library/lotri/R/lotri.rdb': No such file or directory In addition: Warning message: restarting interrupted promise evaluation

although installing the package never had a problem.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/RxODE/issues/75?email_source=notifications&email_token=AAD5VWX4DOMRTVSMVUREWYTPYU45JA5CNFSM4HCGYD5KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODWZ3GCQ#issuecomment-498316042, or mute the thread https://github.com/notifications/unsubscribe-auth/AAD5VWQYPQO3R6VAWSEHM7DPYU45JANCNFSM4HCGYD5A .