Open Anks2030 opened 3 years ago
andrewhooker
commented
3 years ago Hi Ankur,
This can certainly be implemented! You just need to define different groups with different designs in the trial. See for example https://andrewhooker.github.io/PopED/articles/intro-poped.html. You can post example code here if you like and I can try and help.
Best regards, Andy
Anks2030
commented
3 years ago Hi Andrew, Thanks for looking into this. In our trial, we have N=200 in treatment arm until wk-52. Post-52 week we have PK at D1, D2, D3, D4 and D7, D14. And post-52WK the collection is for only N=20. I dont know how to implement this switch of N=200 to N=20 and have an overall influence of design? So far I am just assuming 200 for entire duration and have even tried to evaluate the design only till WK-52 (using just N=200) Please have code for entire study for N=200
PK.2.comp.sc.ode <- function(Time, State, Pars){ with(as.list(c(State, Pars)), {
CL=CL*(WT/76.8)^(WT_CL)
V1=V1*(WT/76.8)^(WT_V1)
dA1 <- -KA*A1
dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
dA3 <- A2* Q/V1-A3* Q/V2
return(list(c(dA1, dA2, dA3)))}) }
ff.PK.2.comp.sc.md.ode <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{ A_ini <- c(A1=0, A2=0, A3=0) times_xt <- drop(xt) dose_times = c(seq(from=0,to=336, by=84),seq(from=336,to=8568,by=TAU)) eventdat <- data.frame(var = c("A1"), time = dose_times, value = c(DOSE), method = c("add")) times <- sort(c(times_xt,dose_times)) times <- unique(times) # remove duplicates out <- ode(A_ini, times, PK.2.comp.sc.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13) y = (out[, "A2"]/(V1/Favail))+BL y=y[match(times_xt,out[,"time"])] y=cbind(y) return(list(y=y,poped.db=poped.db)) }) }
fg <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]exp(b[1]), V1=bpop[2]exp(b[2]), KA=bpop[3], Q= bpop[4], V2=bpop[5], Favail=bpop[6], WT_CL=bpop[7], WT_V1=bpop[8], BL=bpop[9], DOSE=a[1], TAU=a[2], WT= a[3]) return( parameters ) }
feps <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]]
y = y*(1+epsi[,1])+epsi[,2]
return(list( y= y,poped.db =poped.db )) }
poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.sc.md.ode", fg_fun="fg", fError_file="feps",
bpop=c(CL=0.00633,V1=8.75,KA=0.02108,Q= 0.0344, V2= 1.8, Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
notfixed_bpop=c(1,1,1,0,1,0,0,0,0),# decides which tested and which are not
d=c(CL=0.282,V1=0.569), # decides the variances of BSV
notfixed_d = c(1,1),
sigma=c(0.00662,0), # decides the variances of RV prop, additive
notfixed_sigma=c(1,0),
m=1, #number of groups
groupsize=200,
#0,7,28,35,56,59,63,70,77,84
xt= c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904),
minxt=c(0,312,648,1992,2664,3336,6024,7368,8544,8568,8592,8616,8640,8712,8880),
maxxt=c(0,360,696,2040,2712,3384,6072,7416,8592,8616,8640,8664,8688,8760,8928),
bUseGrouped_xt=1,
a=c(DOSE=0.15*76.8,TAU=168, WT=76.8))dat_Original <- model_prediction(poped.db) plot_model_prediction(poped.db,DV=T,
sample.times=F,alpha.sample.times.DV.points = 0.05,
model_num_points = 1000)+# , PI=Tlabs(x = "Time from first dose (Weeks)", y="IgG Conc.(g/L)")+ theme_bw()+ geom_vline(xintercept=336,lty="dashed")+ geom_point(data = dat_Original, aes(Time, PRED ), size=2, alpha=1/2, colour="blue")+#colour="Before Optimization" ggtitle("Original Design")+coord_cartesian(ylim=c(4,12))+ scale_y_continuous(breaks=sort(c(seq(0,12,1))))+ scale_x_continuous(breaks = seq(0,8928,168), labels=c(seq(1,54,1)))+ theme_bw()
(dso <- evaluate_design(poped.db))
andrewhooker
commented
3 years ago Hi Ankur,
Here is am implementation of what you are trying to do.
First define the model:
library(tidyverse)
library(deSolve)
library(PopED)
packageVersion("PopED") #0.5.0
#' Define the ODE system
PK.2.comp.sc.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
CL=CL*(WT/76.8)^(WT_CL)
V1=V1*(WT/76.8)^(WT_V1)
dA1 <- -KA*A1
dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
dA3 <- A2* Q/V1-A3* Q/V2
return(list(c(dA1, dA2, dA3)))
})
}
#' define the initial conditions and the dosing
ff.PK.2.comp.sc.md.ode <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0, A3=0)
times_xt <- drop(xt)
dose_times = c(seq(from=0,to=336, by=84),seq(from=336,to=8568,by=TAU))
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE), method = c("add"))
times <- sort(c(times_xt,dose_times))
times <- unique(times) # remove duplicates
out <- ode(A_ini, times, PK.2.comp.sc.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = (out[, "A2"]/(V1/Favail))+BL
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]),
V1=bpop[2]*exp(b[2]),
KA=bpop[3],
Q= bpop[4],
V2=bpop[5],
Favail=bpop[6],
WT_CL=bpop[7],
WT_V1=bpop[8],
BL=bpop[9],
DOSE=a[1],
TAU=a[2],
WT= a[3])
return( parameters )
}
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional
feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = y*(1+epsi[,1])+epsi[,2]
return(list( y= y,poped.db =poped.db ))
}Add a design with design with extra observations after 52 wks in 20 patients.
poped_db_1 <- create.poped.database(
ff_fun="ff.PK.2.comp.sc.md.ode",
fg_fun="fg",
fError_file="feps",
bpop=c(CL=0.00633,V1=8.75,KA=0.02108,
Q= 0.0344,
V2= 1.8,
Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
notfixed_bpop=c(1,1,1,
0,
1,
0,0,0,0),# decides which tested and which are not
d=c(CL=0.282,V1=0.569), # decides the variances of BSV
notfixed_d = c(1,1),
sigma=c(0.00662,0), # decides the variances of RV prop, additive
notfixed_sigma=c(1,0),
m=2, #number of groups
groupsize=c(20,180),
#0,7,28,35,56,59,63,70,77,84
# Here we have 20 individuals that have sampling like the rest until wk 52 and
# then more samples post wk52.
# 180 have only sampling until wk 52.
xt= list(c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904),
c(0,336,671,2016,2688,3360,6048,7392,8568)),
minxt=list(c(0,312,648,1992,2664,3336,6024,7368,8544,8568,8592,8616,8640,8712,8880),
c(0,312,648,1992,2664,3336,6024,7368,8544)),
maxxt=list(c(0,360,696,2040,2712,3384,6072,7416,8592,8616,8640,8664,8688,8760,8928),
c(0,360,696,2040,2712,3384,6072,7416,8592)),
#bUseGrouped_xt=1,
a=list(c(DOSE=0.15*76.8,TAU=168, WT=76.8),
c(DOSE=0.15*76.8,TAU=168, WT=76.8))
)Evaluate design and plot design.
(p1 <- plot_model_prediction(poped_db_1,separate.groups = T,model_num_points = 1000))
(des1 <- evaluate_design(poped_db_1))Create a design without the extra observations after 52 weeks:
poped_db_2 <- create.poped.database(
ff_fun="ff.PK.2.comp.sc.md.ode",
fg_fun="fg",
fError_file="feps",
bpop=c(CL=0.00633,V1=8.75,KA=0.02108,
Q= 0.0344,
V2= 1.8,
Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
notfixed_bpop=c(1,1,1,
0,
1,
0,0,0,0),# decides which tested and which are not
d=c(CL=0.282,V1=0.569), # decides the variances of BSV
notfixed_d = c(1,1),
sigma=c(0.00662,0), # decides the variances of RV prop, additive
notfixed_sigma=c(1,0),
m=1, #number of groups
groupsize=c(200),
#0,7,28,35,56,59,63,70,77,84
# Here we have 20 individuals that have sampling like the rest until wk 52 and
# then more samples post wk52.
# 180 have only sampling until wk 52.
xt= c(0,336,671,2016,2688,3360,6048,7392,8568),
minxt=c(0,312,648,1992,2664,3336,6024,7368,8544),
maxxt=c(0,360,696,2040,2712,3384,6072,7416,8592),
#bUseGrouped_xt=1,
a=c(DOSE=0.15*76.8,TAU=168, WT=76.8)
)
Evaluate design and plot design.
(p2 <- plot_model_prediction(poped_db_2,separate.groups = T,model_num_points = 1000))
(des2 <- evaluate_design(poped_db_2))There is a big difference in parameter RSE between the two designs
tibble::tibble(names(des1$rse),des1$rse,des2$rse)There is also a big difference in efficiency (52% more individuals in the design without extra observations after 52 wks are needed to achieve same information as the extended design)
efficiency(ofv_init = des2$ofv, ofv_final = des1$ofv, poped_db = poped_db_1)Hi Andy, Thanks for taking time to go over this challenge I have a question. I am trying hard to understand if the code is capturing the design appropriately
I am attaching a single design figure (total duration is 1year) in which we have N=180 only till WK-52 and after WK-52 we have ONLY N=20
I trying to understand if this design is applied in below code you shared?
m=2, #number of groups groupsize=c(20,180),
xt= list(c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904), c(0,336,671,2016,2688,3360,6048,7392,8568)),
The output is Group 1 has N=20 only till WK54 and Group 2 is N=180 till WK-52 I wish to see one group which has overall impact of design (N=180 patients only till WK52 and after WK52 I have only 20 subjects for rich PK)
I am not sure if I am interpreting it correctly please advice
Really appreciate you spreading the word for PopED. Very useful tool. Thanks, -Ankur
andrewhooker
commented
3 years ago Hi
The design I implemented has 180 individuals with measurements between 0 and 52 weeks, and 20 individuals with those measurements PLUS additional measurements after 52 weeks. I assumed that the 20 individuals would be a part of the 200 total that are studies under 52 weeks and then 20 of those individuals are studied for an additional period.
Anks2030
commented
3 years ago Hi Andy, Thank you very much it is clear now. Very useful and I really appreciate for your time and effort. Do you have the analytical solution for the 2-compartment model. I will try to do the optimization for this design and seems it is taking forever for ODE's Thanks, -Ankur
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