Plot optimisation

[AlexandreLeM]

Hello,

First of all, congratulations for your work.

I have started to use your package and for the moment I have coded the palbociclib model which gives similar residual concentration estimates to those estimated by monolix.

I would like to be able to make graphs with just the patient kinetics, observation and a confidence interval. The problem is that I can't remove the population kinetics from the plot(my_est).

I have created a code for a graph but perhaps there is an easier option.

I would also like to simulate scenarios with a new dosage and to make a graph with this new kinetic and its confidence interval. But I have not been able to find a way to do this. Could you tell me if this is possible and if so, show me?

Thanks!

Alexandre

---

Here is the code that allows me to obtain this type of graph

library(mapbayr)
library(mrgsolve)
library(MASS)
library(dplyr)
library(ggplot2)
library(scales)

######Model Palbociclib
code <- "
$PARAM @annotated
TVKA : 0.367 : Abs constant (h-1)
TVVC : 2710 : Central volume (L)
VP : 61300 : Peripheral volume(L)
TVCL : 60.2 : Clearance (L/h)
TVQ : 10.6 : Intercompartmental clearance (L/h)
ALAG : 0.658 : Lag time

ETA1 : 0 : VC (L)
ETA2 : 0 : Q (L)
ETA3 : 0 : CL (L/h)
ETA4 : 0 : Ka (h-1)

$PARAM @annotated @covariates
BW : 73.7 : Body Weight (kg)
AGE : 61 : AGE (years)

$OMEGA
0.091
1.59
0.131
0.699

$SIGMA
0.00 // proportional
0.09 // additive

$CMT @annotated
DEPOT : Depot compartment () [ADM]
CENT  : Central compartment (mg/L) [OBS]
PERIPH: Peripheral compartment ()

$TABLE
double DV = (CENT/VC)* exp(EPS(2)) ;
double PERI = (PERIPH/VP) ;
double gut = DEPOT ;

$MAIN
double BWVC = pow(BW/73.7,0.00906) ;
double AGECL = pow(AGE/61,-0.45) ;
double BWCL = pow(BW/73.7,0.484)  ;

double KA = TVKA * exp(ETA4 + ETA(4));
double VC = TVVC * exp(ETA1 + ETA(1))  * BWVC ;
double Q = TVQ  * exp(ETA2 + ETA(2))  ;
double CL = TVCL * exp(ETA3 + ETA(3)) * AGECL * BWCL ;

double K12 = Q/VC ;
double K21 = Q/VP ;
double KE = CL/VC ;

ALAG_DEPOT = ALAG ;

$ODE @annoted
dxdt_DEPOT = - DEPOT * KA ;
dxdt_CENT = DEPOT * KA + K21 * PERIPH - K12 * CENT - KE * CENT ;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH ;

$CAPTURE DV KA VC VP CL Q ALAG PERI gut
"

my_model <- mcode("Palbociclib_model", code)

######Variable
datfirstdose = as.POSIXct("2021-07-01 19:00:00")
datlastdose = as.POSIXct("2021-07-15 19:00:00")
datprel = as.POSIXct("2021-07-16 15:32:00")
Wt = 52
Age = 52

CPalbo = 97

tpsdose <- as.numeric(as.difftime(difftime(datlastdose,datfirstdose,units="hours"),units="hours"))
tpsprel <- as.numeric(as.difftime(difftime(datprel,datfirstdose,units="hours")))

dose <- 1000 * 125
interdoseinterval = 24

Debut = tpsdose
TpsSimu = 48

BorneInf =60
BorneSup = 62

mycol <- rgb(220,0,0, max = 255, alpha = 150)

######Estimation parameter
my_data <- my_model %>%
  adm_lines(time = 0, amt = dose, ii = interdoseinterval, addl = round(tpsdose/24)+1) %>%
  obs_lines((time = tpsprel), DV = CPalbo) %>%
  add_covariates(covariates = list(BW = Wt, AGE = Age)) %>%
  get_data()

my_est <-  mapbayest(my_model, my_data, verbose = F)

######Simulation for plot
My_simu <- my_est%>%
  augment(start = tpsdose, delta = 1, end = tpsdose+TpsSimu,ci =TRUE,ci_method ="delta",ci_width = 90)

simu <- My_simu$aug_tab
Conc_Pat <- My_simu$aug_tab$value[My_simu$aug_tab$type == "IPRED"]
time <- My_simu$aug_tab$time[My_simu$aug_tab$type == "IPRED"]-min(My_simu$aug_tab$time)
Ci_low <- My_simu$aug_tab$value_low[My_simu$aug_tab$type == "IPRED"]
Ci_up <- My_simu$aug_tab$value_up[My_simu$aug_tab$type == "IPRED"]

Data <- data.frame(time,Conc_Pat,Ci_low,Ci_up)
Data <- Data[Data$time!=(tpsprel-tpsdose),]

Cres <- round(min(Data$Conc_Pat),1)

ggplot(Data, aes(x=time,y=Conc_Pat)) +
  geom_ribbon(aes(ymin = Ci_low, ymax = Ci_up), fill = mycol) +
  geom_point(Data, mapping =  aes(x=tpsprel-tpsdose, y=CPalbo), size=4,color = "black") +
  geom_line(aes(y = Conc_Pat),size = 1,color="blue")+
  scale_x_continuous("Time (h)", limits = c(0, 48),breaks = c(0,12,24,36,48)) +
  scale_y_continuous("Concentrations (ug/L)") +   
  ggtitle(paste0("For a dose of ",dose/1000," mg, every ",interdoseinterval," h."),subtitle =paste0("Cres = ",Cres," ug/L")) +
  theme(panel.background = element_rect(fill = "white", colour = "black"))

[FelicienLL]

Hi Alexandre, thanks for posting here, Glad to here that it works on palbo data vs monolix

To answer your first question, there are no option to remove the population confidence interval. And it should be possible ideally. I'll patch this later. In the meantime, you can just set the value_up and value_low of "PRED" to NA in the aug_tab.

My_simu$aug_tab$value_low[My_simu$aug_tab$type=="PRED"] <- NA
My_simu$aug_tab$value_up[My_simu$aug_tab$type=="PRED"] <- NA

plot(My_simu)

Or just filter the data with IPRED to drop all the "PRED" related data

My_simu$aug_tab <- filter(My_simu$aug_tab, type == "IPRED")
plot(My_simu)

[FelicienLL]

Note that the reported confidence interval is an approximation to several respects:

• first, it is a Normal approximation of the posterior, which maybe is acceptable at the asymptote, but posterior is known to be non-normal especially if data is sparse (for instance see Meier et al, CPT PSP, 2020)
• this Normal approximation is based on the covariance matrix of estimation, calculated from the Hessian. Several methods exists to compute the Hessian. Default is stats::optimHess but one could try nlmixr::nlmixrHess (see https://github.com/FelicienLL/mapbayr/issues/43#issuecomment-1028151709). Both return results somewhat different from NONMEM: sometimes it is rather close (for simple problems), sometimes it is rather far-off (mainly for complex problems).
• The transposition of uncertainty from parameters to observables (i.e. concentrations) is also an approximation. A standard approach would be to simulate (say) 1,000 PK profiles and keep the 5th and 95th percentiles: this is can be achieve with (ci_method = "simulations"). Since it can be rather long, one can use the "delta method" (see Kümmel et al, CPT PSP, 2018), which is an approximation but faster. This is why I set it to default when the user want to plot the results.

For your second question, the best solution is use_posterior(). It turns your estimation object into a mrgsolve model with updated eta and covariates of your patient. Then feel free to simulate any dose regimen within the mrgsolve framework. You can also update_omega = TRUE to replace the OMEGA (prior) matrix into the posterior matrix, thus simulating with uncertainty.

I set a little example below from your code. It shows how to plot and report probability of success of several dosing regimens. Again: it assumes that parameter uncertainty is well described by the covariance matrix computed within the package, which, as you understand, is maybe (maybe not?) doubtful: I would recommend a comprehensive validation before implementing anything. Feel free to come back to me if you need anything else.

Hope it helps

Félicien

scenarios <- expand.ev(
  ID = 1:1000,
  amt = c(75000, 100000, 125000, 150000), 
  addl = 15, ii = 24) %>% 
  mutate(DOSE = amt/1000)

sims <- my_est %>% 
  use_posterior(update_omega = TRUE) %>% 
  ev(scenarios) %>% 
  mrgsim(start = 336, end = 360, obsonly = TRUE, carry_out = "DOSE")

tab <- sims %>% 
  group_by(DOSE, time) %>% 
  summarise(Q05 = quantile(DV, 0.05), 
            Q50 = quantile(DV, 0.50), 
            Q95 = quantile(DV, 0.95)) %>% 
  mutate(DOSE = as.factor(DOSE))
#> `summarise()` has grouped output by 'DOSE'. You can override using the
#> `.groups` argument.

tab %>% 
  ggplot(aes(time)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95, fill = DOSE), alpha = .5) + 
  geom_line(aes(y = Q50)) +
  facet_grid(.~DOSE) +
  theme_bw()

sims %>% 
  filter(time == 360) %>% 
  group_by(DOSE) %>% 
  summarise(
    MEDIAN_CONC = median(DV),
    OK = sum(DV > 97), 
    PERC_OK = scales::percent(OK/1000))
#> # A tibble: 4 x 4
#>    DOSE MEDIAN_CONC    OK PERC_OK
#>   <dbl>       <dbl> <int> <chr>  
#> 1    75        46.1     0 0%     
#> 2   100        61.4    33 3%     
#> 3   125        75.5   167 17%    
#> 4   150        90.0   385 38%

^{Created on 2022-03-30 by the reprex package (v2.0.1)}

[FelicienLL]

N.B.: note that I am currently optimizing the computational performance of the package: it is mainly in order to implement SIR. SIR will be more accurate than this Gaussian approximation when we want to derive a conf interval.