mapbayr is a free and open source package for maximum a posteriori bayesian estimation of PK parameters in R. Thanks to a single function, mapbayest()
, you can estimate individual PK parameters from:
It was designed to be easily wrapped in shiny apps in order to ease model-based Therapeutic Drug Monitoring, also referred to as Model-Informed Prediction Dosing (MIPD).
mapbayr is available on CRAN. You can install the development version from github by executing the following code in R console.
install.packages("devtools") devtools::install_github("FelicienLL/mapbayr")
mapbayr relies on mrgsolve for model implementation and ordinary differential equation solving which requires C++ compilers. If you are a Windows user, you would probably need to install Rtools. Please refer to the installation guide of mrgsolve for additional information.
library(mapbayr) library(mrgsolve)1) Properly code you model
code <- " $PARAM @annotated TVCL: 0.9 : Clearance TVV1: 10.0 : Central volume V2 : 10.0 : Peripheral volume of distribution Q : 1.0 : Intercompartmental clearance ETA1: 0 : Clearance (L/h) ETA2: 0 : Central volume (L) $PARAM @annotated @covariates BW : 70 : Body weight (kg) $OMEGA 0.3 0.3 $SIGMA 0.05 // proportional 0.1 // additive $CMT @annotated CENT : Central compartment (mg/L)[ADM, OBS] PERIPH: Peripheral compartment () $TABLE double DV = (CENT/V1) *(1 + EPS(1)) + EPS(2); $MAIN double CL = TVCL * exp(ETA1 + ETA(1)) * pow(BW / 70, 1.2) ; double V1 = TVV1 * exp(ETA2 + ETA(2)) ; double K12 = Q / V1 ; double K21 = Q / V2 ; double K10 = CL / V1 ; $ODE dxdt_CENT = K21 * PERIPH - (K10 + K12) * CENT ; dxdt_PERIPH = K12 * CENT - K21 * PERIPH ; $CAPTURE DV CL " my_model <- mcode("Example_model", code)
my_data <- data.frame(ID = 1, time = c(0,6,15,24), evid = c(1, rep(0,3)), cmt = 1, amt = c(100, rep(0,3)), rate = c(20, rep(0,3)), DV = c(NA, 3.9, 1.1, 2), mdv = c(1,0,0,1), BW = 90) my_data #> ID time evid cmt amt rate DV mdv BW #> 1 1 0 1 1 100 20 NA 1 90 #> 2 1 6 0 1 0 0 3.9 0 90 #> 3 1 15 0 1 0 0 1.1 0 90 #> 4 1 24 0 1 0 0 2.0 1 90
my_est <- mapbayest(my_model, data = my_data)
As building dataset into a NM-TRAN format can be painful, you can use pipe-friendly obs_rows()
, adm_rows()
and add_covariates()
functions in order to pass administration and observation information, and perform the estimation subsequently.
my_est <- my_model %>% adm_rows(time = 0, amt = 100, rate = 20) %>% obs_rows(time = 6, DV = 3.9) %>% obs_rows(time = 20, DV = 1.1) %>% obs_rows(time = 24, DV = 2, mdv = 1) %>% add_covariates(BW = 90) %>% mapbayest()4) Then, use the estimations
The results are returned in a single object (“mapbayests” S3 class) which includes input (model and data), output (etas and tables) and internal arguments passed to the internal algorithm (useful for debugging). Additional methods are provided to ease visualization and computation of a posteriori outcomes of interest.
print(my_est) #> Model: Example_model #> ID : 1 individual(s). #> OBS: 2 observation(s). #> ETA: 2 parameter(s) to estimate. #> #> Estimates: #> ID ETA1 ETA2 #> 1 1 0.3872104 0.1569604 #> #> Output (4 lines): #> ID time evid cmt amt rate mdv DV IPRED PRED CL BW ETA1 ETA2 #> 1 1 0 1 1 100 20 1 NA 0.000 0.000 1.79 90 0.387 0.157 #> 2 1 6 0 1 0 0 0 3.9 4.162 5.174 1.79 90 0.387 0.157 #> 3 1 15 0 1 0 0 0 1.1 1.087 1.647 1.79 90 0.387 0.157 #> 4 1 24 0 1 0 0 1 2.0 0.556 0.959 1.79 90 0.387 0.157
# Easily extract a posteriori parameter values to compute outcomes of interest get_eta(my_est) #> ETA1 ETA2 #> 0.3872104 0.1569604 get_param(my_est, "CL") #> [1] 1.79217 # The `use_posterior()` functions updates the model object with posterior values and covariates to simulate like with a regular mrgsolve model my_est %>% use_posterior() %>% data_set(expand.ev(amt = c(50, 100, 200, 500), dur = c(5, 24)) %>% mutate(rate = amt/dur)) %>% carry_out(dur) %>% mrgsim() %>% plot(DV~time|factor(dur), scales = "same")
mapbayr is under development. Your feedback for additional feature requests or bug reporting is welcome. Contact us through the issue tracker.
mapbayr is a generalization of the “MAP Bayes estimation” tutorial available on the mrgsolve blog. Additional features are:
mapbayest()
.Reliability of parameter estimation against NONMEM was assessed for a wide variety of models and data. The results of this validation study were published in CPT:Pharmacometrics & System Pharmacology, and materials are available in a dedicated repository. If you observe some discrepancies between mapbayr and NONMEM on your own model and data, feel free to contact us through the issue tracker.
mrgsolve model specificationmapbayr contains a library of example model files (.cpp), accessible with exmodel()
. You are invited to perform MAP-Bayesian estimation with your own models. These model files should be slightly modified in order to be “read” by mapbayr with the subsequent specifications:
$PARAM @annotated ETA1 : 0 : CL (L/h) ETA2 : 0 : VC (L) ETA3 : 0 : F () //do not write ETA(1) //do not write iETA
@covariates
tag to record covariates in the $PARAM
block. Otherwise, you will not be allowed to pass a dataset with covariates columns.$PARAM @annotated @covariates BW : 70 : Body weight (kg) SEX : 0 : Sex (0=Male, 1=Female)
@annotated
tag must be used to record compartments.obs_rows()
to fill the ‘cmt’ column in your dataset.mapbayest()
function. The information is mandatory if you use adm_rows()
to build your dataset in order to automatically set the value of the ‘cmt’ column. Especially useful if you use a model with an absorption from several depot compartment requiring to duplicate administrations lines in the data set.//example: model with dual zero and first order absorption in compartment 1 & 2, respectively, and observation of parent drug + metabolite $CMT @annotated DEPOT: Depot [ADM] CENT_PAR: examplinib central [ADM, OBS] PERIPH : examplinib peripheral CENT_MET : methylexamplinib central [OBS]
$PARAM
.$PARAM
. This cannot be checked by mapbayr !$OMEGA 0.123 0.456 0.789 $OMEGA @block 0.111 0.222 0.333 // reminder: omega values can be recorded in multiple $OMEGA blocks
The definition of the $SIGMA
block may not be as straightforward as other blocks, but we tried to keep it as simple as possible. Keep in mind that mapbayr always expect a pair of sigma values for each type of dependent variable: the first value for proportional error, the second for additive.
Two situations can be distinguished:
$CMT
.Simply write one pair of sigma values to describe proportional and additive error on your concentrations. This error model will be automatically applied to the compartment where observations were recorded in your dataset (i.e. value of CMT when MDV = 0).
$SIGMA 0.111 0 // proportional error
$SIGMA 0 0.222 // (log) additive error
$SIGMA 0.333 0.444 // mixed error
$CMT
.Write as many pairs of sigma values as there are compartments assigned with [OBS] in $CMT
. The order of the pair must respect the order in which compartments were assigned. To put it more clearly, the sigma matrix will be interpreted as such whatever the model :
//example: correlated proportional error between parent and metabolite $SIGMA @block 0.050 // proportional error on parent drug 0.000 0.000 // additive error on parent drug 0.100 0.000 0.200 // proportional error on metabolite 0.000 0.000 0.000 0.000 // additive error on metabolite // reminder: sigma values can be recorded in multiple $SIGMA blocks6.
$TABLE
block or $ERROR
block
DV
. Mind the code, especially if concentrations are observed in multiple compartments.$TABLE double DV = (CENTRAL / VC) * exp(EPS(2)) ;
$TABLE double PAR = (CENT_PAR / V) * (1 + EPS(1)) ; double MET = (CENT_MET / V) * (1 + EPS(3)) ; double DV = PAR ; if(self.cmt == 4) DV = MET ; // reminder: use "self.cmt" to internaly refer to a compartment in a mrgsolve model code.
Note that mapbayr does not strictly rely on this $ERROR
block to define the residual error internally and compute the objective function value, but on information passed in the $SIGMA
block. However, we strongly advise you to properly code your $ERROR
block with EPS(1)
, EPS(2)
etc…, if only to use your code as a regular mrgsolve model code and simulate random effects.
$PARAM
and $OMEGA
. This cannot be checked by mapbayr !$PK double CL = TVCL * exp(ETA1 + ETA(1)) ;
mapbayest()
anyway)mapbayest()
anyway)RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4