popED

PopED example with a covariate

Samer Mouksassi and Kyle Baron 2019-07-03

Setup

library(tidyverse)
library(mrgsolve)
library(PopED)

The model

mod <- mread("pop_ed_model.txt")

[ pkmodel ] cmt="DEPOT CENT PERI", depot=TRUE, trans = 11

[ param ]
CL=20, VC=70, Q=1, VP=30, KA = 0.25, WT= 70

[ main ]
double CLi = CL*pow(WT/70,0.75);
double V2i = VC*(WT/70);
double Qi  = Q*pow(WT/70,0.75);
double V3i = VP*(WT/70);
double KAi = KA;

[ table ]
capture CP = CENT/V2i ;

The important part here is that we have a PK model parameters that are a function of a covariate (WT) and WT is listed as a parameter.

Example simulation

data <- 
  expand.ev(amt=c(4,4,8,15)) %>% 
  mutate(ID=1:4, dose = amt*10, amt=amt*10000, WT=c(35,70,70,70))

out <- 
  mod %>% 
  data_set(data) %>% 
  Req(CP) %>% 
  carry_out(dose) %>% 
  mrgsim(end=48, delta=0.1)

plot(out,logy=TRUE)

The PopED setup

I am using mrgsim_q to get the simulation turned around as quickly as possible, reducing overhead (and dropping features).

ff <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)  
  dose_times <- 0
  time <- sort(unique(c(times_xt,dose_times)))
  is.dose <- time %in% dose_times
  
  data <- data.frame(
    ID = 1,
    time = time,
    amt = ifelse(is.dose,p[["DOSE"]], 0), 
    cmt = ifelse(is.dose, 1, 0)
  )
  
  data[["evid"]] <- data[["cmt"]]
  
  mod <- param(mod, as.list(p))
  
  out <- mrgsim_q(mod,data,recsort=4)
  
  y <- out$CP[match(times_xt,out$time)]
  
  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
}

sfg.mrgsolve.cov <- function(x,a,bpop,b,bocc){
  parameters=c(CL=bpop[1]*((a[3]/70)^0.75) *exp(b[1]),
               VC=bpop[2]*((a[3]/70)^1)*exp(b[2]),
               Q=bpop[3]*((a[3]/70)^0.75),
               VP=bpop[4]*((a[3]/70)^1),
               KA=bpop[5],
               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*exp(epsi[,1])
  return(list( y= y,poped.db =poped.db )) 
}

poped.cov_mrg_poped <- create.poped.database(
  ff_fun=ff,
  fg_fun=sfg.mrgsolve.cov,
  fError_fun=feps,
  bpop=c(CL=20, VC=70,Q=1, VP=30, KA=0.25), 
  notfixed_bpop=c(1,1,1,0,0),
  d=c(CL=0.08,VC=0.1), 
  sigma=c(0.05),
  notfixed_sigma=c(1),
  m=2,
  groupsize=6,
  xt=list(c(3,5,8,18),c(1,6,50,120)),
  minxt=0,
  maxxt=120,
  bUseGrouped_xt=0,
  a = cbind(DOSE=c(250,250), TAU=c(12,12),WT=c(35,70))
)

Plot a model prediction

plot_model_prediction(poped.cov_mrg_poped,model_num_points = 5000) + theme_bw()

Evaluate FIM

FIM <- evaluate.fim(poped.cov_mrg_poped) 
det(FIM)

. [1] 961743530

get_rse(FIM, poped.cov_mrg_poped)

.    bpop[1]    bpop[2]    bpop[3]     D[1,1]     D[2,2] SIGMA[1,1] 
.   8.627307  12.026370   6.253843  44.933063  72.376153  28.601719