andrewhooker / poped Goto Github PK
View Code? Open in Web Editor NEWPopulation Experimental Design (PopED) in R
Home Page: https://andrewhooker.github.io/PopED/
License: GNU Lesser General Public License v3.0
Population Experimental Design (PopED) in R
Home Page: https://andrewhooker.github.io/PopED/
License: GNU Lesser General Public License v3.0
Hi Andy,
Would it be possible to implement parallel computation in plot_efficiency_of_windows()
? I see some commented out references to something called bParallel
in the code, so I wonder if this is something you've already considered.
Thanks,
Tim
Hi Andrew,
It seems that how to implement IV infusion with ODE models is not documented in PopEd or Desolve packages. Could you please shed some light on that?
Thank you,
Sibo
When trying to visualize the optimized sample times on top of the predicted curve, this points seem to be off. To me it looks like a problem with the times.
`library(PopED)
library(deSolve)
sfg <- function(x,a,bpop,b,bocc){
parameters=c(CL=bpop[1]*exp(b[1]),
V=bpop[2]*exp(b[2]),
KA=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1])
return(parameters)
}
system("R CMD SHLIB one_comp_oral_CL.c")
dyn.load(paste("one_comp_oral_CL", .Platform$dynlib.ext, sep = ""))
ff.ODE.compiled <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
A_ini <- c(A1 = DOSE*Favail, A2 = 0)
times <- drop(xt)##xt[,,drop=T]
times <- sort(times)
times <- c(0,times) ## add extra time for start of integration
out <- ode(A_ini, times, func = "derivs", parms = c(CL,V,KA),
#jacfunc = "jac", # not really needed, speed up is minimal if this is defined or not.
dllname = "one_comp_oral_CL",
initfunc = "initmod", nout = 1, outnames = "Sum")
y = out[,"A2"]/(V)
y=y[-1] # remove initial time for start of integration
y = cbind(y) # must be a column matrix
return(list( y= y,poped.db=poped.db))
})
}
poped.db.compiled <- create.poped.database(ff_fun=ff.ODE.compiled,
fg_fun = sfg,
fError_fun = feps.add.prop,
bpop=c(CL=0.15, V=8, KA=1.0, Favail=1),
notfixed_bpop=c(1,1,1,0),
d=c(CL=0.07, V=0.02, KA=0.6),
sigma=c(0.01,0.25),
groupsize=32,
xt=c( 0.5,1,2,6,24,36,72,120),
minxt=0,
maxxt=120,
a=c(DOSE=70),
mina=c(DOSE=0),
maxa=c(DOSE=100))
plot_model_prediction(poped.db.compiled,IPRED=T,DV=T)
output <- poped_optim(poped.db.compiled ,opt_xt=T, opt_a=T, parallel=T, dlls = c('one_comp_oral_CL'))
plot_model_prediction(output$poped.db,IPRED=T,DV=T)`
Hi,
Is it possible to implement the low limit of quantification so that the popED consider the drug concentration below as missing value.
Thank you,
Mark,
Monash University
Hello,
I am having an issue with the following code.
I have two poped.db that differ only by 1 sampling time (xt). poped.db one having an additional sample at t=10^-5.
poped.db$design$xt
obs_1 obs_2 obs_3 obs_4
grp_1 1e-05 10.00001 20.00001 22.50001
grp_2 1e-05 10.00001 20.00001 22.50001
grp_3 1e-05 10.00001 20.00001 22.50001
...
poped.db2$design$xt
obs_1 obs_2 obs_3
grp_1 10.00001 20.00001 22.50001
grp_2 10.00001 20.00001 22.50001
grp_3 10.00001 20.00001 22.50001
...
So in theory, the ofv of poped.db should be equal or better than the ofv of poped.db2 but in practice :
evaluate_design(poped.db)
$ofv
[1] 43.91572
$fim
KNET BMAX INOC ESLOPE_PMB_NR GAMMA_PMB SIGMA[1,1]
KNET 53703.510 6729.2907 2715.7226 -54735.168 71329.543 0.000
BMAX 6729.291 984.2588 390.4884 -6851.221 10780.391 0.000
INOC 2715.723 390.4884 828.8236 -2788.715 3319.783 0.000
ESLOPE_PMB_NR -54735.168 -6851.2210 -2788.7153 55823.678 -70946.687 0.000
GAMMA_PMB 71329.543 10780.3908 3319.7827 -70946.687 203204.863 0.000
SIGMA[1,1] 0.000 0.0000 0.0000 0.000 0.000 6840.097
$rse
KNET BMAX INOC ESLOPE_PMB_NR GAMMA_PMB SIGMA[1,1]
45.5228679 1.8679539 0.5727989 40.7431924 39.1848724 6.4549722
evaluate_design(poped.db2)
$ofv
[1] 67.28574
$fim
KNET BMAX INOC ESLOPE_PMB_NR GAMMA_PMB SIGMA[1,1]
KNET 363225.581 9204.1376 5.843678e+06 -27047.92 71330.730 0.000
BMAX 9204.138 995.3279 4.708744e+04 -2952.73 10780.763 0.000
INOC 5843678.039 47087.4356 1.102244e+08 521117.71 3319.722 0.000
ESLOPE_PMB_NR -27047.917 -2952.7295 5.211177e+05 691485.22 -70947.776 0.000
GAMMA_PMB 71330.730 10780.7627 3.319722e+03 -70947.78 203223.119 0.000
SIGMA[1,1] 0.000 0.0000 0.000000e+00 0.00 0.000 5130.072
$rse
KNET BMAX INOC ESLOPE_PMB_NR GAMMA_PMB SIGMA[1,1]
1.22385846 1.49656707 0.01010109 0.11616444 13.61248583 7.45355992
Could there be a problem with comptation of the fim when xt have a lot of significant digits ?
Full code:
library(PopED)
library(ggplot2)
library(deSolve)
set.seed(1998)
sfg <- function(x,a,bpop,b,bocc){
parameters=c(KA = bpop[1], #no IIV
Vma= bpop[2],
Ama= bpop[3],
CL = bpop[4],
Vme= bpop[5],
Cm = bpop[6],
V1 = bpop[7],
V2 = bpop[8],
Q = bpop[9],
KNET = bpop[10],
BMAX = bpop[11],
INOC = bpop[12],
ESLOPE_PMB_NR = bpop[13],
GAMMA_PMB = bpop[14],
DOSE=a[1],
TAU=a[2])
return(parameters)
}
PKPD.resistant.ip.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
CFREE = A3/V1*0.086 #free concentration, 0.086 as fraction unbound
PLATEAU = 1 - (A5/(10**BMAX))
KILL_PMB_NR = 0
if(CFREE > 0) KILL_PMB_NR = (ESLOPE_PMB_NR * (CFREE**GAMMA_PMB))
dA1 <- -Vma*A1/(A1+Ama)
dA2 <- -KA*A2
dA3 <- Vma*A1/(A1 + Ama) + KA*A2 -
Q/V1*A3 + Q/V2*A4 -
Vme*A3/(A3+(Cm*V1)) - CL/V1*A3
dA4 <- Q/V1*A3 - Q/V2*A4
dA5 <- (KNET*PLATEAU - KILL_PMB_NR) * A5
return(list(c(dA1, dA2, dA3, dA4, dA5)))
})
}
ff.PKPD.resistant.ip.ode <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
# initial conditions of ODE system
A_ini <- c(A1=0, A2=0, A3=0, A4=0, A5=(10**INOC))
#Set up time points to get ODE solutions
times_xt <- drop(xt) # sample times
times_start <- c(0) # add extra time for start of study
times_dose = seq(from=0,to=max(times_xt),by=TAU) # dose times
times <- unique(sort(c(times_start,times_xt,times_dose))) # combine it all
# Dosing
dose_set <- c()
for (i in DOSE)
{
if (i == 0)
{
dose_set <- append(dose_set, 0)
dose_set <- append(dose_set, 0)
} else
if (i > 0 & i < 4.8) {
dose_set <- append(dose_set, i)
dose_set <- append(dose_set, 0)
} else {
dose_set <- append(dose_set, 4.8)
dose_set <- append(dose_set, i - 4.8)
}
}
dose_dat <- data.frame(
var = c("A1","A2"),
time = rep(times_dose, each = 2), #dose administration in two compartment (A1 and A2) at the same time
value = c(dose_set),
method = c("add")
)
out <- ode(A_ini, times, PKPD.resistant.ip.ode, parameters,
events = list(data = dose_dat))#atol=1e-13,rtol=1e-13)
pd <- out[,"A5"]
y = log10(pd)
# extract the time points of the observations
y = y[match(times_xt,out[,"time"])]
y = cbind(y)
return(list(y=y,poped.db=poped.db))
})
}
feps_ODE_compiled <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- ff.PKPD.resistant.ip.ode(model_switch,xt,parameters,poped.db)
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = y + epsi[,1]
return(list(y=y,poped.db=poped.db))
}
poped.db <- create.poped.database(ff_file="ff.PKPD.resistant.ip.ode",
fg_file="sfg",
fError_file="feps_ODE_compiled",
bpop=c(KA = 1.83,
Vma= 11.3,
Ama= 11.5,
CL = 0.178,
Vme= 0.255,
Cm = 0.823,
V1 = 0.325,
V2 = 0.808,
Q = 0.198,
KNET = 1.11,
BMAX = 7.453,
INOC = 6.81,
ESLOPE_PMB_NR = 1.216,
GAMMA_PMB = 0.02626),
notfixed_bpop=c(rep(0, 9),1,1,1,1,1),
sigma=0.4328^2,
#design
groupsize = 4,
m = 30,
a =list(c(DOSE=0, TAU=24),
c(DOSE=0.5, TAU=24), c(DOSE=10, TAU=24),
c(DOSE=20, TAU=24), c(DOSE=30, TAU=24), c(DOSE=45, TAU=24),
c(DOSE=5, TAU=12), c(DOSE=10, TAU=12), c(DOSE=15, TAU=12),
c(DOSE=22.5, TAU=12), c(DOSE=30, TAU=12), c(DOSE=45, TAU=12),
c(DOSE=0.83, TAU=8), c(DOSE=1.67, TAU=8), c(DOSE=3.33, TAU=8),
c(DOSE=5, TAU=8), c(DOSE=6.67, TAU=8), c(DOSE=7.5, TAU=8),
c(DOSE=10, TAU=8), c(DOSE=15, TAU=8), c(DOSE=20, TAU=8),
c(DOSE=25, TAU=8), c(DOSE=30, TAU=8), c(DOSE=40, TAU=8),
c(DOSE=1.67, TAU=4), c(DOSE=3.33, TAU=4), c(DOSE=5, TAU=4),
c(DOSE=7.5, TAU=4), c(DOSE=10, TAU=4), c(DOSE=15, TAU=4)),
mina = c(DOSE=0,TAU=4),
maxa = c(DOSE=45, TAU=24),
#Design Space
xt = c(10^-5, 10+10^-5, 20+10^-5, 22.5+10^-5), discrete_xt = list(10^-5, seq(0.5+10^-5, 24+10^-5, 0.5), seq(0.5+10^-5, 24+10^-5, 0.5),22.5+10^-5),
bUseGrouped_xt=TRUE)
poped.db2 <- create.poped.database(ff_file="ff.PKPD.resistant.ip.ode",
fg_file="sfg",
fError_file="feps_ODE_compiled",
bpop=c(KA = 1.83,
Vma= 11.3,
Ama= 11.5,
CL = 0.178,
Vme= 0.255,
Cm = 0.823,
V1 = 0.325,
V2 = 0.808,
Q = 0.198,
KNET = 1.11,
BMAX = 7.453,
INOC = 6.81,
ESLOPE_PMB_NR = 1.216,
GAMMA_PMB = 0.02626),
notfixed_bpop=c(rep(0, 9),1,1,1,1,1),
sigma=0.4328^2,
#design
groupsize = 4,
m = 30,
a =list(c(DOSE=0, TAU=24),
c(DOSE=0.5, TAU=24), c(DOSE=10, TAU=24),
c(DOSE=20, TAU=24), c(DOSE=30, TAU=24), c(DOSE=45, TAU=24),
c(DOSE=5, TAU=12), c(DOSE=10, TAU=12), c(DOSE=15, TAU=12),
c(DOSE=22.5, TAU=12), c(DOSE=30, TAU=12), c(DOSE=45, TAU=12),
c(DOSE=0.83, TAU=8), c(DOSE=1.67, TAU=8), c(DOSE=3.33, TAU=8),
c(DOSE=5, TAU=8), c(DOSE=6.67, TAU=8), c(DOSE=7.5, TAU=8),
c(DOSE=10, TAU=8), c(DOSE=15, TAU=8), c(DOSE=20, TAU=8),
c(DOSE=25, TAU=8), c(DOSE=30, TAU=8), c(DOSE=40, TAU=8),
c(DOSE=1.67, TAU=4), c(DOSE=3.33, TAU=4), c(DOSE=5, TAU=4),
c(DOSE=7.5, TAU=4), c(DOSE=10, TAU=4), c(DOSE=15, TAU=4)),
mina = c(DOSE=0,TAU=4),
maxa = c(DOSE=45, TAU=24),
#Design Space
xt = c(10+10^-5, 20+10^-5, 22.5+10^-5), discrete_xt = list(seq(10^-5, 24+10^-5, 0.5), seq(10^-5, 24+10^-5, 0.5), seq(10^-5, 24+10^-5, 0.5)),
bUseGrouped_xt=TRUE)
eval_1 <- evaluate_design(poped.db)
poped.db$design$xt
eval_2 <- evaluate_design(poped.db2)
poped.db2$design$xt
Installing the latest dev version from github fails with the following error -
* installing *source* package ‘PopED’ ...
** R
** inst
** tests
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded
**Error in get(name, envir = asNamespace(pkg), inherits = FALSE) :
object 'checkCompilerOptions' not found
Calls: ::: -> get
Execution halted
ERROR: loading failed**
The CRAN version installs fine.
Hi Professor Hooker:
I am new to your package. I was trying to prepare a template using RxODE and encountered two unexpected errors. Could you please help me?
The error is in bold
plot_model_prediction(design$poped.db, model_num_points = 500, PI = T)
Error in model_prediction(poped.db, model_num_points = model_num_points, :
unused argument (PI = T)
In addition, can we change the display of the plot? eg. color, label, legend title, etc
In the output, the parameters were not named as defined in create.poped.database(d = c(KA = 0.09, CL = 0.25^2, V = 0.09)) as in your example
optimized.design$rse
bpop[1] bpop[2] bpop[3] D[1,1] D[2,2] D[3,3]
13.181201 4.887268 9.723204 95.296309 33.934192 51.834468
SIGMA[1,1] SIGMA[2,2]
39.510102 26.998621
The full code is below:
`library(PopED)
library(RxODE)
sessionInfo()
########################################################################################################
########################################################################################################
mod = RxODE({
CP = CENT/V;
d/dt(DEPOT) = -KADEPOT;
d/dt(CENT) = KADEPOT - (CL/V)*CENT;
})
ff = function(model_switch, xt, p, poped.db){
times_xt = drop(xt)
ev = et(time = 0, amt = p[["DOSE"]], ii = p[["TAU"]], until = max(times_xt))%>%
et(times_xt)
out = rxSolve(mod, p, ev)
return(list(y = matrix(out$CP, ncol =1), poped.db = poped.db))
}
sfg = function(x, a, bpop, b, bocc){
parameters = c(
KA = bpop[1]*exp(b[1]),
CL = bpop[2]*exp(b[2]),
V = bpop[3]*exp(b[3]),
DOSE = a[1],
TAU = a[2]
)
return(parameters)
}
feps = function(model_switch, xt, parameters, epsi, poped.db){
y = do.call(poped.db$model$ff_pointer, list(model_switch, xt, parameters, poped.db))[[1]] # simulate to obtain F
y = y*(1+epsi[,1])+epsi[,2]
return(list(y = y, poped.db = poped.db))
}
poped_db = create.poped.database(ff_fun = ff,
fg_fun = sfg,
fError_fun = feps,
bpop = c(KA = 0.25, CL = 3.75, V = 72.8), # population parameter
d = c(KA = 0.09, CL = 0.25^2, V = 0.09), # BSV variability in variance
sigma = c(prop = 0.04, add = 0.0025),
m = 2, # 2 groups
groupsize = 20, # each group has 20 subjects
xt = c(1, 2, 8, 240, 245), # initial design of sampling time
minxt = c(0, 0, 0, 240, 240), # 5 sampling timepoints, minimum for each
maxxt = c(10, 10, 10, 248, 248), # 5 sampling timespoints, maximum for each
bUseGrouped_xt = 1, # set both groups have same sampling time to true
a = cbind(DOSE = c(20, 40), TAU = c(24, 24)), # initial dosing for group 1 (dose 20, tau 24) and group 2 (dose 40, tau 24)
maxa = c(DOSE = 200, TAU = 24), # maximum dose and dosing interval
mina = c(DOSE = 0, TAU = 24)) # minimum dose and dosing interval
pred = model_prediction(poped_db,
model_num_points = 500, # result in 505 observation for each group (+5 xt time points)
include_sample_times = T)
plot_model_prediction(poped_db,
model_num_points = 500,
separate.groups = T)
plot_model_prediction(poped_db,
model_num_points = 500,
separate.groups = T,
PI = T)
intial.design = evaluate_design(poped_db)
intial.design$ofv
intial.design$rse
design = poped_optim(poped_db, opt_xt = T)
plot_model_prediction(design$poped.db, model_num_points = 500)
optimized.design = evaluate_design(design$poped.db)
optimized.design$ofv
optimized.design$rse
`
sessioninfo as below:
Hi Andy,
The sample code d=c(V=0.09,KA=0.09,CL=0.25^2) defines the inter-individual variability well, which is then fed to the desolver to run the simulation. However, is there anyway to define parameter correlation, eg. IIV for Cl and Vd.
Thanks.
Mark
Hi Andy,
I couldn't figure out where the optimization was parallelized. Some of the setup time with RxODE allows threaded parallelization which I was hoping to tap into.
Can you point me in the right direction here?
Matt
Hello Andrew,
Can I specify no autocorrelation to avoid sampling at the same time? According to the online guide, I tried the following script but got no success.
poped_db <- create.poped.database (strAutoCorrelationFile ="”),
I still have a lot of sampling points clustered at certain times.
Thank you,
Jin
PKPD Modeler at AbbVie
Hi Andy,
We're working through using Rmpi under snow, and noted that there's a lot of explicit calls to parallel:: that I think should not be explicit. If they're not explicit, then the overloads from snow:: should work properly.
Would you be willing to consider a pull request to change the behavior to allow Rmpi as the parallel backend?
Warm Regards,
Mike Dodds
library(snow)
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
Hi,
Can you help me to define a delayed first order absorption model in PopED?
Thank you
Alejandra
Hi @andrewhooker
Hope you doing well.
Do you have any experience using popEd if we can come up with an optimal sampling schedule for a trial comparing 3-different drugs.
Any reference or codes would be helpful
Thanks and regards,
-Ankur Sharma
Hi all,
When planning a PK study over a period of 140 days, I aim to reduce the number of blood samples from 16 time points to a smaller number.
I set opt_samps=T in poped_optim(), but then I get an error message:
poped_optim(poped.db.2cmt.RUVPA, method=c("LS"),parallel=T, opt_xt = F, opt_samps =T)
Error in get_par_and_space_optim(poped.db, opt_xt = opt_xt, opt_a = opt_a, :
Sample number optimization is not yet implemented in the R-version of PopED.
Is there a way to optimize the number of samples (or even better optimize both the number of samples and the sampling times at the same time) in popED?
Many thanks for your kind help.
Kind regards,
Eric
Here is my popED script:
cppFunction('List two_comp_ode(double Time, NumericVector A, NumericVector Pars) {
int n = A.size();
NumericVector dA(n);
double CL = Pars[0];
double V1 = Pars[1];
double Q = Pars[2];
double V2 = Pars[3];
double DOSE = Pars[4];
double TINF = Pars[5];
double TAU = Pars[6];
double In;
if ( fmod(Time,TAU) < TINF ){
In = DOSE/TINF;
}
else {
In = 0;
}
dA[0] = In - (CL/V1)*A[0] + (Q/V2)*A[1] - (Q/V1)*A[0];
dA[1] = - (Q/V2)*A[1] + (Q/V1)*A[0];
return List::create(dA);
}')
sfg2 <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]),
V1=bpop[2]*exp(b[2]),
Q=bpop[3],
V2=bpop[4],
DOSE=a[1],
TINF=a[2],
TAU=a[3])
return( parameters )
}
PK_2comp_ode_inf_rcpp_ff <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0)
times_xt <- drop(xt)
times <- c(0,times_xt) ## add extra time for start of integration
times <- sort(times)
times <- unique(times) # remove duplicates
out <- ode(A_ini, times, two_comp_ode, parameters)#, events = list(data = eventdat)) #atol=1e-13,rtol=1e-13)
y = out[, "A1"]/(V1)
y=y[match(times_xt,out[,"time"])]
y=cbind(y)
return(list(y=y,poped.db=poped.db))
})
}
numsubject = 6 # number of subject per dose group
BW = 80 #kg
DOSE_mg = 3 *BW
Infusion_duration = 2 /24 # mg/day
sampletime_original = c(2/24, 4/24, 6/24, 8/24, 12/24, 24/24, 36/24, 2, 3, 4, 6, 14, 28, 42, 56, 84) # original 16 samples
poped.db.2cmt.RUVPA <- create.poped.database(ff_file = "PK_2comp_ode_inf_rcpp_ff",
fError_file="feps.add.prop",
fg_file="sfg2",
groupsize=numsubject,
m=1, #number of groups
sigma=c(0.15,0.1),
notfixed_sigma = c(0,0),
bpop=c(CL=3.2/1000BW,V1=41.5/1000BW,Q=17.3/1000BW,V2=49.4/1000BW),
notfixed_bpop=c(1,1,1,1),
d=c(CL=0.09,V1=0.09),
notfixed_d=c(1,1),
xt=sampletime_original, #original 16 samples
minxt=2/24 , # original minxt 2/24
maxxt=140 , # original maxxt 140
a=c(DOSE=DOSE_mg,TINF=Infusion_duration,TAU=280),
dSeed = 123456)
Hi Andrew,
I found those two arguments called very similar inputs, particularly for compiled code. The HCV example uses 'ff_file=ff_ODE_compiled", while the PK.1.comp.oral example uses the "ff_fun=ff.ode.rcpp". In the online function reference, it seems that ff_file includes one more element "the function name or filname". Could you please clarify that?
Thank you,
Sibo
The term "model _switch" is included in a couple of example files.
In ex1, the "model _switch" is not defined.
In ex3b, it is defined as model_switch=c(1,1,1,1,1,2,2,2,2,2).
Is there any special purpose for this term?
Thanks.
Hello,
Thank you for this great tool ! I was trying to optimize some dosing regimens using an mrgsolve model as input for PopED. My OS is Windows 10. It worked fine while mono-threading but when paralellizing i did get this annoying error message :
There was a problem accessing the model shared object.
Either the model object is corrupted or
the model was not properly compiled and/or loaded.
Check mrgsolve:::funset(mod) for more information.
After some research I stumbled upon this issue in the mrgsolve github (metrumresearchgroup/mrgsolve#471) where the authors suggested to use the loadso()
command to load the model. After trying multiple different places to put this commnad (in the ff.model function for example, makes R crash on my computer), I finally found that by editing the start_parallel script adding these few lines made the parallelization work like a charm :
# load mrgsolve models in workers using loadso
if (!is.null(mrgsolve_model)) {
parallel::clusterCall(cl, loadso, x=mrgsolve_model)
}
I forked the repository (https://github.com/Vincent-AC/PopED) and added this to the package, would you be interested in a pull request ?
Hi Andy
Can a model be defined as a system of ODEs, since the example mentioned is for a one compartment model, analytical solution?
If so, could you please provide an example for a 2 cmgt model, oral absorption.
Thanks in advance
Amit Taneja
It seems that when you specify a value of 0 as minimum of your design space, it gets replaced by 1e-05
poped.db$design_space$minxt (done with the vignette example)
obs_1 obs_2 obs_3 obs_4 obs_5
grp_1 1e-05 1e-05 1e-05 240 240
grp_2 1e-05 1e-05 1e-05 240 240
However if in your initial design you have a time 0, it does not get replaced
poped.db$design$xt
obs_1 obs_2 obs_3 obs_4 obs_5
grp_1 0 2 8 240 245
grp_2 0 2 8 240 245
Apparently, this is not a problem apparently for a continous space, but if instead you provide a discrete space, you get an error because the obs_1 is out of the design space:
poped.db.discrete <- create.poped.database(poped.db,discrete_xt = list(0:248))
Error in eval(substitute(expr), data, enclos = parent.frame()) :
xt value for group 1 (column 1) is not in the design space
Hi @andrewhooker
I am trying to optimize a PIII trial design using popED.
This trial is for one year and the number of subjects stays same (N=200) under week 52, but post-wk-52 the extensive PK is collected only for N=20 subjects
Before wk-52, trough samples are only collected
Do you know if this situation can be implemented in popED for design evaluation and optimization
Will appreciate your response
Thanks,
-Ankur
Hi all
I was wondering if there is an example of a model with measurements in both plasma and urine?
Thanks in advance!
I am having an issue in implementing a multiple dose TMDD model estimated in nonmem to optimize over 4 analytes in PopED. The error message, verbatim pasted below indicates, a needed missing value of an attribute.
Run the attached code. Appreciate any help to identify and fix the needed attribute.
Nitin
Problems inverting the matrix. Results could be misleading.
Error in if (any(ret == 0)) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In sqrt(param_vars[i]) : NaNs produced
I have encountered a problem and need to use the 3-compartment oral model for modeling, but I did not find any example of this model on the Internet, so I would like to ask you if there is any problem with my code. If there is a problem, I hope you can help me modify it, thank you very much.
Here is my code:
#' An implementation of a two compartment model with oral absorption using ODEs.
library(PopED)
library(deSolve)
#' Define the ODE system
PK.3.comp.oral.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
dA1 <- -KA*A1
dA2 <- KA*A1 + A3* Q3/V3 -A2*(CL/V2+Q3/V2+Q4/V2)+A4*Q4/V4
dA3 <- A2* Q3/V2-A3* Q3/V3
dA4 <- A2* Q4/V2-A4* Q4/V4
return(list(c(dA1, dA2, dA3, dA4)))
})
}
#' define the initial conditions and the dosing
ff.PK.3.comp.oral.md.ode <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0, A3=0, A4=0)
times_xt <- drop(xt)
dose_times = seq(from=0,to=max(times_xt),by=TAU)
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE), method = c("add"))
times <- sort(c(times_xt,dose_times))
out <- ode(A_ini, times, PK.3.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = out[, "A2"]/(V2/Favail)
y=y[match(times_xt,out[,"time"])]
y=cbind(y)
return(list(y=y,poped.db=poped.db))
})
}
#' parameter definition function
#' names match parameters in function ff
fg <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]),
V2=bpop[2]*exp(b[2]),
V3=bpop[3]*exp(b[3]),
V4=bpop[4]*exp(b[4]),
KA=bpop[5]*exp(b[5]),
Q3=bpop[6]*exp(b[6]),
Q4=bpop[7]*exp(b[7]),
Favail=bpop[8]*exp(b[8]),
ALAG1=bpop[9]*exp(b[9]),
DOSE=a[1],
TAU=a[2])
return( parameters )
}
#' create poped database
poped.db <- create.poped.database(ff_fun="ff.PK.3.comp.oral.md.ode",
fError_fun="feps.add.prop",
fg_fun="fg",
groupsize=20,
m=1, #number of groups
sigma=c(prop=0.0714^2,add=6.13^2),
bpop=c(CL=42.9, V2=102, V3=2340, V4=421, KA=3.54, Q3=105, Q4=397, Favail=1, ALAG1=0.0626), (bpop)
d=c(CL=0.487^2, V2=0.385^2, V3=0.196^2, V4=0.427^2, KA=0, Q3=0.319^2, Q4=0.450^2, Favail=0, ALAG1=0.101^2),
notfixed_bpop=c(1,1,1,1,1,1,1,0,1),
xt=c(0, 7*24, 7*24+5/60, 7*24+20/60, 7*24+1, 7*24+12, 7*24+24),
minxt=c(-1, 7*24-1, 7*24+0, 7*24+10/60, 7*24+0.5, 7*24+10, 7*24+20),
maxxt=c(0, 7*24, 7*24+10/60, 7*24+0.5, 7*24+1.5, 7*24+14, 7*24+28),
a=c(DOSE=75,TAU=24),
mina=c(DOSE=37.5,TAU=24),
maxa=c(DOSE=150,TAU=24),
discrete_a = list(DOSE=seq(0,150,by=37.5),TAU=24))
#' plot intial design just PRED
plot_model_prediction(poped.db,model_num_points = 500)
I tried running the sample optimization R-script on Windows and took around 3.5 h, so decided to check what happens if I run the same script on Mac.
Does anyone have any experience with that?
@andrewhooker: Will it be possible for you too share your thoughts and how I can avoid such discrepancy? Thanks in advance
Ankur Sharma
Hello Andrew,
I am looking for a best way to include different covariates (e.g. weight, race, etc.). I have seen that one can add such in the "create_design" way but this seems to be good for some specifications per group. I would like to know if one can properly implement, for instance, an specific weight per subject for CL.
Thank you very much for your help in advance.
Roberto
Hi Andrew,
I am having trouble implementing a short IV infusion. The rationale for wanting to do this being that an IV ‘bolus’ is actually given over a period of time e.g. 3 minutes.
Using the code you posted on 11 Jul 2017, with low values of TINF = 0.5, the model prediction graph seems to miss some of the subsequent doses.
If TINF = 1, there doesn’t seem to be this problem.
The issue seems to depend on the DOSE/TINF value (and therefore the amount “In”).
Do you have any suggestions getting this short IV infusion to work?
Best,
Conor
Originally posted by @cohanlon2 in #11 (comment)
Hi Andrew,
I am wondering whether there is a way to implement the covariance of the sigma matrix, something like 'covd' in the function of 'create.poped.database'.
Could you please explain what "notfixed_" means in the function of 'create.poped.database'? Does it mean that those fixed parameters are not taken into account when evaluate/ optimize the sampling design?
Another question, is there any way to evaluate/optimize the sampling design to get good estimates of steady-state Cmax and AUC (in a direct way)? Currently, we optimize the design to have a good estimation of PK parameters (e.g. CL and V). Do you have any suggestion on this?
Thanks a lot.
Hello,
I recently issued a pull request with minor updates to the code written to enable parallelisation of poped_optim with mrgsolve models, and the automatic tests ran by appveyor failed.
After a bit of investigation it is due to the test called "ed_laplace_ofv works", where the following command fails when ran by AppVeyor :
expect_warning(expect_error(poped_optim(poped_db,opt_xt=T,parallel = T, d_switch=0,use_laplace=T)))
However when I download the archive with the tests and ran them locally, the code runs correctly (meaning that I get an error and a warning by the poped_optim function).
I do not know enough about these automated tests to propose a fix, but thought I would let you know anyways !
Dear,
I am using the PopED Matlab Version on Windows 10 / Matlab 2015a. Even running in admin and or compability mode doesnt prevent the app from crashing when starting the example files. I can open them in the PopED GUI but during the run I get an error saying "Matlab closed unexpected" even if it is already open and is opened again by PopED. Have you ever experienced this issue? What could be a solution to it?
Thanks,
Moritz
Hello @andrewhooker,
Thanks for this great package. I had a few questions about evaluate_fim_map()
.
First, the example section of the evaluate_fim_map()
documentation ends with shrinkage(poped.db)
. Shouldn't it be evaluate_fim_map(poped.db)
?
Line 83 in 70320a7
groupsize
is set to 32
for the poped.db
. I understand that it is a default/example value. But, can you confirm that the groupsize is not of concern for evaluate_fim_map()
since we work on a single-individual problem, and any value could be set?
Thanks in advance
Félicien
Is there a place to provide a gradient to speed up the calculation?
Thanks in advance,
Matt
A fix before I fix the code. If you first optimize:
output <- poped_optimize(poped.db,opt_xt=T)
afterwards, run the following lines of code:
output$poped.db$design$xt <- output$poped.db$gxt
output$poped.db$design$x <-output$poped.db$gx
output$poped.db$design$a <-output$poped.db$ga
Then the optimised design is displayed with
plot_model_prediction(output$poped.db)
Hi, I have an issue when trying to set-up a design with 2 groups and multiple responses (3).
sampling.times.full <- c(0,671.5,2015.5,2024,2028,2032,2040,2064,2088,2112,2184,4416,6432,8732) # full sampling group
N.full <- Nfull# number of subjects in the full sampling group
sampling.times.sparse <- c(0,671.5,2064,4416,6432,8732) # sparse sampling group
N.sparse <- 100-Nfull # number of subjects in the sparse groupcreate_design(
- groupsize=c(N.full,N.sparse),
- m=2, #number of groups
- xt=list(rep(sampling.times.full,3),
rep(sampling.times.sparse,3)),
- model_switch = list(sort(rep(c(1,2,3),length(sampling.times.full))),
sort(rep(c(1,2,3),length(sampling.times.sparse))))
- )
Error in test_mat_size(size(xt), model_switch, "model_switch") :
model_switch has dimensions 2x2 and should be 2x42
If for the purpose of testing I remove 'model_switch' (i.e. default), and run create.poped.database I get the following error:
Error in if (mat[i, j] == 0) { : missing value where TRUE/FALSE needed
What am I doing wrong?
Kind regards, Andreas.
Hi Andrew (@andrewhooker),
I am trying to do optimization for PK time-points using a complex dosing regime such as an IV loading dose and then sc maintenance doses once in 28 days x 10 doses. I tried looking in available resources for popED for such complex dosing regimens but didn't got any luck.
Also, how can I assign different combinations of doses (IV+SC) (e.g., in two such combinations) while optimizing time points using different combinations? There are examples of designs for different doses but I wish to learn how to assign different combinations of doses (IV+SC) under #design.
Will it be possible for you to share any codes for assigning the complex dosing regimens in popED. Highly appreciate it!
Thanks and regards,
-Ankur
Hi @andrewhooker,
I have two questions:
sigma
is fixed to one and the additive error is defined in the ff
model. When running evaluate_design
the residual components are not being estimated in the fim
:$fim
tcl tv add.err
tcl 30.2591311 -0.2562828 0
tv -0.2562828 40.0014869 0
add.err 0.0000000 0.0000000 0
$rse
tcl tv add.err
3.76520 30.95288 NA
Warning message:
The following parameters are not estimable:
add.err
Is the design adequate to estimate all parameters?
Is fixing the add.err
component equivalent to what is done?
f
. In dynamic transformation of both sides in NONMEM there is an objective function adjustment. Is that needed?Thanks
The answers to these two questions will help me produce an interface between nlmixr2
and PopED
.
Hi Andrwe,
I am running the wafarin example scripts in the R package.
My Initial parameters(sampling points) are: 0.5 1 2 6 24 36 72 120
The Optimized Sampling Schedule are:
1e-05 1e-05 34.29 34.29 75.92 120 120 120
Since the times of some sampling points are identical. Dose that mean that the sampling points can be reduced from 6 to 3 after removing redundant points ?
Thank you,
Sibo
Hi Andrew,
Is it possible to introduce irregular dosing when creating the PopED database? E.g. switching from a loading dose to a maintenance dose, or having some irregular dosing intervals within the same group/arm.
Thanks a lot in advance.
Ahmed Suleiman
Hello,
I am trying to evaluate and optimize a design using quasi equilibrium TMDD model in PopED. While evaluating the design, I get the warning message "## Problems inverting the matrix. Results could be misleading.". I tried a few different things, such as playing with tolerance, removing and fixing some fixed and random effects via notfixed_bpop and notfixed_d options, but the warning does not go away.
Can someone please share any tips on what may work to remove this warning? What are the implications if we end up moving forward with ignoring the warning message?
Many thanks in advance,
Krina
In the function reference, those two functions both serve for "Optimization main module for PopED". Is there any distinction between them? Which one is typically used for the purpose of optimization?
Thanks.
Hello,
Thank you for your work on PopED.
As with issue #40, I am trying to optimize the number of samples using a function that is not yet available in the R version of PopED.
Before this function gets implemented, I am wondering whether it would be a correct approach to run multiple optimization procedures while testing different designs (e.g. with 4 to 20 samples per subject) and judging the output?
If this is possible, what is the best approach to compare the different designs since more samples will always provide more information? Can this comparison be done by comparing the OFV? Plotting the mean RSE of the parameters over the number of samples? Determining a 'significant' change in the efficiency?
Thank you for your suggestions,
Michiel van Esdonk
The "Examples" vignette mentions that handling derived outputs is "Available in PopED, but not shown in examples". Where can I find more information about this? Can you provide an example showing how to handle derived outputs? Thank you for this amazing package.
hi Andy,
I am not sure how to diagnose the condition where the evaluate.fim
returns a 0
determinant. I know that my ff_file
returns a valid y
vector (checked by running the setup outside a function setting). The poped.db
seems to be have been set up correctly as far as the design goes, this I can confirm by plotting the predictions using the plot_model_prediction
which returns a population prediction at the specified design times. But when I run an evaluate.fim
on this poped.db
the return value is a matrix of zero values and hence the determinant of zero.
I want to go about diagnosing why this happens, and was wondering if you could provide some tips. Here is my create.poped.db
function. Just a note that A
and B
below are parameters that KA is derived from.
poped.db <- create.poped.database(ff_file="ff",
fg_file="sfg",
fError_file="feps",
bpop=c(CL=10,VC=60.1,Q=2,VP=30,TVLAG=0.2,A=0.256,B=0.284),
notfixed_bpop=c(1,1,1,1,0,0,0),
d=c(CL=0.08,VC=0.1),
sigma=c(0.06,6.48),
notfixed_sigma=c(0,0),
m=1,
groupsize=6,
xt=c(3,5,8,10),
minxt=c(3,4,7,9),
maxxt=c(4,7,9,10),
bUseGrouped_xt=0,
a = c(40000,70),
maxa=c(40000,70),
mina=c(40000,70))
Dear @andrewhooker,
Could you please help me navigate through the error message I keep getting when I try optimizing sampling scheme:
Total iterations: 336
Elapsed time: 3794.93 seconds.
Final OFV = 308.99
Parameters: 1103.538 4.108349 1431.097 1025.364 3600
Running BFGS Optimization
initial value -282.706221
Problems inverting the matrix. Results could be misleading.
Error in svd(mat) : infinite or missing values in 'x'
With kind regards,
Ahmed
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.