mrgsolve with poped

[vjd]

this is with reference to a discussion on the poped repo.

how do we add covariate functions in a mrgsolve model that can be used to optimize parameters in poped? I went back and looked at the examples provided in the poped repo, lines 39-40

Based on this, it seems that the ff function is probably the place to include the covariate form. We cannot really define the covariate form in the mrgsolve code, as parameters defined in $PARAM get passed into the ff function via as.list(p), and there is no way to reassign a parameter (defined in $PARAM) in $MAIN to be function of a covariate, like this below

code <- '
$PARAM CL=10, V=2, WT=70

$MAIN
double CL = CL*(WT/70);
'

mrgsolve will throw an error that a variable defined in the parameter/dataset cannot be reassigned in $MAIN I don't know if my understanding is incorrect, but still trying to explore this integration.

[kylebaron]

Here's what I would do:

• Use CL, V in the PopED part (as you have it now)
• Define covariate model in $MAIN, deriving Pi (i.e. CLi, Vi etc) ....
• You only have to do this for parameters that get into the covariate model; if a parameter isn't in the covariate model, just name the parameter in $PARAM the same name that you are using on the PopED side. It will update and you can use that symbol directly.
• Use Pi to advance the system

Let me know if this is confusing or you just don't want to do it this way. This was the least invasive (IMO) that I could think of. There are other ways, but it's all just a naming game.

code <- '

$PARAM CL = 10, V = 2, WT = 70

$MAIN
double CLi = CL*pow(WT/70,0.75);
double Vi = V*(WT/70);

$ODE
dxdt_CENT = -(CLi/VCi)*CENT;

$TABLE
table(CP) = CENT/Vi;
'

This is just a matter of picking a naming convention that will work with everything and is (hopefully) semi-informative.

I think the piece that is confusing here is:

sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                VC=bpop[2]*exp(b[2]),
                Q=bpop[3],
                VP=bpop[4],
                TVLAG=bpop[5],
                DOSE=a[1],
                WT=a[2])
  return( parameters ) 
}

I get that bpop[1] is TVCL, but it would seem like CL now is the individual clearance. But it isn't because we're putting off the covariate mode until later. But that's okay; just pick a new parameter name to derive inside the model for use to advance the system.

[vjd]

I tried out the naming convention you suggested and it seems to work. I missed the point that "we're putting off the covariate model until later". Taking this point into account, the understanding of the poped workflow is a little better.

[smouksassi]

sorry to revive an old thread is it worth it to update the vignette on the wiki ? https://github.com/metrumresearchgroup/mrgsolve/wiki/PopED_vignette . This comes every time someone tries to use mrgsolve with POPED and the old wiki vignette has no covariate example

[kylebaron]

Hi @smouksassi - I'd be happy to revisit the wiki vignette (it does look really old). Will just need to work back through the models to see what the conclusion was for what should be done (not sure what I was talking about in the previous examples).

KTB

[smouksassi]

code that was discussed with parami notation with sfg realizing the cov effect

code_mod_cov <- '
$CMT DEPOT CENT PERI
$PARAM
CL=20, VC=70,Q=1, VP=30, KA = 0.25, WT= 70

$MAIN
double CLi = CL*pow(WT/70,0.75);
double VCi = VC*(WT/70);
double Qi  = Q*pow(WT/70,0.75);
double VPi = VP*(WT/70);

$ODE
dxdt_DEPOT =                                                       - KA*DEPOT;
dxdt_CENT = - (CLi /VCi )*CENT -(Qi /VCi )*CENT  + (Qi /VPi )*PERI + KA*DEPOT;
dxdt_PERI =                     (Qi /VCi )*CENT  - (Qi /VPi )*PERI;
$TABLE double CP = CENT/VCi ;
$CAPTURE CP
'
mod_cov <- mread(code=code_mod_cov, model="excov_mrg_poped")

data <- ev(amt=c(4,4,8,15), cmt=1) %>% as.data.frame %>% as.tbl %>% mutate(ID=1:4, dose = amt*10,amt=amt*10000,
                                                                           WT=c(35,70,70,70))
out <- mod_cov %>% data_set(data) %>% Req(CP) %>% carry.out(dose) %>% mrgsim(end=48, delta=0.1)

plot(out,logy=TRUE)

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 ) 
}

poped.cov_mrg_poped <- create.poped.database(ff_file="ff.mrgsolve",
                                               fg_file="sfg.mrgsolve.cov",
                                               fError_file="feps.log",
                                               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_model_prediction(poped.cov_mrg_poped,model_num_points = 5000)

### Evaluate FIM
FIM <- evaluate.fim(poped.cov_mrg_poped) 
det(FIM)
get_rse(FIM, poped.nocov_mrg_poped)

[kylebaron]

@smouksassi https://github.com/metrumresearchgroup/mrgsolve/wiki/PopED_vignette

[smouksassi]

Looks great thank you we can close the issue and link to new vignette !

[timwaterhouse]

I'm not sure if I'm reading this correctly, but in the vignette on the wiki it looks to me like you're doing the body weight scaling twice here: once in sfg.mrgsolve.cov and again in the mrgsolve model. Should the covariate effect be removed from the fg function (as it is in this PopED vignette)?

[kylebaron]

Thanks, @timwaterhouse. Can you correct the code, please?

#' ---
#' title: ""
#' author: Samer Mouksassi and Kyle Baron
#' output: github_document
#' ---

#+ include = FALSE
knitr::opts_chunk$set(comment = '.', warning = FALSE, message = FALSE)
#+

#' # Setup
library(tidyverse)
library(mrgsolve)
library(PopED)

#' # The model
mod <- mread("pop_ed_model.txt")

#+ c, eval=FALSE, code = mod@code

#+
#' 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)
get_rse(FIM, poped.cov_mrg_poped)

[ 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 ;

[timwaterhouse]

Sure. The model is unchanged, but I've removed the weight effect from sfg.mrgsolve.cov.

#' ---
#' title: ""
#' author: Samer Mouksassi and Kyle Baron
#' output: github_document
#' ---

#+ include = FALSE
knitr::opts_chunk$set(comment = '.', warning = FALSE, message = FALSE)
#+

#' # Setup
library(tidyverse)
library(mrgsolve)
library(PopED)

#' # The model
mod <- mread("pop_ed_model.txt")

#+ c, eval=FALSE, code = mod@code

#+
#' 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]*exp(b[1]),
               VC=bpop[2]*exp(b[2]),
               Q=bpop[3],
               VP=bpop[4],
               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)
get_rse(FIM, poped.cov_mrg_poped)

This makes a pretty big difference to the model prediction for the lower weight group, but only a negligible effect on the RSEs.

[kylebaron]

thanks

[kylebaron]

Let re-open if this isn't correct https://github.com/metrumresearchgroup/mrgsolve/wiki/PopED_vignette