Package 'admixr2'

Description

Creates a control object for nlmixr2(est = "adfo"). The FO estimator linearises model predictions at $\eta = 0$: it is faster than the MC estimator but less accurate for models with large IIV or strongly non-linear individual predictions.

Usage

adfoControl(
  studies = list(),
  grad = c("none", "analytical", "fd", "cfd"),
  algorithm = NULL,
  maxeval = 500L,
  ftol_rel = .Machine$double.eps^(1/2),
  print = 10L,
  seed = 12345L,
  cores = rxode2::rxCores(),
  nDisplayProgress = .Machine$integer.max,
  grad_h = 1e-04,
  grad_bounds = 5,
  cov_h = 0.001,
  cov_h_outer = .Machine$double.eps^(1/5),
  covMethod = c("r", "none"),
  n_restarts = 1L,
  restart_sd = 0.5,
  workers = 1L,
  rxControl = NULL,
  calcTables = FALSE,
  compress = TRUE,
  ci = 0.95,
  sigdig = 4,
  sigdigTable = NULL,
  addProp = c("combined2", "combined1"),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  returnAdmr = FALSE,
  ...
)

Arguments

Value

An adfoControl object (a named list).

Installing memuse

rxode2::rxSolve() estimates free RAM on every call. When the memuse package is not installed its fallback ends up shelling out to vm_stat, a macOS-only command, so on Windows and Linux every solve spawns a process that can only fail. Because the FO estimator issues many small solves, this overhead is measurable (roughly 17% of an FO gradient). Installing memuse makes the fallback unreachable:

install.packages("memuse")

See Also

admControl(), adirmcControl()

Examples

# Inspect defaults
ctl <- adfoControl()
ctl$grad
ctl$maxeval

# Analytical gradient, more evaluations
ctl2 <- adfoControl(grad = "analytical", maxeval = 1000L)


library(rxode2)
library(nlmixr2)

data("examplomycin")
obs    <- examplomycin[examplomycin$EVID == 0, ]
obs    <- obs[order(obs$ID, obs$TIME), ]
times  <- sort(unique(obs$TIME))
ids    <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adfo",
  control = adfoControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100))),
    maxeval = 100L
  )
)
print(fit)

Description

Creates a control object for nlmixr2(est = "adgh"). The GH estimator integrates model predictions against the random-effects prior $\eta \sim N(0, \Omega)$ using a deterministic tensor-product Gauss-Hermite quadrature grid. It is unbiased at any IIV magnitude (unlike FO), noise-free (unlike MC), and much faster than MC for models with up to ~4 etas.

Usage

adghControl(
  studies = list(),
  n_nodes = 5L,
  grad = c("analytical", "fd", "cfd", "none"),
  algorithm = "NLOPT_LN_BOBYQA",
  maxeval = 500L,
  ftol_rel = .Machine$double.eps^(1/2),
  print = 10L,
  seed = 12345L,
  cores = rxode2::rxCores(),
  nDisplayProgress = .Machine$integer.max,
  grad_h = 1e-04,
  grad_bounds = 5,
  cov_h = 0.001,
  cov_h_outer = .Machine$double.eps^(1/4),
  covMethod = c("r", "none"),
  n_restarts = 1L,
  restart_sd = 0.5,
  workers = 1L,
  rxControl = NULL,
  calcTables = FALSE,
  compress = TRUE,
  ci = 0.95,
  sigdig = 4,
  sigdigTable = NULL,
  addProp = c("combined2", "combined1"),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  returnAdmr = FALSE,
  ...
)

Arguments

Value

An adghControl object (a named list).

See Also

admControl(), adfoControl(), adirmcControl()

Examples

ctl <- adghControl()
ctl$n_nodes
ctl$grad

# More nodes for large IIV, analytical gradient
ctl2 <- adghControl(n_nodes = 7L, grad = "analytical", maxeval = 300L)


library(rxode2)
library(nlmixr2)

data("examplomycin")
obs    <- examplomycin[examplomycin$EVID == 0, ]
obs    <- obs[order(obs$ID, obs$TIME), ]
times  <- sort(unique(obs$TIME))
ids    <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adgh",
  control = adghControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100)))
  )
)

Description

Constructs a control object for est = "adirmc", the Iterative Reweighting Monte Carlo estimator.

Usage

adirmcControl(
  studies = list(),
  n_sim = 2500L,
  outer_iter = 50L,
  sampling = c("sobol", "halton", "torus", "lhs", "rnorm"),
  algorithm = NULL,
  maxeval = 5000L,
  ftol_rel = .Machine$double.eps,
  print = 1L,
  omega_expansion = 1,
  seed = 12345L,
  cores = rxode2::rxCores(),
  nDisplayProgress = .Machine$integer.max,
  grad = c("analytical", "none", "fd"),
  kappa_method = c("exact", "linearized", "linearized_gh"),
  kappa_n_nodes = 5L,
  grad_h = 1e-04,
  cov_h = 0.001,
  cov_h_outer = .Machine$double.eps^(1/5),
  phases = c(2, 1, 0.5, 0.01),
  convcrit = 1e-05,
  max_worse = 5L,
  covMethod = c("r", "none"),
  cov_n_sim = 10000L,
  n_restarts = 1L,
  restart_sd = 0.2,
  workers = 1L,
  rxControl = NULL,
  calcTables = FALSE,
  compress = TRUE,
  ci = 0.95,
  sigdig = 4,
  sigdigTable = NULL,
  addProp = c("combined2", "combined1"),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  returnAdmr = FALSE,
  ...
)

Arguments

# independent blocks: per-row variances; optional per-endpoint `n` column
# and per-endpoint `ev` (a list of event tables keyed by endpoint)
list(n = 60L, ev = ev,
     data = data.frame(DVID = c("cp", "cp", "cCSF"), TIME = c(1, 2, 2),
                       E = c(9.1, 7.4, 2.2), V = c(1.2, 0.9, 0.1)))

# joint (same subjects): ONE stacked covariance whose rows/cols align with
# the rows of `data` -- no `cross` blocks to assemble by hand
list(n = 60L, ev = ev, data = data.frame(DVID = ..., TIME = ..., E = ...),
     V = V_joint)

Details

Multi-compartment fits (a study observations list with several observed outputs) are not supported by adirmc; use est = "admc", "adfo", or "adgh" for those. Single-output studies are fit as usual.

Value

An object of class adirmcControl.

Examples

# Inspect defaults
ctl <- adirmcControl()
ctl$phases
ctl$omega_expansion

# Tighter phases, more restarts
ctl2 <- adirmcControl(
  n_sim           = 1000L,
  omega_expansion = 1.5,
  phases          = c(2, 1, 0.5, 0.01),
  n_restarts      = 3L
)


library(rxode2)
library(nlmixr2)

data("examplomycin")
obs   <- examplomycin[examplomycin$EVID == 0, ]
obs   <- obs[order(obs$ID, obs$TIME), ]
times <- sort(unique(obs$TIME))
ids   <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- diag(diag(cov.wt(dv_mat, method = "ML")$cov))

pk_model <- function() {
  ini({
    tcl <- log(5);  tv1 <- log(12); tv2 <- log(25)
    tq  <- log(12); tka <- log(1.2)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v1 ~ 0.09; eta.v2 ~ 0.09
    eta.q  ~ 0.09; eta.ka ~ 0.09
  })
  model({
    cl <- exp(tcl + eta.cl); v1 <- exp(tv1 + eta.v1)
    v2 <- exp(tv2 + eta.v2); q  <- exp(tq  + eta.q)
    ka <- exp(tka + eta.ka)
    d/dt(depot)      <- -ka * depot
    d/dt(central)    <- ka * depot - (cl/v1 + q/v1) * central + (q/v2) * peripheral
    d/dt(peripheral) <- (q/v1) * central - (q/v2) * peripheral
    cp <- central / v1
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adirmc",
  control = adirmcControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100))),
    n_sim   = 500L
  )
)
print(fit)

Description

Constructs a control object for est = "admc", the Monte Carlo aggregate data modelling estimator.

Usage

admControl(
  studies = list(),
  n_sim = 5000L,
  sampling = c("sobol", "halton", "torus", "lhs", "rnorm"),
  algorithm = NULL,
  maxeval = 500L,
  ftol_rel = .Machine$double.eps^2,
  print = 10L,
  seed = 12345L,
  cores = rxode2::rxCores(),
  nDisplayProgress = .Machine$integer.max,
  grad = c("sens", "fd", "cfd", "none"),
  grad_h = 1e-04,
  cov_h = 0.001,
  cov_h_outer = .Machine$double.eps^(1/5),
  grad_bounds = 5,
  covMethod = c("r", "none"),
  cov_n_sim = 10000L,
  n_restarts = 1L,
  restart_sd = 0.5,
  workers = 1L,
  rxControl = NULL,
  calcTables = FALSE,
  compress = TRUE,
  ci = 0.95,
  sigdig = 4,
  sigdigTable = NULL,
  addProp = c("combined2", "combined1"),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  returnAdmr = FALSE,
  ...
)

Arguments

# independent blocks: per-row variances; optional per-endpoint `n` column
# and per-endpoint `ev` (a list of event tables keyed by endpoint)
list(n = 60L, ev = ev,
     data = data.frame(DVID = c("cp", "cp", "cCSF"), TIME = c(1, 2, 2),
                       E = c(9.1, 7.4, 2.2), V = c(1.2, 0.9, 0.1)))

# joint (same subjects): ONE stacked covariance whose rows/cols align with
# the rows of `data` -- no `cross` blocks to assemble by hand
list(n = 60L, ev = ev, data = data.frame(DVID = ..., TIME = ..., E = ...),
     V = V_joint)

Value

An object of class admControl.

Examples

# Minimal control object -- inspect defaults
ctl <- admControl()
ctl$n_sim
ctl$algorithm

# Override key settings without fitting
ctl2 <- admControl(
  n_sim    = 2000L,
  maxeval  = 300L,
  grad     = "fd",
  seed     = 42L
)


library(rxode2)
library(nlmixr2)

data("examplomycin")
obs   <- examplomycin[examplomycin$EVID == 0, ]
obs   <- obs[order(obs$ID, obs$TIME), ]
times <- sort(unique(obs$TIME))
ids   <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5);  tv1 <- log(12); tv2 <- log(25)
    tq  <- log(12); tka <- log(1.2)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v1 ~ 0.09; eta.v2 ~ 0.09
    eta.q  ~ 0.09; eta.ka ~ 0.09
  })
  model({
    cl <- exp(tcl + eta.cl); v1 <- exp(tv1 + eta.v1)
    v2 <- exp(tv2 + eta.v2); q  <- exp(tq  + eta.q)
    ka <- exp(tka + eta.ka)
    d/dt(depot)      <- -ka * depot
    d/dt(central)    <- ka * depot - (cl/v1 + q/v1) * central + (q/v2) * peripheral
    d/dt(peripheral) <- (q/v1) * central - (q/v2) * peripheral
    cp <- central / v1
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "admc",
  control = admControl(
    studies  = list(study1 = list(E = E, V = V, n = length(ids),
                                  times = times, ev = et(amt = 100))),
    n_sim    = 1000L,
    maxeval  = 200L
  )
)
print(fit)

Description

Returns a minimal NONMEM-style data frame that satisfies nlmixr2's data argument requirement. The single observation row carries a non-NA placeholder DV (1); the dose row keeps DV = NA. The placeholder is purely for dispatch and output construction and never enters the reported objective – each estimator overwrites fit$env$objective (and OBJF/logLik/AIC/BIC) with its own aggregate -2LL. A non-NA observation is required because nlmixr2's post-fit output construction (nlmixr2CreateOutputFromUi) solves the model over this frame, and rxode2's event-table translation rejects a dataset with no non-NA observation rows ("no rows in event table or input data"), the same reason the multi-endpoint frame below uses a placeholder DV.

Usage

admData(outputs = NULL)

Arguments

Value

A data frame with columns ID, TIME, DV, AMT, EVID, CMT (single-endpoint), plus a DVID endpoint column when outputs is given.

Examples

admData()
admData(c("cp", "cCSF"))

Description

Stops any worker processes (mirai daemons) started by a parallel-restart fit (admControl(workers = N)). Workers are stopped automatically after the restart phase completes, so this function is only needed if a fit was interrupted before cleanup could run.

Usage

admStopWorkers()

Value

NULL, invisibly.

Examples

# Safe to call at any time; no-op if no workers are running
admStopWorkers()

Description

Generates population mean vectors (E) and covariance matrices (V) for each study by integrating over the IIV distribution – either by Monte Carlo (the default) or by a deterministic First-Order expansion (method = "fo", see datagenControl()). Each study may specify its own PK/PD model (as would be the case when digitising data from several published studies, each fit with a different structural model). True parameter values are taken from the ini() block of each study's model. Each element of the returned list is ready to supply directly to admControl(studies = ...).

Usage

datagen(studies, model = NULL, control = datagenControl())

Arguments

Details

With control = datagenControl(method = "mc") (the default) population moments are computed via the same Monte Carlo engine as est = "admc":

$E_t = \bar{f}_s(\hat\theta_s, \eta_i, t)$

$V_{ts} = \widehat{\mathrm{Cov}}_\eta[f_{s,t}, f_{s,s'}] + \Sigma_s$

where $f_s$ and $\hat\theta_s$ are the model and initial estimates from the ini() block of study $s$, the sample covariance uses the ML denominator n_sim, and $\Sigma_s$ is diagonal with entries determined by that study model's residual error type (additive, proportional, or log-normal).

With method = "fo" the moments are instead the deterministic First-Order expansion used by est = "adfo":

$E = f_s(\hat\theta_s, 0)$

$V = J \Omega_s J^\top + \Sigma_s, \quad J_{tj} = \partial f_{s,t}/\partial \eta_j |_{\eta = 0}$

with the Jacobian $J$ obtained from the sensitivity model (or finite differences if that is unavailable). This is the natural choice for design evaluation and optimal design: the moments are fast and reproducible, and because the data-generating and data-analytic models coincide, the FO Hessian of the log-likelihood (the expected information matrix) is evaluated at the true maximum rather than at a point that is not an MLE of the generated data. Note est = "adfo" always adds $\Sigma$ to its predicted covariance, so for a consistent FIM keep the residual error in the generating model; omit it only when residual-free (IIV-only) moments are genuinely what you want.

With method = "gh" the moments are computed by deterministic Gauss-Hermite quadrature over the random-effects prior $\eta \sim N(0, \Omega)$:

$E = \sum_q w_q f(\hat\theta, \eta_q), \quad V = \sum_q w_q (f_q - E)(f_q - E)^\top + \Sigma$

where $(\eta_q, w_q)$ are the Cholesky-scaled tensor-product GH nodes and weights. Unlike FO this is unbiased at any IIV magnitude; unlike MC the result is noise-free and exactly reproducible. Matching the moments of est = "adgh" makes method = "gh" the natural choice for optimal design with that estimator.

Models are compiled and cached on first use (keyed by model expression digest), so repeated calls or multiple studies sharing the same model incur only a single compilation.

Value

A named list with one element per study. Each element contains:

E

Population mean vector at times.

V

Population covariance matrix (length(times) x length(times); ML denominator n_sim for method = "mc", the analytical FO covariance for method = "fo", or the GH weighted covariance for method = "gh"). The diagonal carries the model's residual-error variance; to generate residual-free (IIV-only) moments, omit the error term from the model.

n

Sample size (NA_integer_ if not supplied).

times

Observation times.

ev

Dosing event table.

samples

Raw ⁠n_sim x length(times)⁠ prediction matrix (only when control$return_samples = TRUE).

See Also

datagenControl(), admControl()

Examples

library(rxode2)

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

study_data <- datagen(
  studies = list(
    study1 = list(times = c(1, 2, 4, 8, 12, 24),
                  ev = rxode2::et(amt = 100), n = 200L)
  ),
  model   = pk_model,
  control = datagenControl(n_sim = 2000L)
)

# E and V plug directly into admControl(studies = ...)
round(study_data$study1$E, 2)

Control parameters for datagen()

Description

Control parameters for datagen()

Usage

datagenControl(
  method = c("mc", "fo", "gh"),
  n_sim = 5000L,
  n_nodes = 5L,
  sampling = c("sobol", "halton", "torus", "lhs", "rnorm"),
  seed = 12345L,
  cores = 1L,
  return_samples = FALSE
)

Arguments

Value

A list of class "datagenControl".

See Also

datagen()

Examples

ctrl <- datagenControl(n_sim = 2000L)
ctrl$sampling  # "sobol"

# Deterministic FO moments for design evaluation:
datagenControl(method = "fo")$method  # "fo"

# GH quadrature moments (unbiased, noise-free):
datagenControl(method = "gh", n_nodes = 5L)$n_nodes

Description

A simulated pharmacokinetic dataset for the fictional drug examplomycin, intended as a worked example for aggregate data modelling with admixr2. The dataset contains 500 subjects, each with 9 observation time points, generated from a two-compartment model with first-order absorption.

Usage

examplomycin

Format

A data frame with 5000 rows and 6 columns:

• ID: Subject identifier (integer, 1–500).

• TIME: Time after dose (hours).

• DV: Observed plasma concentration (mg/L).

• AMT: Dose amount (mg); 100 for dosing records, 0 otherwise.

• EVID: Event type (101 = dose, 0 = observation).

• CMT: Compartment (1 = depot, 2 = central).

Details

True population parameters:

Single oral dose of 100 mg; sampling at 0.1, 0.25, 0.5, 1, 2, 3, 5, 8, and 12 hours post-dose.

Source

Generated from a two-compartment PK model using rxode2::rxSolve(). See vignette("admixr2") for a full modelling example.

Examples

data("examplomycin")
head(examplomycin)

# Compute aggregate statistics
obs <- examplomycin[examplomycin$EVID == 0, ]
obs <- obs[order(obs$ID, obs$TIME), ]
times <- sort(unique(obs$TIME))
E <- sapply(times, function(t) mean(obs$DV[obs$TIME == t]))
round(E, 3)

Description

Called automatically by nlmixr2(model, admData(), est = "adfo", control = adfoControl(...)). Not typically called directly.

Usage

## S3 method for class 'adfo'
nlmixr2Est(env, ...)

Arguments

Value

An admFit nlmixr2 fit object.

Description

Called automatically by nlmixr2(model, admData(), est = "adgh", control = adghControl(...)). Not typically called directly.

Usage

## S3 method for class 'adgh'
nlmixr2Est(env, ...)

Arguments

Value

An admFit nlmixr2 fit object.

Description

Called automatically by nlmixr2(model, admData(), est = "adirmc", control = adirmcControl(...)). Not typically called directly.

Usage

## S3 method for class 'adirmc'
nlmixr2Est(env, ...)

Arguments

Value

An admFit nlmixr2 fit object.

Description

Called automatically by nlmixr2(model, admData(), est = "admc", control = admControl(...)). Not typically called directly.

Usage

## S3 method for class 'admc'
nlmixr2Est(env, ...)

Arguments

Value

An admFit nlmixr2 fit object.

Description

Generates up to four diagnostic panels:

Usage

## S3 method for class 'admFit'
plot(x, which = c("mean", "cov", "nll", "par"), n_sim = NULL, seed = 1L, ...)

Arguments

Details

1. "mean" – Observed vs predicted mean per study (2x2 grid). Upper row: observed and predicted mean lines with +/-1 SD ribbon on a shared y scale (black throughout). Lower row: raw residual lollipop with +/-2 SE band and standardised residual z-scores with +/-1.96 reference lines.

2. "cov" – Observed vs predicted (co)variance heatmaps per study (2x2 grid). Upper row shares a common colour scale (blue-white-red). Lower row uses distinct diverging scales: residual (red-white-green) and standardised residual (gold-white-purple). Significance stars overlaid on the standardised residual panel.

3. "nll" – NLL trace per restart over optimizer evaluations. Restarts coloured with the Okabe-Ito palette.

4. "par" – Parameter trace per restart on the natural scale (struct thetas back-transformed, sigma as SD, omega diagonal as variance labelled V(eta.x)). Facets ordered as in the model ini() block. Restarts coloured with the Okabe-Ito palette.

Value

A named list of ggplot2 objects, invisibly. Prints each selected top-level panel. For the "mean" and "cov" panels the returned list also contains each sub-panel individually so a single panel (or a few) can be extracted in code without reprinting the whole grid. Elements can be pulled out by name – plot(fit, which = "mean")$mean_study1_pred or plot(fit, which = "cov")$cov_study1_std_resid – or by position, with the combined 2x2 grid stored first per study (plot(fit, which = "mean")[[1]] is the full grid, ⁠[1]⁠ the length-1 named sub-list). The sub-panel keys are ⁠<type>_<study>_obs⁠, ⁠_pred⁠, ⁠_resid⁠, and ⁠_std_resid⁠; the combined grid stays under ⁠<type>_<study>⁠. The extra sub-panel keys are not printed on their own.

Aggregate data slot

Every admixr2 fit also carries the observed and predicted aggregate data in fit$env$aggData, a named list with one entry per study. Each entry holds the observation times, the study n, and two moment sets – obs (from the data) and pred (predicted at the fitted parameters) – each a list with the mean vector E and the (co)variance matrix V:

  fit$env$aggData$study1$obs$E    # observed mean vector
  fit$env$aggData$study1$obs$V    # observed covariance matrix
  fit$env$aggData$study1$pred$E   # predicted mean vector
  fit$env$aggData$study1$pred$V   # predicted covariance matrix

The predicted moments are computed by one MC simulation at the fitted parameters using the fit's own n_sim and a fixed seed, so they match the default plot(fit) mean/cov panels. The slot is absent only when the fit cannot be simulated (no simulation model available).

nlmixr2 traceplot()

admixr2 fits also plug into the nlmixr2 traceplot() generic. During fitting the parameter iteration history of the best restart is stored on the fit in the standard parHistData slot (natural scale), so traceplot(fit) produces the familiar per-parameter, free-y facetted trace used elsewhere in the nlmixr2 ecosystem. There is no burn-in marker (admixr2 records optimizer evaluations, not SAEM iterations), and only the best restart is shown – the per-restart overlay and the NLL trace remain available via plot(fit, which = c("par", "nll")). The trace stores only improving evaluations (steps that lowered the best NLL), so the iter axis indexes those improvement steps rather than raw optimizer iterations.

Examples

library(rxode2)
library(nlmixr2)

data("examplomycin")
obs    <- examplomycin[examplomycin$EVID == 0, ]
obs    <- obs[order(obs$ID, obs$TIME), ]
times  <- sort(unique(obs$TIME))
ids    <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adfo",
  control = adfoControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100))),
    maxeval = 100L
  )
)
plot(fit)

Description

Delegates to print.nlmixr2FitCore for the standard nlmixr2 coloured output. admFit class is kept on the object during the call so that head.admFit intercepts any head(fit) calls that arise in the paged- output path (R Markdown / notebooks), preventing the ⁠[.data.frame(.subset2(env, integer))⁠ crash that occurs when an environment-backed fit is subscripted like a plain list.

Usage

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

Arguments

Value

x, invisibly.

Examples

library(rxode2)
library(nlmixr2)

data("examplomycin")
obs    <- examplomycin[examplomycin$EVID == 0, ]
obs    <- obs[order(obs$ID, obs$TIME), ]
times  <- sort(unique(obs$TIME))
ids    <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adfo",
  control = adfoControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100))),
    maxeval = 100L
  )
)
print(fit)