Sparse Sampling for Population Models using PopED
Heller
2025-02-28
#libraries
rm(list=ls())
#install.packages ("PopED")
library(PopED)## Warning: package 'PopED' was built under R version 4.4.2#packageVersion("PopED")##Define the Model #one-compartment pharmacokinetic model with linear absorption
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
N = floor(xt/TAU)+1
f=(DOSE*Favail/V)*(KA/(KA - CL/V)) *
(exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
return(list( f=f,poped.db=poped.db))
})
}#Define Parameters
sfg <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU=a[2])
}#Define RUV
feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- ff(model_switch,xt,parameters,poped.db)
f <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = f*(1+epsi[,1])+epsi[,2]
return(list(y=y,poped.db=poped.db))
}#Create poped database
poped.db <- create.poped.database(
# Model
ff_fun=ff,
fg_fun=sfg,
fError_fun=feps,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(prop=0.04,add=5e-6),
notfixed_sigma=c(1,0),
# Design
m=2,
groupsize=20,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24),
xt=c( 1,2,8,240,245),
# Design space
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
bUseGrouped_xt=TRUE)#Plot model prediction
plot_model_prediction(poped.db, model_num_points = 300)
#plot the expected prediction interval
plot_model_prediction(poped.db,
PI=TRUE,
separate.groups=T,
model_num_points = 300,
sample.times = FALSE)
#predictions numerically
dat <- model_prediction(poped.db,DV=TRUE)
head(dat,n=5);tail(dat,n=5)## ID Time DV IPRED PRED Group Model DOSE TAU
## 1 1 1 0.03021539 0.03065119 0.05325024 1 1 20 24
## 2 1 2 0.04233682 0.05439378 0.09204804 1 1 20 24
## 3 1 8 0.09548180 0.11126030 0.16409609 1 1 20 24
## 4 1 240 0.10143200 0.12658722 0.12671376 1 1 20 24
## 5 1 245 0.22143036 0.20042526 0.24980320 1 1 20 24## ID Time DV IPRED PRED Group Model DOSE TAU
## 196 40 1 0.1201239 0.1013835 0.1065005 2 1 40 24
## 197 40 2 0.1930922 0.1750751 0.1840961 2 1 40 24
## 198 40 8 0.2605835 0.3001336 0.3281922 2 1 40 24
## 199 40 240 0.1371124 0.1428533 0.2534275 2 1 40 24
## 200 40 245 0.4112760 0.3820690 0.4996064 2 1 40 24#evaluate the initial design
(ds1 <- evaluate_design(poped.db))## $ofv
## [1] 39.309
##
## $fim
## V KA CL d_V d_KA d_CL
## V 0.05336692 -8.683963 -0.05863412 0.000000 0.000000 0.000000
## KA -8.68396266 2999.851007 -14.43058560 0.000000 0.000000 0.000000
## CL -0.05863412 -14.430586 37.15243290 0.000000 0.000000 0.000000
## d_V 0.00000000 0.000000 0.00000000 999.953587 312.240246 3.202847
## d_KA 0.00000000 0.000000 0.00000000 312.240246 439.412556 2.287838
## d_CL 0.00000000 0.000000 0.00000000 3.202847 2.287838 3412.005199
## sig_prop 0.00000000 0.000000 0.00000000 575.347261 638.581909 1182.325475
## sig_prop
## V 0.0000
## KA 0.0000
## CL 0.0000
## d_V 575.3473
## d_KA 638.5819
## d_CL 1182.3255
## sig_prop 33864.3226
##
## $rse
## V KA CL d_V d_KA d_CL sig_prop
## 8.215338 10.090955 4.400304 39.844763 60.655110 27.562541 13.865357#Try 5 samples with four early and one late - model “db_new1”
poped.db_new1 <- create.poped.database(
# Model
ff_fun=ff,
fg_fun=sfg,
fError_fun=feps,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(prop=0.04,add=5e-6),
notfixed_sigma=c(1,0),
# Design
m=2,
groupsize=20,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24),
xt=c( 1,2,8,10,248),
# Design space
minxt=c(0,0,0,0,240),
maxxt=c(10,10,10,10,248),
bUseGrouped_xt=TRUE)#plot the expected prediction interval for db_new1
plot_model_prediction(poped.db_new1,
PI=TRUE,
separate.groups=T,
model_num_points = 300,
sample.times = FALSE)
#evaluate the design db_new1
(ds2 <- evaluate_design(poped.db_new1))## $ofv
## [1] 38.6089
##
## $fim
## V KA CL d_V d_KA d_CL
## V 0.05603522 -7.708389 0.2139121 0.00000 0.00000 0.00000
## KA -7.70838947 3247.802593 32.0070606 0.00000 0.00000 0.00000
## CL 0.21391207 32.007061 23.8488710 0.00000 0.00000 0.00000
## d_V 0.00000000 0.000000 0.0000000 1102.44707 246.02562 42.62909
## d_KA 0.00000000 0.000000 0.0000000 246.02562 515.05406 11.25545
## d_CL 0.00000000 0.000000 0.0000000 42.62909 11.25545 1405.95318
## sig_prop 0.00000000 0.000000 0.0000000 608.89790 795.37576 1870.41520
## sig_prop
## V 0.0000
## KA 0.0000
## CL 0.0000
## d_V 608.8979
## d_KA 795.3758
## d_CL 1870.4152
## sig_prop 35181.7085
##
## $rse
## V KA CL d_V d_KA d_CL sig_prop
## 7.431050 8.892071 5.777082 35.440487 52.564091 44.308695 14.090858##Try 4 samples with three early and one late - model “db_new2”
poped.db_new2 <- create.poped.database(
# Model
ff_fun=ff,
fg_fun=sfg,
fError_fun=feps,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(prop=0.04,add=5e-6),
notfixed_sigma=c(1,0),
# Design
m=2,
groupsize=20,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24),
xt=c( 1,2,8,248),
# Design space
minxt=c(0,0,0,240),
maxxt=c(10,10,10,248),
bUseGrouped_xt=TRUE)#plot the expected prediction interval for db_new2
plot_model_prediction(poped.db_new2,
PI=TRUE,
separate.groups=T,
model_num_points = 300,
sample.times = FALSE)
#evaluate the design db_new2
(ds3 <- evaluate_design(poped.db_new2))## $ofv
## [1] 37.02264
##
## $fim
## V KA CL d_V d_KA d_CL
## V 0.05093387 -8.914773 0.1572673 0.0000 0.00000 0.00000
## KA -8.91477347 2962.509660 18.6114517 0.0000 0.00000 0.00000
## CL 0.15726732 18.611452 23.2198932 0.0000 0.00000 0.00000
## d_V 0.00000000 0.000000 0.0000000 910.8544 329.05905 23.04160
## d_KA 0.00000000 0.000000 0.0000000 329.0590 428.54106 3.80578
## d_CL 0.00000000 0.000000 0.0000000 23.0416 3.80578 1332.77128
## sig_prop 0.00000000 0.000000 0.0000000 545.2588 594.11904 1935.44440
## sig_prop
## V 0.0000
## KA 0.0000
## CL 0.0000
## d_V 545.2588
## d_KA 594.1190
## d_CL 1935.4444
## sig_prop 24679.9751
##
## $rse
## V KA CL d_V d_KA d_CL sig_prop
## 9.231920 11.057673 5.789418 43.325191 63.926813 46.650646 17.228914#Let’s try to optimize ds3. This is design db_new2
output <- poped_optim(poped.db_new2, opt_xt=TRUE)## ===============================================================================
## Initial design evaluation
##
## Initial OFV = 37.0226
##
## Initial design
## expected relative standard error
## (%RSE, rounded to nearest integer)
## Parameter Values RSE_0
## V 72.8 9
## KA 0.25 11
## CL 3.75 6
## d_V 0.09 43
## d_KA 0.09 64
## d_CL 0.0625 47
## sig_prop 0.04 17
##
## ==============================================================================
## Optimization of design parameters
##
## * Optimize Sampling Schedule
##
## ************* Iteration 1 for all optimization methods***********************
##
## *******************************************
## Running Adaptive Random Search Optimization
## *******************************************
## Initial OFV = 37.0226
## It. 5 | OFV = 37.3105
## It. 10 | OFV = 37.3105
## It. 15 | OFV = 37.3105
## It. 20 | OFV = 37.776
## It. 25 | OFV = 37.776
## It. 30 | OFV = 37.776
## It. 35 | OFV = 37.776
## It. 40 | OFV = 37.776
## It. 45 | OFV = 37.776
## It. 50 | OFV = 37.8718
## It. 55 | OFV = 37.8718
## It. 60 | OFV = 37.8718
## It. 65 | OFV = 37.8718
## It. 70 | OFV = 37.8718
## It. 75 | OFV = 37.8718
## It. 80 | OFV = 37.8718
## It. 85 | OFV = 37.8718
## It. 90 | OFV = 37.8718
## It. 95 | OFV = 37.8718
## It. 100 | OFV = 37.9063
## It. 105 | OFV = 37.9063
## It. 110 | OFV = 39.2877
## It. 115 | OFV = 39.2877
## It. 120 | OFV = 39.2877
## It. 125 | OFV = 39.2877
## It. 130 | OFV = 39.2877
## It. 135 | OFV = 39.2877
## It. 140 | OFV = 39.4671
## It. 145 | OFV = 39.4671
## It. 150 | OFV = 39.4671
## It. 155 | OFV = 39.6128
## It. 160 | OFV = 39.6128
## It. 165 | OFV = 39.6128
## It. 170 | OFV = 39.6128
## It. 175 | OFV = 39.6128
## It. 180 | OFV = 39.6128
## It. 185 | OFV = 39.6128
## It. 190 | OFV = 39.6128
## It. 195 | OFV = 39.6128
## It. 200 | OFV = 39.6128
## It. 205 | OFV = 39.6128
## It. 210 | OFV = 39.6128
## It. 215 | OFV = 39.6128
## It. 220 | OFV = 39.6128
## It. 225 | OFV = 39.6128
## It. 230 | OFV = 39.6128
## It. 235 | OFV = 39.6128
## It. 240 | OFV = 39.6128
## It. 245 | OFV = 39.6128
## It. 250 | OFV = 39.6128
## It. 255 | OFV = 39.6128
## It. 260 | OFV = 39.6128
## It. 265 | OFV = 39.6128
## It. 270 | OFV = 39.6128
## It. 275 | OFV = 39.6128
## It. 280 | OFV = 39.6128
## It. 285 | OFV = 39.6128
## It. 290 | OFV = 39.6128
## It. 295 | OFV = 39.6128
## It. 300 | OFV = 39.6128
## It. 305 | OFV = 39.6128
## It. 310 | OFV = 39.6128
## It. 315 | OFV = 39.6128
## It. 320 | OFV = 39.6128
## It. 325 | OFV = 39.6128
## It. 330 | OFV = 39.6128
## It. 335 | OFV = 39.6128
## It. 340 | OFV = 39.6128
## It. 345 | OFV = 39.6128
## It. 350 | OFV = 39.6128
## Maximum number of identical optimal values reached (max_run=200), optimization stopped.
##
## Total iterations: 354
## Elapsed time: 6.9 seconds.
##
## Final OFV = 39.61277
## Parameters: 0.3926222 10 10 240
##
## *******************************************
## Running BFGS Optimization
## *******************************************
## initial value -39.612772
## final value -39.621173
## converged
##
## *******************************************
## Running Line Search Optimization
## *******************************************
##
## Initial parameters: 0.4648692 10 10 240
## Initial OFV: 39.62117
##
## Searching parameter 2
## No change
## Searching parameter 1
## No change
## Searching parameter 4
## Changed from 240 to 240 ; OFV = 39.6212
## Searching parameter 3
## No change
##
## Elapsed time: 3.59 seconds.
##
## Final OFV = 39.62117
## Parameters: 0.4648692 10 10 240
##
## *******************************************
## Stopping criteria testing
## (Compare between start of iteration and end of iteration)
## *******************************************
## Difference in OFV: 2.6
## Relative difference in OFV: 7.02%
## Efficiency:
## ((exp(ofv_final) / exp(ofv_init))^(1/n_parameters)) = 1.4495
##
## Efficiency stopping criteria:
## Is (1.4495 <= 1.001)? No.
## Stopping criteria NOT achieved.
##
## Stopping criteria NOT achieved.
##
## ************* Iteration 2 for all optimization methods***********************
##
## *******************************************
## Running Adaptive Random Search Optimization
## *******************************************
## Initial OFV = 39.6212
## It. 5 | OFV = 39.6212
## It. 10 | OFV = 39.6212
## It. 15 | OFV = 39.6212
## It. 20 | OFV = 39.6212
## It. 25 | OFV = 39.6212
## It. 30 | OFV = 39.6212
## It. 35 | OFV = 39.6212
## It. 40 | OFV = 39.6212
## It. 45 | OFV = 39.6212
## It. 50 | OFV = 39.6212
## It. 55 | OFV = 39.6212
## It. 60 | OFV = 39.6212
## It. 65 | OFV = 39.6212
## It. 70 | OFV = 39.6212
## It. 75 | OFV = 39.6212
## It. 80 | OFV = 39.6212
## It. 85 | OFV = 39.6212
## It. 90 | OFV = 39.6212
## It. 95 | OFV = 39.6212
## It. 100 | OFV = 39.6212
## It. 105 | OFV = 39.6212
## It. 110 | OFV = 39.6212
## It. 115 | OFV = 39.6212
## It. 120 | OFV = 39.6212
## It. 125 | OFV = 39.6212
## It. 130 | OFV = 39.6212
## It. 135 | OFV = 39.6212
## It. 140 | OFV = 39.6212
## It. 145 | OFV = 39.6212
## It. 150 | OFV = 39.6212
## It. 155 | OFV = 39.6212
## It. 160 | OFV = 39.6212
## It. 165 | OFV = 39.6212
## It. 170 | OFV = 39.6212
## It. 175 | OFV = 39.6212
## It. 180 | OFV = 39.6212
## It. 185 | OFV = 39.6212
## It. 190 | OFV = 39.6212
## It. 195 | OFV = 39.6212
## It. 200 | OFV = 39.6212
## Maximum number of identical optimal values reached (max_run=200), optimization stopped.
##
## Total iterations: 200
## Elapsed time: 3.75 seconds.
##
## Final OFV = 39.62117
## Parameters: 0.4648692 10 10 240
##
## *******************************************
## Running BFGS Optimization
## *******************************************
## initial value -39.621173
## final value -39.621173
## converged
##
## *******************************************
## Running Line Search Optimization
## *******************************************
##
## Initial parameters: 0.4648692 10 10 240
## Initial OFV: 39.62117
##
## Searching parameter 4
## Changed from 240 to 240 ; OFV = 39.6212
## Searching parameter 2
## Changed from 10 to 10 ; OFV = 39.6212
## Searching parameter 1
## No change
## Searching parameter 3
## Changed from 10 to 10 ; OFV = 39.6212
##
## Elapsed time: 3.73 seconds.
##
## Final OFV = 39.62117
## Parameters: 0.4648692 10 10 240
##
## *******************************************
## Stopping criteria testing
## (Compare between start of iteration and end of iteration)
## *******************************************
## Difference in OFV: -1.3e-09
## Relative difference in OFV: -3.29e-09%
## Efficiency:
## ((exp(ofv_final) / exp(ofv_init))^(1/n_parameters)) = 1
##
## Efficiency stopping criteria:
## Is (1 <= 1.001)? Yes.
## Stopping criteria achieved.
##
## Stopping criteria achieved.
##
## ===============================================================================
## FINAL RESULTS
## Optimized Sampling Schedule
## Group 1: 0.4649 10 10 240
## Group 2: 0.4649 10 10 240
##
## OFV = 39.6212
##
## Efficiency:
## ((exp(ofv_final) / exp(ofv_init))^(1/n_parameters)) = 1.4495
##
## Expected relative standard error
## (%RSE, rounded to nearest integer):
## Parameter Values RSE_0 RSE
## V 72.8 9 6
## KA 0.25 11 8
## CL 3.75 6 4
## d_V 0.09 43 33
## d_KA 0.09 64 53
## d_CL 0.0625 47 29
## sig_prop 0.04 17 19
##
## Total running time: 19.72 seconds#Let’s do some efficiency comparisons #Compare efficiency of ds2 vs. ds1; this is 4 early, 1 late vs. original 3 early, 2 late
efficiency(ds2$ofv,ds1$ofv,poped.db)## [1] 1.105187
## attr(,"description")
## [1] "(exp(ofv_final) / exp(ofv_init))^(1/n_parameters)"#My ds2 design requires ~ 10% more subjects to get the same info as the original design
#Compare efficiency of ds3 vs. ds1; this is 3 early, 1 late vs. original 3 early, 2 late
efficiency(ds3$ofv,ds1$ofv,poped.db)## [1] 1.386279
## attr(,"description")
## [1] "(exp(ofv_final) / exp(ofv_init))^(1/n_parameters)"#my ds3 design (4 samples) requires about 40% more subjects to get the same info as the original design.
#Compare efficiency of ds3 vs. ds2; This is 3 early, 1 late vs. 4 early, 1 late
efficiency(ds3$ofv,ds2$ofv,poped.db)## [1] 1.254339
## attr(,"description")
## [1] "(exp(ofv_final) / exp(ofv_init))^(1/n_parameters)"#My ds3 design (four samples) requires about 25% more subjects to get the same information as myds2 design (five samples).