Title:	MAP-Bayesian Estimation of PK Parameters
Version:	0.10.0
Description:	Performs maximum a posteriori Bayesian estimation of individual pharmacokinetic parameters from a model defined in 'mrgsolve', typically for model-based therapeutic drug monitoring. Internally computes an objective function value from model and data, performs optimization and returns predictions in a convenient format. The performance of the package was described by Le Louedec et al (2021) <doi:10.1002/psp4.12689>.
License:	GPL-3
Encoding:	UTF-8
RoxygenNote:	7.2.3
Imports:	dplyr, ggplot2 (≥ 3.4.0), magrittr, mrgsolve (≥ 1.0.8), progress, purrr, rlang, stringr, tidyr
URL:	https://github.com/FelicienLL/mapbayr
BugReports:	https://github.com/FelicienLL/mapbayr/issues
Suggests:	knitr, lubridate, minqa, scales, testthat (≥ 3.0.0), tibble
Depends:	R (≥ 2.10)
LazyData:	true
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2023-07-17 07:54:03 UTC; le_louedec
Author:	Felicien Le Louedec [aut, cre], Kyle T Baron [ctb], Laura Morvan [ctb]
Maintainer:	Felicien Le Louedec <felicienlelouedec@proton.me>
Repository:	CRAN
Date/Publication:	2023-07-17 08:20:02 UTC

Add covariate columns to a dataset

Description

The add_covariates() function adds an one or several covariate columns to a dataset provided as a proper data.frame or within a 'mrgsolve' model. Used in combination with adm_rows() and obs_rows(), it helps the creation of datasets in the proper format for simulations with 'mrgsolve' or parameter estimation with 'mapbayr', as explained in data_helpers.

Usage

add_covariates(x, ...)

## S3 method for class 'data.frame'
add_covariates(x, ..., covariates = list(), AOLA = FALSE, TOLA = FALSE)

## S3 method for class 'mrgmod'
add_covariates(x, ..., covariates = list(), AOLA = NULL, TOLA = NULL)

Arguments

Value

a data.frame, or a 'mrgsolve' model with a dataset in the ⁠@args$data⁠ slot (accessible with get_data()).

See Also

data_helpers

Examples

# Cannot start from scratch
## Not run: 
add_covariates(BW = 90, SEX = 0)

## End(Not run)

library(magrittr)
adm_rows(time = c(0, 24, 48), cmt = 1, amt = c(100, 200, 300)) %>%
  add_covariates(BW = c(90, 85, 80), SEX = 0)

# If covariates are stored in a list, use `covariates = `
adm_rows(time = c(0, 24, 48), cmt = 1, amt = c(100, 200, 300)) %>%
  add_covariates(covariates = list(BW = c(90, 85, 80), SEX = 0))

# Missing values are filled with the "next observation carried backward" rule
adm_rows(time = c(0, 24, 48), cmt = 1, amt = c(100, 200, 300)) %>%
  add_covariates(BW = c(90, 85, 80), SEX = 0) %>%
  obs_rows(time = 36, DV = .0123, cmt = 2)
# Always verify the output in case of time-varying covariates

# Possibility to add Time and Amount of last administration as covariates
adm_rows(time = c(0, 24, 48), amt = c(100, 200, 300), cmt = 1) %>%
  obs_rows(time = c(8, 16, 32, 40), cmt = 2, DV = runif(4)) %>%
  add_covariates(TOLA = TRUE, AOLA = TRUE) %>%
  obs_rows(time = 72, cmt = 2, DV = .123) # AOLA/TOLA re-updated afterwards

# Automatic inclusion of `TOLA`/`AOLA` if they are covariates of the model
library(mrgsolve)
model <- mcode("model", "
$PARAM @annotated @covariates
TOLA : 0 : Time Last Adm
AOLA : 0 : Amount Last Adm
", compile = FALSE)

model %>%
  adm_rows(time = c(0, 24, 48), amt = c(100, 200, 300), cmt = 1) %>%
  add_covariates() %>%
  get_data()

Add administration lines to a dataset

Description

The adm_rows() function adds an one or several administration lines to a dataset provided as a proper data.frame or within a 'mrgsolve' model. Used in combination with obs_rows() and add_covariates(), it helps the creation of datasets in the proper format for simulations with 'mrgsolve' or parameter estimation with 'mapbayr', as explained in data_helpers.

Usage

adm_rows(x, ...)

## S3 method for class 'data.frame'
adm_rows(
  x,
  ID = NULL,
  time = NULL,
  evid = 1L,
  cmt,
  amt,
  mdv = 1L,
  addl = NULL,
  ss = NULL,
  ii = NULL,
  rate = NULL,
  .datehour = NULL,
  ...
)

## S3 method for class 'missing'
adm_rows(...)

## S3 method for class 'mrgmod'
adm_rows(x, cmt = adm_cmt(x), rate = NULL, ...)

Arguments

Value

a data.frame, or a 'mrgsolve' model with a dataset in the ⁠@args$data⁠ slot (accessible with get_data()).

See Also

data_helpers

Examples

# Create a dataset from scratch
adm_rows(amt = 100, cmt = 1)

# Pipe-friendly addition of administration record to a pre-existing dataset
library(magrittr)
adm_rows(amt = 100, cmt = 1) %>%
  adm_rows(time = 3, amt = 200, cmt = 1, addl = 3, ii = 1)

# Inform times using the `.datehour` argument:
adm_rows(.datehour = "2020-01-01 11:11", amt = 100, cmt = 1) %>%
  adm_rows(.datehour = "2020-01-02 22:22", amt = 200, cmt = 1) %>%
  adm_rows(time = 48, amt = 300, cmt = 1)

# Start from a 'mrgsolve' model
library(mrgsolve)
house() %>%
  adm_rows(amt = 100, cmt = 1) %>%
  adm_rows(time = 3, amt = 200, cmt = 1, addl = 3, ii = 1) %>%
  mrgsim(delta = 1)

# Default administration compartments
# Set default administration compartments in the code with `[ADM]`
model <- mcode("model", "
$CMT @annotated
DEPOT : Depot [ADM]
CENTR : Central
", compile = FALSE)
adm_cmt(model)

# Thus, no need to manually specify `cmt = 1` anymore.
model %>%
  adm_rows(amt = 100) %>%
  adm_rows(time = 3, amt = 200, addl = 3, ii = 1) %>%
  get_data()

# Automatic lines duplication if multiple depot compartments
# Automatic `rate = -2` increment if model with 0-order absorption
model <- mcode("model", "
$PARAM DUR = 1.0
$CMT @annotated
DEPOT : Depot [ADM]
CENTR : Central [ADM]
$MAIN
D_CENTR = DUR ;
", compile = FALSE)
adm_cmt(model)

model %>%
  adm_rows(amt = 100) %>%
  adm_rows(time = 3, amt = 200, addl = 3, ii = 1) %>%
  get_data()

Return the mapbay_tab as a data.frame

Description

Return the mapbay_tab as a data.frame

Usage

## S3 method for class 'mapbayests'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

a data.frame (the mapbay_tab from estimation)

Compute full PK profile prediction from mapbayr estimates.

Description

Compute full PK profile prediction from mapbayr estimates.

Usage

augment(x, ...)

Arguments

Value

an augmented object (depending on the object passed).

Compute full PK profile prediction from mapbayr estimates.

Description

Compute full PK profile prediction from mapbayr estimates.

Usage

## S3 method for class 'mapbayests'
augment(
  x,
  data = NULL,
  start = NULL,
  end = NULL,
  delta = NULL,
  ci = FALSE,
  ci_width = 90,
  ci_method = "delta",
  ci_sims = 500,
  ...
)

Arguments

Details

This function is called in the background by plot() in order to simulate the full PK profile, and return a mapbayests object with an additional aug_tab data.frame inside. The latter is used with by the plot method. The time grid, for each PK profile (i.e. patient) is defaulted with the minimum time in the dataset for start and the maximum time in the dataset +20% for end. delta is a power of 10 (e.g. 0.1, 1, 10 etc...), automatically chosen to render visually appealing graphs with a reasonable computing time (about 200 time points). Additional arguments can be passed to mrgsim() through .... Note that recsort is set to 3 (see mrgsolve documentation for more details).

Value

a mapbayests object, augmented of an aug_tab data.frame.

Examples

#x is the result of `mapbayest()`.
#Default plot is returned by:
# plot(x)
#Argument passed to `plot()` are passed to `augment()` in the background:
# plot(x, end = 240, ci = TRUE)
#Save the augmented object if simulation time is long
# x2 <- augment(x, ci = TRUE, ci_method = "simulations", ci_sims = 10000) %>%
# plot(x2)

Check if model is valid for 'mapbayr'

Description

Checks that the model respects points related exclusively to 'mapbayr'. Useful at the time you wish to convert a "regular" 'mrgsolve' model you used for simulation into a model to perform MAP-Bayesian estimation. Note that some elements cannot be checked:

• In ⁠$MAIN⁠ block, make sure that you added ⁠ETA1, ETA2...⁠ in the code. For instance: ⁠double CL = TVCL * exp(ETA(1) + ETA1) ;⁠.

• In ⁠$OMEGA⁠ block, make sure the order of the (diagonal) values is the same as for ETAs in ⁠$PARAM⁠. For instance, if ETA1 corresponds to clearance, the first value in ⁠$OMEGA⁠ must be the variance of clearance.

• In ⁠$SIGMA⁠ block, make sure the order is respected: proportional error first, and additive error secondly.

Usage

check_mapbayr_model(x, check_compile = TRUE)

Arguments

Value

TRUE (invisibly) if checks are passed, errors otherwise.

Examples

library(mapbayr)
library(mrgsolve)
## Not run: check_mapbayr_model(house())

Compute the objective function value

Description

Compute the objective function value

Usage

compute_ofv(
  eta,
  qmod,
  sigma,
  omega_inv,
  all_cmt,
  log_transformation,
  lambda = 1,
  idvaliddata,
  idDV,
  idcmt,
  idblq = NULL,
  idlloq = NULL,
  ...
)

do_compute_ofv(eta, argofv, ...)

Arguments

Details

This function is called iteratively by the optimization function. Arguments should not be passed directly, but generated by the pre-processing functions (see preprocess.ofv).

Value

a single numeric value (the objective function value)

Data helpers: functions to build the dataset

Description

Use adm_rows(), obs_rows() and add_covariates() to create or modify a dataset from scratch, from a pre-existing dataset, or from a dataset stored into a 'mrgsolve' model.

Details

Instead of importing a '.csv' file, or painfully build a data set with a call to data.frame() and mind how to format the data, you can pass information about:

• administrations with adm_rows(),

• observations with obs_rows(),

• covariates with add_covariates(),

all being called jointly with a pipe (⁠%>%⁠ or ⁠|>⁠). These functions can be used to create or modify a dataset as a proper data.frame, or to create or modify a dataset within a 'mrgsolve' model (⁠@args$data⁠ slot). The latter is particularly useful in order to:

• automatically use default administration and observation compartments,

• automatically duplicate rows if there are several depot compartments,

• automatically set rate = -2 if model has zero-order absorption pathways,

• automatically duplicate rows if concentrations of Parent drug and Metabolite are observed together,

• automatically add "Amount Of Last Administration" and "Time Of Last Administration" variables if these are covariates,

• subsequently call mrgsim() or mapbayest().

Examples

library(magrittr)
# First option: work with a data.frame

adm_rows(amt = 1000, cmt = 1, addl = 4, ii = 12) %>%
  obs_rows(time = c(12.3, 45.6), DV = c(.111, .222), cmt = 2) %>%
  obs_rows(time = 48, cmt = 2) %>%
  add_covariates(BW = 90, SEX = 0, TOLA = TRUE)

# You can even inform "time" using date and hours:
adm_rows(amt = 1000, cmt = 1, addl = 4, ii = 12, .datehour = "2022-01-01 11:11:11") %>%
  obs_rows(.datehour = "2022-01-02 22:22:22", DV = 0.111, cmt = 2)

# Second option: work with a dataset within a 'mrgsolve' model

mod <- exmodel(add_exdata = FALSE)
# call `mrgsolve::see(mod)` to see how default compartment were coded
adm_cmt(mod)
obs_cmt(mod)

mod %>%
  adm_rows(amt = 10000) %>%
  obs_rows(time = c(1.5, 4.4, 7.5, 24.6), DV = c(91.2904, 110.826, 79.384, 20.6671)) %>%
  # get_data() # for curiosity, you can extract the data set at this step
  mapbayest()

Deprecated functions

Description

Deprecated functions

Usage

mbrest(...)

adm_lines(...)

obs_lines(...)

Arguments

Details

• mbrest() is now mapbayest()

• adm_lines() is now adm_rows()

• obs_lines() is now obs_rows()

Simulate with mapbayr

Description

A wrapper around mrgsolve::mrgsim() for results generated from mapbayest(). Exported for the purpose of utility but might be prone to changes.

Usage

do_mapbayr_sim(
  x,
  data,
  recsort = 3,
  output = "df",
  ...,
  eta = NULL,
  nrep = NULL,
  new_omega = NULL,
  new_sigma = NULL
)

Arguments

Value

An output from mrgsolve::mrgsim().

Examples

library(mrgsolve)
mod <- exmodel(1, exdata = FALSE)
dat <- exdata(ID = c(1,2))

# Classic framework
set.seed(123)
do_mapbayr_sim(x = mod, data = dat, Request = "DV")

# No random effect
do_mapbayr_sim(x = zero_re(mod), data = dat)
do_mapbayr_sim(x = mod, data = dat, new_omega = "zero_re")

# New random effects
## New omega matrix
do_mapbayr_sim(x = mod, data = dat, new_omega = dmat(0.1, 0.03, 0.01), nrep = 10)

## Matrix with "eta" as mean and "new_omega" as variance covariance matrix
etamat <- get_eta(est001, output = "num")[1:2,]

do_mapbayr_sim(
  x = mod, data = dat, nrep = 10,
  eta = etamat, new_omega = dmat(0.1, 0.03, 0.01)
)

Estimation object

Description

An example of mapbayests object, corresponding to the parameter estimation of the 8 subjects from model 1. Note that the model object within is not associated to a shared object, which make some features unavailable. This object can be re-generated by executing est001 <- mapbayest(exmodel(ID = 1:8))

Usage

est001

Format

An object of class mapbayests of length 9.

See Also

mapbayest

Generate a vector of "ETA"

Description

Generate a vector of "ETA" values. If x is a mrgsolve model, these will be extracted from values defined in ⁠$PARAM⁠. Otherwise, any numeric values passed to x and ... as vector(s) or list(s) will be coerced as a single vector. Alternatively, if x and ... are missing, generate a vector of ETA equal to val of length n.

Usage

eta(x, ..., n, val = 0)

Arguments

Value

a single named vector of numeric

Examples

# Extract ETA from the model
mod <- exmodel()
eta(mod)

# Coerce numeric values
eta(0.1, 0.2, c(0.3, 0.4), list(0.5, 0.6))
eta(rnorm(4))

# Generate a sequence from scratch
eta(n = 3)
eta(n = 3, val = 0.001)

Example model and data

Description

A collection of example models and corresponding data to test and explore mapbayr.

Usage

exmodel(
  num = 1,
  add_exdata = TRUE,
  cache = TRUE,
  quiet = getOption("mrgsolve_mread_quiet", TRUE),
  ...,
  ID = 1,
  clean_data = TRUE
)

exdata(num = 1, ID = 1, clean_data = TRUE)

Arguments

Details

Available models are:

• 1: Base model. A simple monocompartmental PK model with inter-individual variability on absorption constant (KA), volume of distribution (VC) and clearance (CL). The residual error model is proportional.

• 6: Complex absorption model. Dual 0- and 1st orders absorption phenomenons.

• 301: Time-varying covariates. A continuous covariate (body weight "BW") and a categorical one (sex "SEX") influence the clearance parameter. In the corresponding dataset, the values randomly changes from one record to another within a single individual.

• 401: Metabolite. The PK model of both a parent drug and its metabolite.

An example dataset of eight (simulated) individuals is available for each model. Individuals differ in terms of sampling times (sparse or rich) and dosing regimen (single or multiple dosing).

Model code and data files are stored at the location given by system.file("exmodel", package = "mapbayr").

These models and data were created for the validation study of mapbayr published in CPT:Pharmacometrics & System Pharmacology. More models and full datasets can be accessed in a dedicated repository

Value

exmodel() reads and compiles code, and returns a (mrgmod) model object. exdata() returns a data.frame.

Source

https://github.com/FelicienLL/mapbayr-CPTPSP-2021

Examples

# Models can be loaded with data (the default), ready for parameter estimation
est <- mapbayest(exmodel())

# Number of subjects in dataset can be chosen up to 8 individuals
exdata(301, ID = c(5,8))

Filter a dataset within a mrgmod

Description

Filter a dataset within a mrgmod

Usage

## S3 method for class 'mrgmod'
filter(.data, ..., .preserve = FALSE)

Arguments

Value

a mrgmod

Examples

library(magrittr)
mod <- mrgsolve::mcode("mod", "$CMT FOO", compile = FALSE)
mod %>%
  adm_rows(amt = c(100, 200, 300), cmt = 1) %>%
  filter(amt != 200) %>%
  get_data()

Get content from object

Description

Helpful functions to get content from a mrgmod object (i.e. data) or from a mapbayests object (data, eta, cov, param, phi).

Usage

get_data(x, ...)

## S3 method for class 'mrgmod'
get_data(x, ...)

## S3 method for class 'mapbayests'
get_data(x, ..., output = "df")

get_eta(x, ...)

## S3 method for class 'mapbayests'
get_eta(x, ..., output = NULL)

get_cov(x, ...)

## S3 method for class 'mapbayests'
get_cov(x, ..., simplify = TRUE)

get_param(x, ...)

## S3 method for class 'mapbayests'
get_param(x, ..., output = NULL, keep_ID = NULL, keep_names = NULL)

get_phi(x, ...)

## S3 method for class 'mapbayests'
get_phi(x, ...)

Arguments

Value

the class of the object returned depends on the function, and on their arguments. Typically, a data.frame or a vector if the output can be reduced to one line.

Examples

# From a model object (mrgmod)
mod <- exmodel(ID = 1:2, cache = FALSE, capture = "CL")
get_data(mod)

# From an estimation object (mapbayests)
est <- mapbayest(mod)
get_data(est)
get_data(est, output = "list")

get_eta(est)
get_eta(est, output = "list")

get_cov(est)

get_param(est)

get_phi(est)

Plot posterior distribution of bayesian estimates

Description

Plot posterior distribution of bayesian estimates

Usage

## S3 method for class 'mapbayests'
hist(x, select_eta = x$arg.optim$select_eta, shk = c("sd", "var", NA), ...)

Arguments

Details

Use this function to plot the results of the estimations, in the form of histograms with the a priori distribution in the background. For every parameter, the inter-individual variability is displayed, as well as the percentile of the patient in the corresponding distribution (if n = 1 patient). For additional modifications, you can add extra ⁠+function(...)⁠ in order to modify the plot as a regular ggplot2 object.

Value

a ggplot object.

Examples

est <- mapbayest(exmodel(ID = 1))

# Default Method
h <- hist(est)

# Can be modified with `ggplot2`
h +
  ggplot2::labs(title = "Awesome estimations")

# Select the ETAs
hist(est, select_eta = c(1,3))

Estimate parameters (maximum a posteriori)

Description

The main function of the mapbayr package. Performs a maximum a posteriori Bayesian estimation of parameters, from a mrgsolve model object and a dataset containing information about administrations and observed concentrations.

Usage

mapbayest(
  x,
  data = NULL,
  method = c("L-BFGS-B", "newuoa"),
  hessian = stats::optimHess,
  select_eta = NULL,
  lambda = 1,
  lloq = NULL,
  force_initial_eta = NULL,
  quantile_bound = 0.001,
  control = list(),
  check = TRUE,
  verbose = TRUE,
  progress = TRUE,
  reset = 50,
  output = NULL,
  ...
)

Arguments

Value

a mapbayests object. Basically a list containing:

• model: the model object

• arg.ofv.optim, arg.ofv.fix, arg.ofv.id: arguments passed to the optimization function. Useful for debugging but not relevant for a basic usage. Access to the data with get_data(x)

• opt.value: the original output of the optimization function

• final_eta: a list of individual vectors of final estimates. Access it with x$final_eta or get_eta(x).

• covariance: a list of individual variance-covariance matrix of estimation. Access it with x$covariance or get_cov(x).

• mapbay_tab: an output table containing the results of your estimations (data, IPRED, PRED, covariates, captured items, ETA etc...). Access it with x$mapbay_tab, as.data.frame(x) or as_tibble(x).

• information: run times and package versions.

See Also

hist.mapbayests

plot.mapbayests

use_posterior

Examples

# First, code a model
code1 <- "$PARAM ETA1 = 0, ETA2 = 0,
KA = 0.5, TVCL = 1.1, TVV = 23.3
$OMEGA 0.41 0.32
$SIGMA 0.04 0
$CMT DEPOT CENT
$PK
double CL=TVCL*exp(ETA1+ETA(1));
double V=TVV*exp(ETA2+ETA(2)) ;
$ERROR
double DV=CENT/V*(1+EPS(1))+EPS(2);
$PKMODEL ncmt = 1, depot = TRUE
$CAPTURE DV CL
"

my_model <- mrgsolve::mcode("my_model", code1)
# Then, import your data
my_data <- data.frame(ID = 1, TIME = c(0, 1.1, 5.2, 12.3), EVID = c(1,0,0,0), AMT = c(500, 0,0,0),
 CMT = c(1,2,2,2), DV = c(0, 15.1, 29.5, 22.3))
print(my_data)

# And estimate
my_est <- mapbayest(x = my_model, data = my_data)
print(my_est)
# see also plot(my_est) and hist(my_est)

# Use your estimation
get_eta(my_est)
get_param(my_est)
as.data.frame(my_est)
use_posterior(my_est)

Make mapbayr plot

Description

Make mapbayr plot

Usage

mapbayr_plot(
  aug_tab,
  obs_tab = NULL,
  PREDICTION = c("IPRED", "PRED"),
  MODEL_color = NULL
)

Arguments

Value

a ggplot object

Examples

aug <- data.frame(
  ID = 1, name = factor("DV"), cmt = 2, time = rep(c(0,8,16,24), each = 2),
  type = rep(c("PRED", "IPRED"), 4), value = c(0,0, 1, 2, 4, 8, 2, 4)
  )

obs <- data.frame(
  ID = 1, time = c(6, 20), evid = 0,
  mdv = c(0,1), DV = c(0.5, 5), cmt = 2
  )

mapbayr_plot(aug, obs)
mapbayr_plot(aug, obs, PREDICTION = "IPRED")

aug2 <- dplyr::bind_rows(
  FOO = aug,
  BAZ = dplyr::mutate(aug, value = value * 2),
  BAR = dplyr::mutate(aug, value = value * 3),
  .id = "MODEL"
  )

mapbayr_plot(aug2, obs)
mapbayr_plot(aug2, obs, PREDICTION = "IPRED")
mapbayr_plot(aug2, obs, PREDICTION = "IPRED", MODEL_color = c(FOO = "black"))

Visual Predicted Checks

Description

Visual Predicted Checks

Usage

mapbayr_vpc(
  x,
  data = NULL,
  nrep = 500,
  pcvpc = TRUE,
  idv = "time",
  stratify_on = NULL,
  start = NULL,
  end = NULL,
  delta = 1,
  ...
)

Arguments

Details

• Prediction-corrected VPC

By default, VPC are prediction corrected (Bergstrand et al (2011) doi:10.1208/s12248-011-9255-z). This correction is advised if several levels of doses or covariates are in the dataset for instance. Note that the implemented correction formula does not take into account the 'lower bound' term (lbij), nor the log-transformed variables.

Value

a ggplot object, results of the VPC. The median and the 50%, 80% and 90% prediction intervals of the simulated distributions are reported.

Examples

library(mrgsolve)
library(magrittr)
# Define a model. Adding variability to house model because default is 0.
mod <- house() %>%
  omat(dmat(rep(0.2,4)))

# Creating dataset for the example
# Same concentration, but different dose (ID 2) and covariate (ID 3)
data <- adm_rows(ID = 1, amt = 1000, cmt = 1, addl = 6, ii = 12) %>%
  obs_rows(DV = 50, cmt = 2, time = 7 * 12) %>%
  adm_rows(ID = 2, time = 0, amt = 2000, cmt = 1, addl = 6, ii = 12) %>%
  obs_rows(DV = 50, cmt = 2, time = 7 * 12) %>%
  adm_rows(ID = 3, time = 0, amt = 1000, cmt = 1, addl = 6, ii = 12) %>%
  obs_rows(DV = 50, cmt = 2, time = 7 * 12) %>%
  add_covariates(SEX = c(0,0,0,0,1,1))

mapbayr_vpc(mod, data, nrep = 30) # prediction-corrected by default
mapbayr_vpc(mod, data, idv = "tad", start = 72, nrep = 30)
mapbayr_vpc(mod, data, pcvpc = FALSE, nrep = 30)
mapbayr_vpc(mod, data, stratify_on = "SEX", nrep = 30)

Average predictions from multiple models

Description

Model Averaging consists in analyzing the same data with different models and to average their predictions. In order to perform weighted means of clearance predictions, (or concentrations, or any metric of interest), it is necessary to compute the "weight" of each estimation. It is informed by the likelihood of estimation. Two weighting scheme are currently implemented, one based on the log- likelihood ("LL", the default), the other on the Akaike criterion ("AIC"). The method was previously described by Uster et al (2021) doi:10.1002/cpt.2065.

Usage

model_averaging(
  ...,
  output_function = as.data.frame,
  scheme = c("LL", "AIC"),
  estlist = NULL
)

compute_weights(..., scheme = c("LL", "AIC"), estlist = NULL)

do_model_averaging(list_of_tabs, weights_matrix)

Arguments

Value

• model_averaging() and do_model_averaging(): a data.frame of the same dimensions and attributes as the outputs

• compute_weights(): a matrix with IDs as rows and estimation weights as columns

Examples

library(magrittr)

# Three different models: A, B, and C.
modA <- exmodel(1, add_exdata = FALSE)
modB <- mrgsolve::param(modA, TVCL = 2, TVVC = 30)
modC <- mrgsolve::param(modA, TVCL = 10)

# A common dataset that has 2 patients (ID 2 & 9)
data <- adm_rows(ID = 2, time = 0, amt = 200, addl = 3, ii = 24, cmt = 1) %>%
    obs_rows(ID = 2, time = 84, DV = 1.5, cmt = 2) %>%
    adm_rows(ID = 9, time = 0, amt = 100, addl = 3, ii = 24, cmt = 1) %>%
    obs_rows(ID = 9, time = 96, DV = 1, cmt = 2)

# Three different estimation objects: A, B and C.
estA <- mapbayest(modA, data)
as.data.frame(estA)
plot(estA) # Fit is pretty good

estB <- mapbayest(modB, data)
as.data.frame(estB)
plot(estB) # Excellent fit

estC <- mapbayest(modC, data)
as.data.frame(estC)
plot(estC) # Fit is worst

# Model averaging
model_averaging(A = estA, B = estB, C = estC)
# Weighted average of the table returned by as.data.frame(est))

# Internally, it first computes the "weight" of each model such as:
W <- compute_weights(A = estA, B = estB, C = estC)

# Then multiply the prediction table with each weight such as:
do_model_averaging(
  list_of_tabs = list(
    A = as.data.frame(estA),
    B = as.data.frame(estB),
    C = as.data.frame(estC)
    ),
  weights_matrix = W
  )

# If you do not want to perform an average of the full table, you can specify
# a function that takes the estimation object as an input and returns
# value(s) of interest: a single prediction, a clearance value, a full
# table of augmented predictions... as long as the structure of the final
# object is the same whatever the model.

reframe <- function(est){
  # From any estimation object, return a table with ID, time and predictions
  as.data.frame(est)[,c("ID", "time", "DV", "IPRED")]
}

model_averaging(A = estA, B = estB, C = estC, output_function = reframe)

# Make a plot that compares predictions
List_aug_tab <- lapply(
  X = list(A = estA, B = estB, C = estC),
  FUN = \(x) augment(x)$aug_tab
)
List_aug_tab$.AVERAGE <- do_model_averaging(List_aug_tab, W)

mapbayr_plot(
 aug_tab = dplyr::bind_rows(List_aug_tab, .id = "MODEL"),
 obs_tab = data,
 PREDICTION = "IPRED",
 MODEL_color = c(.AVERAGE = "black")
)

Add observation lines to a dataset

Description

The obs_rows() function adds an one or several observation lines to a dataset provided as a proper data.frame or within a 'mrgsolve' model. Used in combination with adm_rows() and add_covariates(), it helps the creation of datasets in the proper format for simulations with 'mrgsolve' or parameter estimation with 'mapbayr', as explained in data_helpers.

Usage

obs_rows(x, ...)

## S3 method for class 'data.frame'
obs_rows(
  x,
  ID = NULL,
  time = NULL,
  evid = 0L,
  cmt,
  DV = NA_real_,
  mdv = NULL,
  .datehour = NULL,
  ...
)

## S3 method for class 'missing'
obs_rows(...)

## S3 method for class 'mrgmod'
obs_rows(x, cmt = NULL, DV = NA_real_, DVmet = NULL, ...)

Arguments

Value

a data.frame, or a 'mrgsolve' model with a dataset in the ⁠@args$data⁠ slot (accessible with get_data()).

See Also

data_helpers

Examples

# Create a dataset from scratch
obs_rows(time = 12, DV = 0.12, cmt = 2)

# Pipe-friendly addition of observation record to a pre-existing dataset
library(magrittr)
obs_rows(time = 12, DV = 0.12, cmt = 2) %>%
  obs_rows(time = c(24, 36, 48), DV = c(0.34, 0.56, 0.78), mdv = c(0, 1, 0), cmt = 2)

# Inform times using the `.datehour` argument:
obs_rows(.datehour = "2020-01-01 11:11", DV = 0.12, cmt = 1) %>%
  obs_rows(.datehour = "2020-01-02 22:22", DV = 0.34, cmt = 1) %>%
  obs_rows(time = 48, DV = 0.56, cmt = 1)

# Start from a 'mrgsolve' model
library(mrgsolve)
house() %>%
  obs_rows(time = 12, DV = 0.12, cmt = 2) %>%
  obs_rows(time = c(24, 36, 48), DV = c(0.34, 0.56, 0.78), mdv = c(0, 1, 0), cmt = 2) %>%
  mrgsim()

# Default observation compartments
# Set default observation compartments in the code with `[OBS]`
model <- mcode("model", "
$CMT @annotated
DEPOT : Depot
CENTR : Central [OBS]
", compile = FALSE)
obs_cmt(model)

# Thus, no need to manually specify `cmt = 2` anymore.
model %>%
  obs_rows(time = 12, DV = 0.12) %>%
  obs_rows(time = c(24, 36, 48), DV = c(0.34, 0.56, 0.78), mdv = c(0, 1, 0)) %>%
  get_data()

# Automatic lines duplication if parent + metabolite defined in the model
model <- mcode("model", "
$CMT @annotated
DEPOT : Depot
CENTR : Central [OBS]
PERIPH : Periph
METABO : Metabo [OBS]
", compile = FALSE)
obs_cmt(model)

model %>%
  obs_rows(time = 12, DV = 0.12, DVmet = 120) %>%
  obs_rows(
    time = c(24, 36, 48), DV = c(0.34, 0.56, 0.78),
    mdv = c(0, 1, 0), DVmet = c(340, 560, 780)
  ) %>%
  get_data()

Parse value to "POSIXct"

Description

A wrapper around functions of lubridate, mainly in order to transform characters into a date-time ("POSIXct") format.

Usage

parse_datehour(
  x,
  orders = getOption("mapbayr.datehour", default = c("Ymd HMS", "Ymd HM", "dmY HMS",
    "dmY HM"))
)

Arguments

Value

a POSIXct

Examples

# POSITct are returned as is.
parse_datehour(x = as.POSIXct("2022-02-02 22:22:22", tz = "UTC"))

# Numerics are passed to `lubridate::as_datetime()`.
parse_datehour(1643840542)

# Characters are passed to `lubridate::parse_date_time()`.
# The format used will be the one defined in `orders`
parse_datehour(x = "2022-02-02 22:22:22", orders = "Ymd HMS")
parse_datehour(x = "02-02-2022 22:22", orders = "dmY HM")

# By default, the following formats will be subsequently tried:
# "Ymd HMS", "Ymd HM", "dmY HMS", "dmY HM"

# Alternatively, set a format through `options(mapbayr.datehour)`.
# Convenient for the use `.datehour` in  `adm_rows()` and `obs_rows()`.

# Following format will return NA:
adm_rows(.datehour = "22:22 02-02-2022", amt = 100, cmt = 1)

options(mapbayr.datehour = "HM dmY")
adm_rows(.datehour = "22:22 02-02-2022", amt = 100, cmt = 1)
options(mapbayr.datehour = NULL)

Plot predictions from mapbayests object

Description

Plot predictions from mapbayests object

Usage

## S3 method for class 'mapbayests'
plot(x, ..., PREDICTION = c("IPRED", "PRED"))

Arguments

Details

Use this function to plot the results of the estimations, in the form of concentration vs time profiles for every patient of the data set. For additional modifications, you can:

• see augment.mapbayests to modify the simulation output.

• add extra ⁠+function(...)⁠ in order to modify the plot as a regular ggplot2 object.

Value

a ggplot object.

Examples

est <- mapbayest(exmodel(ID = 1))
plot(est, end = 48) +
  ggplot2::labs(title = "Awesome prediction")

Preprocess model and data for ofv computation

Description

Functions to generate arguments passed to compute_ofv. Arguments that are fixed between individuals are created once (preprocess.ofv.fix), while others are specific of each individual (preprocess.ofv.id).

Usage

preprocess.ofv.fix(x, data, select_eta = seq_along(eta(x)), lambda = 1)

preprocess.ofv.id(x, iddata)

Arguments

Value

A list of arguments used to compute the objective function value.

The following arguments are fixed between individuals:

• qmod: model object, modified to simulate without random effects and with controlled outputs

• sigma: a single matrix object

• log_transformation: a logical, whether predictions need to be log-transformed for ofv computation

• omega_inv: a single matrix object

• all_cmt: a vector of compartment numbers where observations can be expected

The following arguments differs between individuals:

• idvaliddata: a matrix, individual data set (with administrations and covariates), validated with valid_data_set

• idDV: a vector of (possibly log-transformed) observations

• idcmt: a vector of compartments where observations belong to

• idblq,idlloq: optional, a logical and numerical vector indicating if the observation is below the lower limit of quantification, and the LLOQ value, respectively

Examples

mod <- exmodel(add_exdata = FALSE, compile = FALSE)
dat <- exdata(ID = c(1,4))

preprocess.ofv.fix(x = mod, data = dat)
preprocess.ofv.id(x = mod, iddata = dat[dat$ID == 1,])
preprocess.ofv.id(x = mod, iddata = dat[dat$ID == 4,])

Pre-process: arguments for optimization function

Description

Pre-process: arguments for optimization function

Usage

preprocess.optim(
  x,
  method = c("L-BFGS-B", "newuoa"),
  select_eta = NULL,
  control = list(),
  force_initial_eta = NULL,
  quantile_bound = 0.001
)

Arguments

Value

a list of named arguments passed to optimizer (i.e. arg.optim)

Print a mapbayests object

Description

Print a mapbayests object

Usage

## S3 method for class 'mapbayests'
print(x, ...)

Arguments

Value

print the results of the estimation to the console, and returns it invisibly.

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

dplyr

filter

Use posterior estimation

Description

Use posterior estimation

Usage

use_posterior(
  x,
  update_omega = FALSE,
  update_cov = TRUE,
  update_eta = TRUE,
  .zero_re = NULL,
  simplify = TRUE
)

Arguments

Details

This function takes the results of an estimation (i.e. a mapbayests object) and return a modified mrgmod in order to perform a posteriori simulations. Modifications are:

• If update_eta is TRUE, the values of ETA are updated to the estimated values (instead of 0) in $PARAM.

• If update_cov is TRUE, the covariates values are updated to the values of the individual (instead of default model values) in $PARAM.

• If update_omega is TRUE, the values of OMEGA are updated with the variance-covariance matrix of estimation (i.e. an approximation of the a posteriori distribution) instead of the inter-individual variability (i.e. the a priori distribution). Use this command in order to derive a confidence interval of concentrations that reflects the uncertainty about parameter estimation when a large number of profiles are simulated. Note that if inter-individual variability was initially defined in multiple $OMEGA blocks in the model, they will be collapsed to a single full matrix (this is irreversible).

• Depending on the values of .zero_re, the elements of $OMEGA or $SIGMA can be set to zero, whether you want to simulate one profile, or several in order to derive confidence/prediction intervals. It does not handle time-varying covariates: only the first value will be used as the individual value.

Value

a mrgmod, or a list of mrgmod if there is more than 1 ID

Examples

library(magrittr)
est <- mapbayest(exmodel())
est %>%
  use_posterior() %>%
  mrgsolve::ev(amt = 50000) %>%
  mrgsolve::mrgsim()

Compare results to NONMEM .phi

Description

Compare results to NONMEM .phi

Usage

read_nmphi(x)

merge_phi(mapbayr_phi, nonmem_phi)

plot_phi(merged_phi, only_ETA = TRUE)

summarise_phi(
  merged_phi,
  group,
  only_ETA = TRUE,
  levels = c(Excellent = 0, Acceptable = 0.001, Discordant = 0.1)
)

bar_phi(summarised_phi, xaxis = NULL, facet = NULL)

Arguments

Details

These functions were made to easily compare the results of mapbayr to NONMEM. For instance, it could be useful in the case of the transposition of a pre-existing NONMEM model into mapbayr. For this, you need to code your model in both mapbayr and NONMEM, and perform the MAP-Bayesian estimation on the same dataset. Ideally, the latter contains a substantial number of patients. NONMEM returns the estimations results into a .phi file.

Use read_nmphi() to parse the NONMEM .phi file into a convenient tibble data.frame with the columns:

• SUBJECT_NO, ID: Subject identification.

• ETA1, ETA2, ..., ETAn: Point estimates of eta.

• ETC1_1, ETC2_1, ETC2_2, ..., ETCn_n: Variance-covariance matrix of estimation.

• OBJ: objective function value

Use get_phi() to access to the estimations of mapbayr with the same "phi" format.

Use merge_phi() to combine mapbayr and NONMEM "phi files" into a single long-form data.frame with the columns:

• SUBJECT_NO, ID: Subject identification.

• variable name and its type: ETA (point estimate), VARIANCE (on-diagonal element of the matrix), COVARIANCE (off-diagonal), and OBJ.

• mapbayr and nonmem: corresponding values

• adiff: absolute difference between mapbayr and nonmem values.

Use plot_phi() to graphically represent adiff vs variable. Alternatively, the table returned by merge_phi() is easy to play with in order to derive performance statistics or the graphical plot of your choice.

Use summarise_phi() to classify the estimation as "Excellent", "Acceptable" or "Discordant", over the whole dataset or by group.

Use bar_phi() to graphically represent the proportion of the aforementioned classification as bar plot.

Value

• read_nmphi: a tibble data.frame with a format close to the original .phi file

• merge_phi: a long-form tibble data.frame with results of mapbayr and NONMEM

• summarise_phi: a summarized tibble data.frame classifying the performance of mapbayr and NONMEM

• plot_phi, bar_phi: a ggplot2 object

Examples

library(mapbayr)
nmphi <- read_nmphi(system.file("nm001", "run001.phi", package = "mapbayr"))
mapbayrphi <- get_phi(est001)

merged <- merge_phi(mapbayrphi, nmphi)
plot_phi(merged)

summarised <- summarise_phi(merged)
bar_phi(summarised)


# Analyse the results of multiple runs simultaneously

#Example dataset that represents 3 runs
merge3 <- dplyr::bind_rows(merged, merged, merged, .id = "RUN")
merge3$adiff <- 10 ^ runif(nrow(merge3), -8, 0)

summarised3 <- summarise_phi(merge3, group = RUN)
bar_phi(summarised3, xaxis = "RUN")

Read compartment options in a model

Description

Read compartment options in a model

Usage

adm_cmt(x)

obs_cmt(x)

Arguments

Details

In a mrgsolve model, it is possible to specify options in ⁠$CMT⁠. If ⁠[ADM]⁠ or ⁠[OBS]⁠ are set, mapbayr will interpret these as defaults administration and observation compartments, respectively.

Value

a vector of compartment identified as default "administration" or "observation" compartments.

Examples

#Administration:  Both 1st and 0- order
model <- exmodel(6, compile = FALSE)
mrgsolve::see(model)
adm_cmt(model)

#Observation: Both parent drug and metabolite
model <- exmodel(401, compile = FALSE)
mrgsolve::see(model)
obs_cmt(model)