diff --git a/R/poped.R b/R/poped.R index f032ad66..87beb62f 100644 --- a/R/poped.R +++ b/R/poped.R @@ -1252,6 +1252,8 @@ attr(rxUiGet.popedNotfixedSigma, "desc") <- "PopED database $notfixed_sigma" .env$mt <- -Inf .wcmt <- which(.nd == "cmt") .wdvid <- which(.nd == "dvid") + .wg_xt <- which(.nd == "g_xt") + .G_xt <- NULL .multipleEndpoint <- FALSE .poped$uid <- unique(.data[[.wid]]) if (length(ui$predDf$cond) > 1L) { @@ -1267,9 +1269,17 @@ attr(rxUiGet.popedNotfixedSigma, "desc") <- "PopED database $notfixed_sigma" .env$mt <- max(c(.time, .env$mt)) if (length(.wdvid) == 1L) { .wd <- which(.data[[.wdvid]] == i) - if (length(.wd) == 0) .wd <- which(.data[[.wdvid]] == ui$predDf$cond[i]) + if (length(.wd) == 0) { + .wd <- which(.data[[.wdvid]] == + ui$predDf$cond[i]) + } if (length(.wd) > 0) { .time <- .time[.wd] + if (length(.wg_xt) == 1L) { + .g_xt <- .data[[.wg_xt]] + .g_xt <- .g_xt[.wd] + return(time=.time, dvid=i, G_xt=.g_xt) + } return(data.frame(time=.time, dvid=i)) } } @@ -1281,6 +1291,11 @@ attr(rxUiGet.popedNotfixedSigma, "desc") <- "PopED database $notfixed_sigma" } if (length(.wd) > 0) { .time <- .time[.wd] + if (length(.wg_xt) == 1L) { + .g_xt <- .data[[.wg_xt]] + .g_xt <- .g_xt[.wd] + return(time=.time, dvid=i, G_xt=.g_xt) + } return(data.frame(time=.time, dvid=i)) } } @@ -1288,6 +1303,9 @@ attr(rxUiGet.popedNotfixedSigma, "desc") <- "PopED database $notfixed_sigma" call.=FALSE) })) }) + if (length(.wg_xt) == 1L) { + .G_xt <- .xt$G_xt + } } else { .xt <- lapply(.poped$uid, function(id) { @@ -1373,6 +1391,7 @@ attr(rxUiGet.popedNotfixedSigma, "desc") <- "PopED database $notfixed_sigma" }) if (.single) .maxxt <- .maxxt[[1]] } + .poped$G_xt <- .G_xt .poped$discrete_a <- discrete_a .poped$discrete_xt <- discrete_xt .designSpace1 <- list(maxni = maxni, @@ -1822,7 +1841,8 @@ rxUiGet.popedSettings <- function(x, ...) { G_xt=.env$G_xt, a=.a, discrete_xt=.poped$discrete_xt, - discrete_a=.poped$discrete_a) + discrete_a=.poped$discrete_a, + G_xt=.poped$G_xt) return(.appendPopedProps(.ret, .ctl)) } else { .w <- which(.nd == "evid") @@ -1912,6 +1932,7 @@ rxUiGet.popedSettings <- function(x, ...) { " fError_fun=fepsFun,", paste0(" discrete_xt=", deparse1(.poped$discrete_xt), ","), paste0(" discrete_a=", deparse1(.poped$discrete_a), ","), + paste0(" G_xt=", deparse1(.poped$G_xt), ","), paste0(" bUseGrouped_xt=", deparse1(rxode2::rxGetControl(ui, "bUseGrouped_xt", FALSE)) , ", "), paste0(" m=", length(.a), ", #number of groups"), paste0(" groupsize=", .groupsize, ", #group size"), @@ -2034,7 +2055,8 @@ rxUiGet.popedSettings <- function(x, ...) { fError_fun=.err, bUseGrouped_xt=rxode2::rxGetControl(ui, "bUseGrouped_xt", FALSE), discrete_xt=.poped$discrete_xt, - discrete_a=.poped$discrete_a) + discrete_a=.poped$discrete_a, + G_xt=.poped$G_xt) } else { .ln <- tolower(names(data)) @@ -2101,6 +2123,7 @@ rxUiGet.popedSettings <- function(x, ...) { " fError_fun=fepsFun, ", paste0(" discrete_xt=", deparse1(.poped$discrete_xt), ","), paste0(" discrete_a=", deparse1(.poped$discrete_a), ","), + paste0(" G_xt=", deparse1(.poped$G_xt), ","), paste0(" bUseGrouped_xt=", deparse1(rxode2::rxGetControl(ui, "bUseGrouped_xt", FALSE)) ,")"), "", "# Plot the model", diff --git a/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R new file mode 100644 index 00000000..9e62f452 --- /dev/null +++ b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R @@ -0,0 +1,98 @@ + +library(babelmixr2) + +library(PopED) + +##-- Model: One comp first order absorption + inhibitory imax +## -- works for both mutiple and single dosing +f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tFavail <- fix(0.9) + tE0 <- 1120 + tImax <- 0.87 + tIC50 <- 0.0993 + + eta.v ~ 0.09 + eta.ka ~ 0.09 + eta.cl ~ 0.25^2 + eta.e0 ~ 0.09 + + conc.prop.sd <- fix(sqrt(0.04)) + conc.sd <- fix(sqrt(5e-6)) + + eff.prop.sd <- fix(sqrt(0.09)) + eff.sd <- fix(sqrt(100)) + }) + model({ + V<- tV*exp(eta.v) + KA <- tKa*exp(eta.ka) + CL <- tCl*exp(eta.cl) + Favail <- tFavail + E0 <- tE0*exp(eta.e0) + IMAX <- tImax + IC50 <- tIC50 + # PK + N <- floor(time/TAU)+1 + CONC <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (time - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + # PD model + EFF <- E0*(1 - CONC*IMAX/(IC50 + CONC)) + + CONC ~ add(conc.sd) + prop(conc.prop.sd) + EFF ~ add(eff.sd) + prop(eff.prop.sd) + + }) +} + +# Note that design point 240 is repeated +e1 <- et(c( 1,2,8,240, 240)) %>% + as.data.frame() %>% + dplyr::mutate(dvid=1) + +e1$low <- c(0,0,0,240, 240) +e1$high <- c(10,10,10,248, 248) +# Since the design point is repeated, there needs to be a grouping +# variable which is defined in the dataset as G_xt since it is defined +# in PopED as G_xt +e1$G_xt <- seq_along(e1$low) + +e2 <- e1 +e2$dvid <- 2 +e <- rbind(e1, e2) + +babel.db <- nlmixr2(f, e, "poped", + popedControl( + groupsize=20, + discrete_xt = list(0:248), + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24), + c(DOSE=40, TAU=24), + c(DOSE=0, TAU=24)), + maxa=c(DOSE=200,TAU=40), + mina=c(DOSE=0,TAU=2), + ourzero=0 + )) + + +## create plot of model and design +plot_model_prediction(babel.db,facet_scales="free") + +plot_model_prediction(babel.db,IPRED=T,DV=T,facet_scales="free",separate.groups=T) + +## evaluate initial design +evaluate_design(babel.db) +shrinkage(babel.db) + +# Optimization +output <- poped_optim(babel.db, opt_xt = T, parallel = T) + +summary(output) + +get_rse(output$FIM,output$poped.db) + +plot_model_prediction(output$poped.db,facet_scales="free") diff --git a/tests/testthat/test-poped.R b/tests/testthat/test-poped.R index 69ec1547..93130e40 100644 --- a/tests/testthat/test-poped.R +++ b/tests/testthat/test-poped.R @@ -295,23 +295,6 @@ if (requireNamespace("PopED", quietly=TRUE)) { sample.times = FALSE)) - dat <- model_prediction(db,DV=TRUE) - - expect_equal(head(dat, n=4), - data.frame( - ID = factor(c(1L, 1L, 1L, 1L), levels=paste(1:40)), - Time = c(1, 2, 8, 240), - DV = c(0.0353822273010824, 0.0721765325175048, - 0.142203020963518, 0.121570466918341), - IPRED = c(0.0379965748703227, 0.0654575999147953, - 0.118727861585151, 0.15387388187677), - PRED = c(0.0532502332862765, 0.0920480197661157, - 0.164096088998621, 0.126713764327394), - Group = factor(c(1L, 1L, 1L, 1L), levels = c("1", "2")), - Model = factor(c(1L, 1L, 1L, 1L), levels = "1"), - DOSE = c(20, 20, 20, 20)), - tolerance=1e-4) - expect_equal(evaluate_design(db), list(ofv = 39.3090057657525, fim = lotri({ @@ -334,111 +317,278 @@ if (requireNamespace("PopED", quietly=TRUE)) { }) - test_that("single endpoint solve", { + test_that("example 1", { + + library(PopED) + + f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + eta.v ~ 0.09 + eta.ka ~ 0.09 + eta.cl ~0.25^2 + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + Favail <- tF + N <- floor(time/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (time - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + y ~ prop(prop.sd) + add(add.sd) + }) + } + + # minxt, maxxt + e <- et(list(c(0, 10), + c(0, 10), + c(0, 10), + c(240, 248), + c(240, 248))) %>% + as.data.frame() + + #xt + e$time <- c(1,2,8,240,245) + - ff_analytic <- function(model_switch,xt,parameters,poped.db){ + babel.db <- nlmixr2(f, e, "poped", + popedControl(groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24), + c(DOSE=40, TAU=24)), + maxa=c(DOSE=200,TAU=24), + mina=c(DOSE=0,TAU=24))) + + babelDesign <- evaluate_design(babel.db) + + ##-- Model: One comp first order absorption + ## -- Analytic solution for both mutiple and single dosing + ff <- function(model_switch,xt,parameters,poped.db) { with(as.list(parameters),{ y=xt N = floor(xt/TAU)+1 - f=(DOSE/V)*(KA/(KA - CL/V)) * + y=(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)) + return(list( y=y,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]) + ## -- parameter definition function + ## -- names match parameters in function ff + 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]) return( parameters ) } + ## -- Residual unexplained variablity (RUV) function + ## -- Additive + Proportional feps <- function(model_switch,xt,parameters,epsi,poped.db){ - f <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))[[1]] - y = f*(1+epsi[,1])+epsi[,2] - return(list(y=y,poped.db=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_analytic <- create.poped.database( - ff_fun =ff_analytic, - fg_fun =sfg, - fError_fun=feps, - bpop=c(KA=0.25,CL=3.75,V=72.8), - d=c(KA=0.09,CL=0.25^2,V=0.09), - sigma=c(prop=0.04,add=0.0025), - m=2, - groupsize=20, - xt=c( 1,2,8,240,245), - minxt=c(0,0,0,240,240), - maxxt=c(10,10,10,248,248), - bUseGrouped_xt=1, - a=cbind(DOSE=c(20,40),TAU=c(24,24)), - maxa=c(DOSE=200,TAU=24), - mina=c(DOSE=0,TAU=24)) - - - e <- et(amt=1, ii=24, until=250) %>% - et(list(c(0, 10), - c(0, 10), - c(0, 10), - c(240, 248), - c(240, 248))) %>% - dplyr::mutate(time =c(0, 1,2,8,240,245)) - - # model + ## -- Define design and design space + poped.db <- create.poped.database(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(0.04,5e-6), + notfixed_sigma=c(0,0), + m=2, + groupsize=20, + xt=c( 1,2,8,240,245), + minxt=c(0,0,0,240,240), + maxxt=c(10,10,10,248,248), + bUseGrouped_xt=1, + a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), + maxa=c(DOSE=200,TAU=24), + mina=c(DOSE=0,TAU=24)) + + popedDesign <- evaluate_design(poped.db) + + expect_equal(babelDesign$ofv, popedDesign$ofv) + + }) + + test_that("example 2", { + + library(PopED) + f <- function() { ini({ - tKA <- 0.25 - tCL <- 3.75 - tV <- 72.8 - Favail <- fix(0.9) - eta.ka ~ 0.09 - eta.cl ~ 0.25 ^ 2 - eta.v ~ 0.09 - prop.sd <- sqrt(0.04) - add.sd <- sqrt(5e-6) + tCl <- 0.15 + tV <- 8 + tKA <- 1.0 + tFavail <- fix(1) + eta.cl ~ 0.07 + eta.v ~ 0.02 + eta.ka ~ 0.6 + prop.sd <- sqrt(0.01) + }) + model({ + CL <- tCl*exp(eta.cl) + V <- tV*exp(eta.v) + KA <- tKA*exp(eta.ka) + Favail <- tFavail + y <- (DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*time)-exp(-KA*time)) + y ~ prop(prop.sd) + }) + } + + e <- et(c(0.5, 1,2,6,24,36,72,120)) %>% + as.data.frame() + + ## -- Define initial design and design space + babel.db <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=32, + minxt=0, + maxxt=120, + a=70)) + + babelPred <- model_prediction(babel.db) + + babelED <- withr::with_seed(42, { + evaluate_design(babel.db,d_switch=FALSE,ED_samp_size=20) + }) + + + # Example changes the model to add+prop error model + + f <- function() { + ini({ + tCl <- 0.15 + tV <- 8 + tKA <- 1.0 + tFavail <- fix(1) + eta.cl ~ 0.07 + eta.v ~ 0.02 + eta.ka ~ 0.6 + prop.sd <- sqrt(0.01) + add.sd <- sqrt(0.25) }) model({ - ka <- tKA * exp(eta.ka) - v <- tV * exp(eta.v) - cl <- tCL * exp(eta.cl) - d/dt(depot) <- -ka * depot - d/dt(central) <- ka * depot - cl / v * central - cp <- central / v - f(depot) <- DOSE * Favail - cp ~ add(add.sd) + prop(prop.sd) + CL <- tCl*exp(eta.cl) + V <- tV*exp(eta.v) + KA <- tKA*exp(eta.ka) + Favail <- tFavail + y <- (DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*time)-exp(-KA*time)) + y ~ prop(prop.sd) + add(add.sd) + }) + } + + e <- et(c(0.5, 1,2,6,24,36,72,120)) %>% + as.data.frame() + + babel.db <- nlmixr2(f, e, "poped", + popedControl(groupsize=32, + minxt=0, + maxxt=120, + a=70, + mina=0, + maxa=100, + # selecting important/unimportant + # parameters assumes Ds optimal design. + important=c("tCl", "tV", "tKa"))) + + + babelDs <- withr::with_seed(42, evaluate_design(babel.db)) + + ff <- function(model_switch,xt,parameters,poped.db){ + ##-- Model: One comp first order absorption + with(as.list(parameters),{ + y=xt + y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt)) + return(list(y=y,poped.db=poped.db)) }) } - poped_db_ode_babelmixr2 <- nlmixr(f, e, - popedControl(a=list(c(DOSE=20), - c(DOSE=40)), - maxa=c(DOSE=200), - mina=c(DOSE=0))) - - eb <- evaluate_design(poped_db_ode_babelmixr2) - - expect_equal(list(ofv = 58.595188415962, - fim=lotri({ - tKA + tCL + tV ~ c(2999.85286248872, - -14.4306063906335, 37.1524325291608, - -8.68396393002655, -0.0586341331887462, 0.0533669124790745) - d_eta.ka + d_eta.cl + d_eta.v + sig_prop.sd + sig_add.sd ~ - c(439.413101383167, - 2.2878439908508, 3412.00514559432, - 312.24028527664, 3.20285380277308, 999.953201635217, - 638.582017761863, 1182.32547832173, 575.347350448117, 33864.3205676804, - 82065.2676506916, 38107.6501681295, 21309.2735081853, 2327835.71273584, 404178779.571231) - }), - rse = c(tKA = 10.0909508715706, tCL = 4.40030409644082, tV = 8.21533739750444, - d_eta.ka = 61.371504736694, d_eta.cl = 27.592223215635, - d_eta.v = 40.0643564351632, - sig_prop.sd = 17.6873376767579, sig_add.sd = 1297.44406776262)), - eb) + sfg <- function(x,a,bpop,b,bocc){ + ## -- parameter definition function + 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) + } + + feps <- function(model_switch,xt,parameters,epsi,poped.db){ + ## -- Residual Error function + ## -- Proportional + 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]) + return(list(y=y,poped.db=poped.db)) + } + + ## -- Define initial design and design space + poped.db <- create.poped.database(ff_file="ff", + fg_file="sfg", + fError_file="feps", + 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=0.01, + groupsize=32, + xt=c( 0.5,1,2,6,24,36,72,120), + minxt=0, + maxxt=120, + a=70) + + popedPred <- model_prediction(poped.db) + + popedED <- withr::with_seed(42, { + evaluate_design(poped.db,d_switch=FALSE,ED_samp_size=20) + }) + + expect_equal(popedPred$ofv, babelPred$ofv) + + expect_equal(popedED$ofv, babelED$ofv) + + poped.db <- create.poped.database(ff_fun=ff, + 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=70, + mina=0, + maxa=100, + ds_index=c(0,0,0,1,1,1,1,1), + ofv_calc_type=6) + popedDs <- withr::with_seed(42, evaluate_design(poped.db)) + + expect_equal(popedDs$ofv, babelDs$ofv) + + }) + + test_that("example 3", { }) }