search

metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
129 stars 36 forks source link

MAP Bayes estimation with multiple analytes (parent drug and metabolites) #149

Closed MarioGS closed 7 years ago

MarioGS commented 7 years ago

Hi,

I am trying to do MAP Bayes estimation with mrgsolve of a complex PK model. The model has 7 compartments. Compartments 1-3 are from parent drug (CPT-11). Compartments 4 and 5 are for one metabolite (SN-38) and compartments 6 and 7 for a second metabolite (SN-38G). The model is the one proposed by Xie in 2002 and published in Clinical Pharmacology & Therapeutics. I implemented the model in mrgsolve as follows:

code <- '
  $PARAM TVCL = 16.9, TVQ2 = 61, TVQ3 = 4.75, TVV1 = 36.7, TVV2 = 35.9, TVV3 = 67.9,
  TVCLSN38=712, TVQ5=1530, TVV4=408, TVV5=71600,
  TVCLSN38G=66.8, TVQ7=23.9, TVV6=37, TVV7=48.4, BSA=1.87
  ETA1 = 0, ETA2 = 0, ETA3 = 0, ETA4 = 0, ETA5 = 0, ETA6 = 0 
  ETA7 = 0, ETA8 = 0, ETA9 = 0, ETA10 = 0, 
  ETA11 = 0, ETA12 = 0, ETA13 = 0, ETA14 = 0 
  $CMT A1 A2 A3 A4 A5 A6 A7
  $MAIN
  double CL = BSA*TVCL*exp(ETA1 + ETA(1));
  double Q2 = BSA*TVQ2*exp(ETA2 + ETA(2));
  double Q3 = BSA*TVQ3*exp(ETA3 + ETA(3));
  double V1 = BSA*TVV1*exp(ETA4 + ETA(4));
  double V2 = BSA*TVV2*exp(ETA5 + ETA(5));
  double V3 = BSA*TVV3*exp(ETA6 + ETA(6));
  double CLSN38 = TVCLSN38*exp(ETA7 + ETA(7));
  double Q5 = TVQ5*exp(ETA8 + ETA(8));
  double V4 = TVV4*exp(ETA9 + ETA(9));
  double V5 = TVV5*exp(ETA10 + ETA(10));
  double CLSN38G = TVCLSN38G*exp(ETA11 + ETA(11));
  double Q7 = TVQ7*exp(ETA12 + ETA(12));
  double V6 = TVV4*exp(ETA13 + ETA(13));
  double V7 = TVV5*exp(ETA14 + ETA(14));
  $ODE
  dxdt_A1 = -(CL/V1)*A1 - (Q2/V1)*A1 - (Q3/V1)*A1 + (Q2/V2)*A2 + (Q3/V3)*A3;
  dxdt_A2 =               (Q2/V1)*A1              - (Q2/V2)*A2;
  dxdt_A3 =                            (Q3/V1)*A1              - (Q3/V3)*A3;
  dxdt_A4 =  (CL/V1)*A1 - (CLSN38/V4)*A4 -  (Q5/V4)*A4 + (Q5/V5)*A5;
  dxdt_A5 =                                 (Q5/V4)*A4 - (Q5/V5)*A5;
  dxdt_A6 =               (CLSN38/V4)*A4 - (CLSN38G/V6)*A6 - (Q7/V6)*A6 + (Q7/V7)*A7;
  dxdt_A7 =                                                  (Q7/V6)*A6 - (Q7/V7)*A7;
  $OMEGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  $SIGMA 0 0 0 0 0
  $TABLE
  double DV = (A1/V1)*(1+EPS(1))+EPS(2);
  double DV4 = (A4/V4)*(1+EPS(3))+EPS(4);
  double DV6 = (A6/V6)+EPS(5);
  double ET1 = ETA(1);
  double ET2 = ETA(2);
  double ET3 = ETA(3);
  double ET4 = ETA(4);
  double ET5 = ETA(5);
  double ET6 = ETA(6);
  double ET7 = ETA(7);
  double ET8 = ETA(8);
  double ET9 = ETA(9);
  double ET10 = ETA(10);
  double ET11 = ETA(11);
  double ET12 = ETA(12);
  double ET13 = ETA(13);
  double ET14 = ETA(14);
  $CAPTURE DV DV4 DV6 ET1 ET2 ET3 ET4 ET5 ET6 ET7 ET8 ET9 ET10 ET11 ET12 ET13 ET14
  $SET delta=0.1, end=48      
  '

However, since I have observations in compartments 1, 4 and 6, I am hesitating about how should I define sigma (because of the proportional + additive error) and how the should I modify mapbayes function: 

    mapbayes <- function(eta,y,d,m,pred=FALSE) {
      sig2 <- as.numeric(sigma)
      eta %<>% as.list
      names(eta) <- names(init)
      eta_m <- eta %>% unlist %>% matrix(nrow=1)
      m %<>% param(eta)
      if(!pred) m %<>% obsonly
      out <- m %>% drop.re() %>% data_set(d) %>% mrgsim
      if(pred) return(out %>% as.tbl)
      # http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/
      sig2j <- out$DV^2*sig2
      sqwres <- log(sig2j) + (1/sig2j)*(y-out$DV)^2
      nOn <- diag(eta_m %*% omega.inv %*% t(eta_m))
      return(sum(sqwres) + nOn)
    }

in order to make it compatible with such a model.

Any suggestions?

Thanks in advance,

Mario

kylebaron commented 7 years ago

@MarioGS -

Here are some ideas that I've tried before. Are you going to be at ACoP? Can we meet there and discuss?

Additive and proportional error

You can modify sig2j <- out$DV^2*sig2 to be sig2j <- out$DV^2*sig2prop + sig2add. With 3 species I imagine that there are three sets of error variances?

It's a matter of bookkeeping ... here's one way:

When I have done multiple endpoints before, I ended up switching to long format. mrgsolve is wide (you get predictions for all the compartments on every record). But you can do it long too:

$TABLE
double DV = A2;
if(type==2) DV = A3;
if(type==3) DV = A4;

$CAPTURE DV type

That makes the bookkeeping a little easier, IMO: you do a little work a head of time, but when it comes time to calculate likelihood, you can still pick off out$DV and get the predictions in the same order as the observed values.

One approach

One way I've done it is to get an integer vector of the rows in the input data for the different species. Something like spec1 <- which(data$type==1), spec2 <- which(data$type==2) if data is the input data set. Then when you take out$DV[spec1] for the first species and out$DV[spec2] for second species etc ... then, just calculate the likelihoods separately (for each species) with simple scalar variances and add them all up at the end.

This approach may not require the data to be long (not sure ... ). I wrote this up last but it might be the easiest way to go. We had a model with 4 endpoints recently and took this approach.

Another approach

In the vignette online, sig2 was a simple scalar. Another way to get the right variance with the right species in this example, I'd make type to be c(1,2,3) and make sig2prop to be c(sig2prop1,sig2prop2,sig2prop3) and similar vector for additive error. Then you can get the vector of variances that matches the observation by sig2prop[out$type] and sig2add[out$type]. This is assuming that you can filter the simulated output down to observations only so that out$type is always 1 or 2 or 3.

Happy to look at code or try something out with you if you can keep sharing the progress.

Kyle

kylebaron commented 7 years ago

@MarioGS Also ... (you might already be doing this). Make a data set with all of the observation records at their appropriate times. Don't let mrgsolve pad in observations for you; just get them where you have data. Maybe already doing that but I saw $SET delta=0.1, end=48 and wasn't sure.

MarioGS commented 7 years ago

Thanks for your quick reply Kyle.

Unfortunately, I will miss this ACoP :(

I am not sure how should I "play" with the omega and sigma when I filtering by "TYPE" and using mapbayes. Here is the code I tried:

### Packages
library(ggplot2)
library(plyr)
library(reshape2)
library(compiler)
library(DT)
library(dplyr)
library(minqa)
library(mrgsolve)
library(magrittr)

########################### Generate data starts here
## ATENTION: dose in ug, volume in L, concentration in ug/L = ng/mL
dose = c(300000) #"Dose" (ug)
dur  = c(1) #Infusion duration (hours)
time = c(0,rep(0.25,times=3),rep(1,times=3),rep(1.5,times=3),rep(2,times=3),
         rep(3,times=3),rep(7,times=3),rep(27,times=3),rep(46,times=3)) #Sampling times (hours)
dat = as.data.frame(expand.grid("dose"=dose,"dur"=dur,"time"=time)) #Expand grid
dat = arrange(dat,time)
dat$TYPE = c(0,rep(seq(1:3),8))
dat$amt = ifelse(dat$time==0,dat$dose,0) #Generate AMT
dat$DV = c(0,1686.090652,12.5616153,43.45757775,1590.931119,16.88053955,77.58447629,
           1497.367194,13.42229193,91.09384167,1260.858318,10.75658997,92.01476183,
           914.2540691,7.307412688,82.5660438,492.590467,6.095285892,60.06956412,
           351.4107986,2.698844692,21.16848336,106.7711819,2.148240981,13.69989089
           )
dat$evid = ifelse(dat$amt>0,1,0)  
dat$ID = 1
dat$cmt = ifelse(dat$amt>0,1,0) 
dat$BSA = 1.8
########################### Generate data end here

########################### Model definition starts here
code <- '
$PARAM TVCL = 16.9, TVQ2 = 61, TVQ3 = 4.75, TVV1 = 36.7, TVV2 = 35.9, TVV3 = 67.9,
TVCLSN38=712, TVQ5=1530, TVV4=408, TVV5=71600,
TVCLSN38G=66.8, TVQ7=23.9, TVV6=37, TVV7=48.4, BSA=1.87, TYPE=1,
ETA1 = 0, ETA2 = 0, ETA3 = 0, ETA4 = 0, ETA5 = 0, ETA6 = 0 
ETA7 = 0, ETA8 = 0, ETA9 = 0, ETA10 = 0, 
ETA11 = 0, ETA12 = 0, ETA13 = 0, ETA14 = 0 
$CMT A1 A2 A3 A4 A5 A6 A7
$MAIN
double CL = BSA*TVCL*exp(ETA1 + ETA(1));
double Q2 = BSA*TVQ2*exp(ETA2 + ETA(2));
double Q3 = BSA*TVQ3*exp(ETA3 + ETA(3));
double V1 = BSA*TVV1*exp(ETA4 + ETA(4));
double V2 = BSA*TVV2*exp(ETA5 + ETA(5));
double V3 = BSA*TVV3*exp(ETA6 + ETA(6));
double CLSN38 = TVCLSN38*exp(ETA7 + ETA(7));
double Q5 = TVQ5*exp(ETA8 + ETA(8));
double V4 = TVV4*exp(ETA9 + ETA(9));
double V5 = TVV5*exp(ETA10 + ETA(10));
double CLSN38G = TVCLSN38G*exp(ETA11 + ETA(11));
double Q7 = TVQ7*exp(ETA12 + ETA(12));
double V6 = TVV4*exp(ETA13 + ETA(13));
double V7 = TVV5*exp(ETA14 + ETA(14));
$ODE
dxdt_A1 = -(CL/V1)*A1 - (Q2/V1)*A1 - (Q3/V1)*A1 + (Q2/V2)*A2 + (Q3/V3)*A3;
dxdt_A2 =               (Q2/V1)*A1              - (Q2/V2)*A2;
dxdt_A3 =                            (Q3/V1)*A1              - (Q3/V3)*A3;
dxdt_A4 =  (CL/V1)*A1 - (CLSN38/V4)*A4 -  (Q5/V4)*A4 + (Q5/V5)*A5;
dxdt_A5 =                                 (Q5/V4)*A4 - (Q5/V5)*A5;
dxdt_A6 =               (CLSN38/V4)*A4 - (CLSN38G/V6)*A6 - (Q7/V6)*A6 + (Q7/V7)*A7;
dxdt_A7 =                                                  (Q7/V6)*A6 - (Q7/V7)*A7;
$OMEGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0
$SIGMA 0 0 0 0 0
$TABLE
double DV = (A1/V1)*(1+EPS(1))+EPS(2);
if(TYPE==2) DV = (A4/V4)*(1+EPS(3))+EPS(4);
if(TYPE==3) DV = (A6/V6)+EPS(5);
double ET1 = ETA(1);
double ET2 = ETA(2);
double ET3 = ETA(3);
double ET4 = ETA(4);
double ET5 = ETA(5);
double ET6 = ETA(6);
double ET7 = ETA(7);
double ET8 = ETA(8);
double ET9 = ETA(9);
double ET10 = ETA(10);
double ET11 = ETA(11);
double ET12 = ETA(12);
double ET13 = ETA(13);
double ET14 = ETA(14);
$CAPTURE DV TYPE ET1 ET2 ET3 ET4 ET5 ET6 ET7 ET8 ET9 ET10 ET11 ET12 ET13 ET14
'
########################### Model definition ends here

#compile model
mod <- mcode("map", code) 

#Define diagonal matrix
omega <- dmat(0.32^2,0.74^2,0.47^2,0.26^2,0.27^2,0.45^2, #IIV CPT-11 as (%CV)
              0.64^2,0.35^2,0.29^2,0.31^2,               #IIV SN38 as (%CV)
              0.80^2,0.97^2,0.75^2,0.59^2)               #IIV SN38G as (%CV)
omega.inv <- solve(omega) 

#Define proportional error for cpt-11
sigma1prop <- matrix(0.167^2) #(16.7 %CV) proportional error for cpt-11 compartment 1
#Define additive error for cpt-11
sigma1add <- matrix(1.3) #(ng/ml) additive error for cpt-11 compartment 1

#Define proportional error for sn-38
sigma4prop <- matrix(0.268^2) #(26.8 %CV) proportional error for sn38 compartment 4
#Define additive error for sn-38
sigma4add <- matrix(0.56) #(ng/ml) additive error for sn38 compartment 4

#Define additive error for sn-38g
sigma6add <- matrix(34.5) #(ng/ml) additive error for sn38g compartment 6

#Generate y for mapbayes function
y1 <- dat %>% filter(evid==0) %>% filter(TYPE==1) %>% select(DV) %>% unlist %>% unname
y1
y4 <- dat %>% filter(evid==0) %>% filter(TYPE==2) %>% select(DV) %>% unlist %>% unname
y4
y6 <- dat %>% filter(evid==0) %>% filter(TYPE==3) %>% select(DV) %>% unlist %>% unname
y6

#Generate data for mapbayes function
data = dat
data1 = filter(data,TYPE<=1)
data4 = filter(data,TYPE==2 | TYPE==0)
data6 = filter(data,TYPE==3 | TYPE==0)

#mapbayes for cpt-11 [PROP + ADD RV CPT-11]
mapbayes1 <- function(eta,y,d,m,pred=FALSE) {
  sig2prop <- as.numeric(sigma1prop)
  sig2add <- as.numeric(sigma1add)
  eta %<>% as.list
  names(eta) <- names(init)
  eta_m <- eta %>% unlist %>% matrix(nrow=1)
  m %<>% param(eta)
  if(!pred) m %<>% obsonly
  out <- m %>% drop.re() %>% data_set(d) %>% mrgsim
  if(pred) return(out %>% as.tbl)
  # http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/
  #sig2j <- out$DV^2*sig2
  sig2j <- out$DV^2*sig2prop + sig2add
  sqwres <- log(sig2j) + (1/sig2j)*(y-out$DV)^2
  nOn <- diag(eta_m %*% omega.inv %*% t(eta_m))
  return(sum(sqwres) + nOn)
}

#mapbayes for SN-38 [PROP + ADD RV SN-38]
mapbayes4 <- function(eta,y,d,m,pred=FALSE) {
  sig2prop <- as.numeric(sigma4prop)
  sig2add <- as.numeric(sigma4add)
  eta %<>% as.list
  names(eta) <- names(init)
  eta_m <- eta %>% unlist %>% matrix(nrow=1)
  m %<>% param(eta)
  if(!pred) m %<>% obsonly
  out <- m %>% drop.re() %>% data_set(d) %>% mrgsim
  if(pred) return(out %>% as.tbl)
  # http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/
  #sig2j <- out$DV^2*sig2
  sig2j <- out$DV^2*sig2prop + sig2add
  sqwres <- log(sig2j) + (1/sig2j)*(y-out$DV)^2
  nOn <- diag(eta_m %*% omega.inv %*% t(eta_m))
  return(sum(sqwres) + nOn)
}

#mapbayes for SN-38G [ADD RV sn38G]
mapbayes6 <- function(eta,y,d,m,pred=FALSE) {
  sig2add <- as.numeric(sigma6add)
  eta %<>% as.list
  names(eta) <- names(init)
  eta_m <- eta %>% unlist %>% matrix(nrow=1)
  m %<>% param(eta)
  if(!pred) m %<>% obsonly
  out <- m %>% drop.re() %>% data_set(d) %>% mrgsim
  if(pred) return(out %>% as.tbl)
  # http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339294/
  #sig2j <- out$DV^2*sig2
  sig2j <- out$DV^2 + sig2add
  sqwres <- log(sig2j) + (1/sig2j)*(y-out$DV)^2
  nOn <- diag(eta_m %*% omega.inv %*% t(eta_m))
  return(sum(sqwres) + nOn)
}

#Bayesian estimation for CPT-11
# I do not care about ETA>7
init1 <- c(ETA1=-0.3, ETA2=0.2,ETA3=-0.3, ETA4=0.2,ETA5=-0.3, ETA6=0.2,
          ETA7=0, ETA8=0,ETA9=0, ETA10=0,
          ETA11=0, ETA12=0,ETA13=0, ETA14=0)
fit1 <- newuoa(init1,mapbayes1,y=y1,m=mod,d=data1)
fit1$par[1] #IIV CL
fit1$par[2] #IIV Q2
fit1$par[3] #IIV Q3
fit1$par[4] #IIV V1
fit1$par[5] #IIV V2
fit1$par[6] #IIV V3

#Bayesian estimation for SN-38
# I do not care about ETA<7 & ETA>10
init4 <- c(ETA1=0.3, ETA2=0.3, ETA3=0.1, ETA4=0.1, ETA5=0.1, ETA6=0.1,
           ETA7=0.1, ETA8=0.1,ETA9=0.1, ETA10=0.1,
           ETA11=0, ETA12=0,ETA13=0, ETA14=0)
fit4 <- newuoa(init4,mapbayes4,y=y4,m=mod,d=data4)
fit4$par[7] #IIV CLSN38
fit4$par[8] #IIV Q5
fit4$par[9] #IIV V4
fit4$par[10] #IIV V5

#Bayesian estimation for SN-38G
# I do not care about ETA<10
init6 <- c(ETA1=0, ETA2=0,ETA3=0, ETA4=0,ETA5=0, ETA6=0,
           ETA7=0, ETA8=0,ETA9=0, ETA10=0,
           ETA11=0.1, ETA12=0.1,ETA13=0.1, ETA14=0.1)
fit6 <- newuoa(init6,mapbayes6,y=y6,m=mod,d=data6)
fit6$par[11] #IIV CLSN38G
fit6$par[12] #IIV Q7
fit6$par[13] #IIV V6
fit6$par[14] #IIV V7

#Compute individual parameters
data$CL = data$BSA[1]*16.9*exp(fit1$par[1])
data$Q2 = data$BSA[1]*61*exp(fit1$par[2])
data$Q3 = data$BSA[1]*4.75*exp(fit1$par[3])
data$V1 = data$BSA[1]*36.7*exp(fit1$par[4])
data$V2 = data$BSA[1]*35.9*exp(fit1$par[5])
data$V3 = data$BSA[1]*67.9*exp(fit1$par[6])
data$CLSN38 = 712*exp(fit4$par[7])
data$Q5 = 1530*exp(fit4$par[8])
data$V4 = 408*exp(fit4$par[9])
data$V5 = 71600*exp(fit4$par[10])
data$CLSN38G = 66.8*exp(fit6$par[7])
data$Q7 = 23.9*exp(fit6$par[8])
data$V6 = 37*exp(fit6$par[9])
data$V7 = 48.4*exp(fit6$par[10])

#Update parameters to perform simulation
mod_alt <- param(mod,
                 TVCL = data$CL[1],
                 TVQ2 = data$Q2[1],
                 TVQ3 = data$Q3[1],
                 TVV1 = data$V1[1],
                 TVV2 = data$V2[1],
                 TVV3 = data$V3[1],
                 TVCLSN38 = data$CLSN38[1],
                 TVQ5 = data$Q5[1],
                 TVV4 = data$V4[1],
                 TVV5 = data$V5[1],
                 TVCLSN38G = data$CLSN38G[1],
                 TVQ7 = data$Q7[1],
                 TVV6 = data$V6[1],
                 TVV7 = data$V7[1])

DOSE = data$dose[1]
RATE = data$dose[1]/data$dur

#Generate tbl with results from the final simulation
final = mod_alt %>%
  ev(amt=DOSE,rate=RATE) %>%
  mrgsim(nid=1, start=0, end=50) %>%
  as.tbl()

final$DV1 = final$A1/data$V1[1]
final$DV4 = final$A4/data$V4[1]
final$DV6 = final$A6/data$V6[1]

#Plot CPT-11
ggplot() +
  geom_point(data=subset(dat,TYPE==1),aes(time,DV))+
  geom_line(data=final,aes(time,DV1))+
  scale_y_log10()+
  annotation_logticks(sides = "l")+
  ylab("CPT-11 [ug/l]")+
  xlab("Time [h]")

#Plot SN38
ggplot() +
  geom_point(data=subset(dat,TYPE==2),aes(time,DV))+
  geom_line(data=final,aes(time,DV4))+
  scale_y_log10()+
  annotation_logticks(sides = "l")+
  ylab("SN-38 [ug/l]")+
  xlab("Time [h]")

#Plot SN38G
ggplot() +
  geom_point(data=subset(dat,TYPE==3),aes(time,DV))+
  geom_line(data=final,aes(time,DV6))+
  scale_y_log10()+
  annotation_logticks(sides = "l")+
  ylab("SN-38G [ug/l]")+
  xlab("Time [h]")

Any additional ideas?

Thanks again,

Mario