poso_estim_mcmc {posologyr}

Estimate the posterior distribution of individual parameters by MCMC

Description

Estimates the posterior distribution of individual parameters by Markov Chain Monte Carlo (using a Metropolis-Hastings algorithm)

Usage

poso_estim_mcmc(
  dat = NULL,
  prior_model = NULL,
  return_model = TRUE,
  burn_in = 50,
  n_iter = 1000,
  n_chains = 4,
  nocb = FALSE,
  control = list(n_kernel = c(2, 2, 2), stepsize_rw = 0.4, proba_mcmc = 0.3, nb_max = 3)
)

Arguments

dat

Dataframe. An individual subject dataset following the structure of NONMEM/rxode2 event records.

prior_model

A posologyr prior population pharmacokinetics model, a list of six objects.

return_model

A boolean. Returns a rxode2 model using the estimated ETAs if set to TRUE.

burn_in

Number of burn-in iterations for the Metropolis-Hastings algorithm.

n_iter

Total number of iterations (following the burn-in iterations) for each Markov chain of the Metropolis-Hastings algorithm.

n_chains

Number of Markov chains

nocb

A boolean. for time-varying covariates: the next observation carried backward (nocb) interpolation style, similar to NONMEM. If FALSE, the last observation carried forward (locf) style will be used. Defaults to FALSE.

control

A list of parameters controlling the Metropolis-Hastings algorithm.

Value

If return_model is set to FALSE, a list of one element: a dataframe ⁠$eta⁠ of ETAs from the posterior distribution, estimated by Markov Chain Monte Carlo. If return_model is set to TRUE, a list of the dataframe of the posterior distribution of ETA, and a rxode2 model using the estimated distributions of ETAs.

Author(s)

Emmanuelle Comets, Audrey Lavenu, Marc Lavielle, Cyril Leven

References

Comets E, Lavenu A, Lavielle M. Parameter estimation in nonlinear mixed effect models using saemix, an R implementation of the SAEM algorithm. Journal of Statistical Software 80, 3 (2017), 1-41.

Examples

# model
mod_run001 <- list(
ppk_model = rxode2::rxode({
  centr(0) = 0;
  depot(0) = 0;

  TVCl = THETA_Cl;
  TVVc = THETA_Vc;
  TVKa = THETA_Ka;

  Cl = TVCl*exp(ETA_Cl);
  Vc = TVVc*exp(ETA_Vc);
  Ka = TVKa*exp(ETA_Ka);

  K20 = Cl/Vc;
  Cc = centr/Vc;

  d/dt(depot) = -Ka*depot;
  d/dt(centr) = Ka*depot - K20*centr;
  d/dt(AUC) = Cc;
}),
error_model = function(f,sigma) {
  dv <- cbind(f,1)
  g <- diag(dv%*%sigma%*%t(dv))
  return(sqrt(g))
},
theta = c(THETA_Cl=4.0, THETA_Vc=70.0, THETA_Ka=1.0),
omega = lotri::lotri({ETA_Cl + ETA_Vc + ETA_Ka ~
    c(0.2,
      0, 0.2,
      0, 0, 0.2)}),
sigma = lotri::lotri({prop + add ~ c(0.05,0.0,0.00)}))
# df_patient01: event table for Patient01, following a 30 minutes intravenous
# infusion
df_patient01 <- data.frame(ID=1,
                        TIME=c(0.0,1.0,14.0),
                        DV=c(NA,25.0,5.5),
                        AMT=c(2000,0,0),
                        EVID=c(1,0,0),
                        DUR=c(0.5,NA,NA))
# estimate the posterior distribution of population parameters
poso_estim_mcmc(dat=df_patient01,prior_model=mod_run001,
n_iter=50,n_chains=2)