---
title: "Vancomycin (Jung 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vancomycin (Jung 2024)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Jung_2024_vancomycin"))$reference`
* Description: `r rxode2::rxode(readModelDb("Jung_2024_vancomycin"))$description`
* Article: <https://doi.org/10.1097/CCE.0000000000001159>
* Supplement (eMethods, eTables 1-4, eFigures 1-10): <http://links.lww.com/CCX/B406>

Jung 2024 is a prospective, two-site population PK study of IV vancomycin in
critically ill children. The paper has three quantitative products, and this
vignette exercises all three:

1. **A two-compartment population PK model** (Table 2) - the model packaged here
   as `Jung_2024_vancomycin`.
2. **A nonlinear regression predicting AUC24 from a single trough
   concentration** (eTable 3, tabulated in Table 3). This is a regression, not
   an ODE model, so it is not packaged as a separate model file; it is
   reproduced below as a plain function and used as an independent cross-check.
3. **A target-attainment simulation** of four dosing regimens against
   AUC24/MIC >= 400 (Table 4, with the underlying AUC24 distributions in
   eTable 4) - reproduced below from the packaged model.

# Population

The model was fit to 301 critically ill children (302 enrolled between
2018-04-15 and 2020-02-06; one excluded for unavailable PK data) aged 90 days to
under 18 years, admitted to the PICU at the University of Maryland Children's
Hospital or Texas Children's Hospital and prescribed IV vancomycin for any
indication (Jung 2024 Table 1). Median age was 6.0 years (IQR 1.6-13.0) and
median body weight 20.1 kg (IQR 10.7-41.6); 58.3% were male; 72.8% White, 16.9%
Black or African American, 10.3% other, and 40.7% Hispanic or Latino. Median
bedside-Schwartz eGFR was 138.3 mL/min/1.73 m^2 (IQR 101.5-179.8), with 81.8% of
subjects in the normal renal-function stratum and only 1.0% severe - this is a
cohort whose renal function skews *high*, which is what makes the paper's
"adjust the initial dose upwards for supranormal eGFR" conclusion follow.

Patients on any form of extracorporeal therapy (ECMO, CRRT, extracorporeal liver
support), those with cardiopulmonary bypass within 7 days of starting
vancomycin, and those on chronic dialysis were excluded; a separate 33-subject
extracorporeal-therapy cohort was recruited but not analysed in this paper.
Dosing was chosen by the treating clinicians, giving actual regimens of 5-31
mg/kg q6h, most frequently 1-hour infusions of 15 mg/kg every 6-9 hours. A
median of 4 samples per subject (IQR 3-6) gave 1027 serum concentration records,
assayed by validated LC-MS/MS.

The same information is available programmatically via the model's `population`
metadata (`readModelDb("Jung_2024_vancomycin")()$population`).

# Source trace

Per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Jung_2024_vancomycin.R`. The table below
collects them in one place.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (CL at WT 20 kg, eGFR 141) | 2.763 L/h (RSE 2.746%) | Jung 2024 Table 2, clearance row |
| `lvc` (Vc at WT 20 kg) | 8.63 L (RSE 14.39%) | Jung 2024 Table 2, central volume row |
| `lq` (Q at WT 20 kg) | 3.808 L/h (RSE 20.45%) | Jung 2024 Table 2, intercompartmental clearance row |
| `lvp` (Vp at WT 20 kg) | 9.566 L (RSE 14.46%) | Jung 2024 Table 2, peripheral volume row |
| `e_wt_cl_q` | 0.75, fixed | Jung 2024 Table 2: `(WT/20)^0.75` printed on both CL and Q, no RSE; eTable 2 names the reference model "2-compartment model with allometric scaling" |
| `e_wt_vc_vp` | 1, fixed | Jung 2024 Table 2: `(WT/20)` printed on both Vc and Vp with no exponent, i.e. linear |
| `e_crcl_cl` | 0.5259 (RSE 14.92%) | Jung 2024 Table 2: `(eGFR/141)^0.5259` |
| Reference weight | 20 kg | Jung 2024 Table 2 and Results ("typical values for a subject weighing 20 kg") |
| Reference eGFR | 141 mL/min/1.73 m^2 | Jung 2024 Table 2 and Results ("CL 2.76 L/h for eGFR of 141 mL/min/1.73 m^2") |
| eGFR definition | `0.413 * height (cm) / Scr (mg/dL)` | Jung 2024 Table 1 footnote a (bedside Schwartz) |
| `etalcl` | omega^2 = 0.135678 | Jung 2024 Table 2: 38.12% CV (RSE 15.94%); supplement eMethods 4.1 gives `%CV = sqrt(exp(omega^2) - 1) * 100`, so omega^2 = log(CV^2 + 1) |
| `etalvc` | omega^2 = 0.719144 | Jung 2024 Table 2: 102.6% CV (RSE 28.64%); same conversion |
| BSV on Q, Vp | absent | Jung 2024 Table 2: reported as "NE" (not estimated) |
| `expSd` (`lnorm`) | 0.3310 | Jung 2024 Table 2: residual variability (log scale) variance 0.1096 (RSE 12.69%), sd 0.3310; supplement eMethods 4.2.4 specifies the log error model `log(Y) = log(Yhat) + eps` |
| Two-compartment linear-elimination ODEs | n/a | Jung 2024 Results and Conclusions: "A two-compartment model allometrically scaled with body weight with linear elimination" |
| IV infusion, no depot / bioavailability | n/a | Jung 2024 Methods: IV vancomycin, 1-hour infusions (Table 4 footnote) |
| AUC24 identity `Dose * (24/tau) / CL` | n/a | Jung 2024 Methods and supplement eMethods "Target Attainment Analysis" |
| AUC24-from-trough regression coefficients | 366, 54.9, 0.0894, 0.495 | Jung 2024 eTable 3 |

Two source details deserve emphasis because they are easy to get wrong and both
are resolved only by the supplement, not the main text:

* **The BSV column is an exact log-normal CV, not `100 * omega`.** eMethods 4.1
  states `%CV = sqrt(exp(omega^2) - 1) * 100`. Treating 102.6% as `100 * omega`
  would have set `omega^2 = 1.0527` instead of the correct 0.7191 - a 46%
  overstatement of the Vc variance.
* **The residual error is a log error model, not a proportional one.**
  eMethods 4.2.4 gives `log(Y) = log(Yhat) + eps`, which is nlmixr2's `lnorm()`,
  and the "Estimate" column of Table 2 is the *variance*
  (`sqrt(0.1096) = 0.3311`, matching the separately-printed 0.3310 sd).

# The packaged model reproduces Table 2 exactly

```{r table2-typical}
mod <- readModelDb("Jung_2024_vancomycin")

typical <- rxode2::rxSolve(
  mod,
  events = rxode2::et(amt = 0, cmt = "central") |> rxode2::et(0, cmt = "central"),
  params = c(WT = 20, CRCL = 141),
  omega  = NA,
  returnType = "data.frame"
)

table2 <- tibble::tibble(
  Parameter = c("Clearance (L/h)", "Central volume (L)",
                "Intercompartmental clearance (L/h)", "Peripheral volume (L)"),
  Published = c(2.763, 8.63, 3.808, 9.566),
  Model     = round(unlist(typical[1, c("cl", "vc", "q", "vp")]), 4)
)
stopifnot(max(abs(table2$Model - table2$Published)) < 1e-8)

knitr::kable(
  table2,
  caption = paste("Typical values at the Jung 2024 reference subject (WT 20 kg,",
                  "eGFR 141 mL/min/1.73 m^2) versus Jung 2024 Table 2."),
  align = c("l", "r", "r")
)
```

# Virtual cohort

Original subject-level data are not publicly available. Jung 2024 built its
virtual population of 2250 subjects "by randomly sampling with replacement the
baseline covariate vector of subjects from the study"; only the marginal median
and IQR of each covariate are published (Table 1). The cohort below therefore
uses log-normal marginals matched exactly to the published median and IQR of
body weight and eGFR, coupled by a seeded Latin-hypercube permutation, with the
two etas placed on their own exact quantiles. This is deterministic - no random
draw - so the numbers below are reproducible and cannot be moved by reseeding
(see Assumptions and deviations for the independence caveat).

```{r cohort}
n_per_arm <- 200L # cohort cap; 200/arm is ample and keeps the render fast
set.seed(20240908)

quantile_grid <- (seq_len(n_per_arm) - 0.5) / n_per_arm

# Log-normal quantile function pinned to a published median and IQR.
lognormal_at <- function(p, median, q25, q75) {
  sdlog <- (log(q75) - log(q25)) / (2 * stats::qnorm(0.75))
  median * exp(stats::qnorm(p) * sdlog)
}

base_cohort <- tibble::tibble(
  # Jung 2024 Table 1: WT median 20.1 kg (IQR 10.7-41.6)
  WT     = lognormal_at(quantile_grid[sample(n_per_arm)], 20.1, 10.7, 41.6),
  # Jung 2024 Table 1: eGFR median 138.3 (IQR 101.5-179.8) mL/min/1.73 m^2
  CRCL   = lognormal_at(quantile_grid[sample(n_per_arm)], 138.3, 101.5, 179.8),
  # etas on their own exact normal quantiles (omega^2 from the model file)
  etalcl = stats::qnorm(quantile_grid[sample(n_per_arm)], 0, sqrt(0.135678)),
  etalvc = stats::qnorm(quantile_grid[sample(n_per_arm)], 0, sqrt(0.719144))
)

tibble::tibble(
  Covariate = c("Body weight (kg)", "eGFR (mL/min/1.73 m^2)"),
  `Published median (IQR)` = c("20.1 (10.7-41.6)", "138.3 (101.5-179.8)"),
  `Cohort median (IQR)` = c(
    sprintf("%.1f (%.1f-%.1f)", median(base_cohort$WT),
            quantile(base_cohort$WT, 0.25), quantile(base_cohort$WT, 0.75)),
    sprintf("%.1f (%.1f-%.1f)", median(base_cohort$CRCL),
            quantile(base_cohort$CRCL, 0.25), quantile(base_cohort$CRCL, 0.75))
  )
) |>
  knitr::kable(caption = "Virtual cohort covariates versus Jung 2024 Table 1.")
```

# Simulation of the four dosing regimens

The four regimens are the ones Jung 2024 simulated (Table 4): 15 mg/kg q6h,
15 mg/kg q8h, 20 mg/kg q6h, and 30 mg/kg q8h, each as a 1-hour infusion
(Table 4 footnote). Steady state is imposed exactly with `ss = 1` on the first
dose rather than by simulating a burn-in, so the AUC over the 0-24 h window is
the steady-state AUC24 the paper's target-attainment analysis used. Each arm
gets a disjoint block of subject IDs.

```{r events}
regimens <- tibble::tibble(
  treatment = c("15 mg/kg q6h", "15 mg/kg q8h", "20 mg/kg q6h", "30 mg/kg q8h"),
  mg_per_kg = c(15, 15, 20, 30),
  tau       = c(6, 8, 6, 8),
  id_offset = c(0L, 200L, 400L, 600L)
)

make_arm <- function(mg_per_kg, tau, id_offset, treatment) {
  subj <- base_cohort |>
    dplyr::mutate(
      id        = id_offset + dplyr::row_number(),
      treatment = treatment,
      dose_mg   = mg_per_kg * WT
    )

  doses <- subj |>
    tidyr::expand_grid(dose_time = seq(0, 24 - tau, by = tau)) |>
    dplyr::transmute(
      id, treatment, WT, CRCL, etalcl, etalvc,
      time = dose_time,
      amt  = dose_mg,
      rate = dose_mg / 1, # 1-hour infusion (Jung 2024 Table 4 footnote)
      evid = 1L,
      cmt  = "central",
      # ss = 1 on the time-zero dose imposes exact steady state for the
      # q-tau regimen; later doses inside the window are ordinary doses.
      ss   = ifelse(dose_time == 0, 1L, 0L),
      ii   = ifelse(dose_time == 0, tau, 0)
    )

  # Observation grid: 0.25 h throughout, refined to 0.05 h over the infusion
  # and the immediate post-infusion peak so that trapezoidal / lin-up-log-down
  # AUC reproduces the analytic integral. Records are also placed immediately
  # BEFORE each dose: an observation sharing a dose's timestamp is evaluated
  # after that dose, so without these the interval minimum would be the
  # concentration a quarter-hour early rather than the true trough.
  peak_times <- as.numeric(outer(seq(0, 24 - tau, by = tau),
                                 seq(0, 1.5, by = 0.05), "+"))
  pre_dose_times <- seq(tau, 24, by = tau) - 1e-4
  grid <- sort(unique(c(seq(0, 24, by = 0.25), pre_dose_times,
                        peak_times[peak_times <= 24])))

  obs <- subj |>
    tidyr::expand_grid(time = grid) |>
    dplyr::transmute(
      id, treatment, WT, CRCL, etalcl, etalvc, time,
      amt = NA_real_, rate = 0, evid = 0L, cmt = "central", ss = 0L, ii = 0
    )

  dplyr::bind_rows(doses, obs) |> dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(lapply(
  seq_len(nrow(regimens)),
  function(i) {
    with(regimens[i, ], make_arm(mg_per_kg, tau, id_offset, treatment))
  }
))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
stopifnot(dplyr::n_distinct(events$id) == 4L * n_per_arm)
```

```{r simulate}
# etalcl / etalvc are supplied as columns, so omega = NA gives a fully
# deterministic solve (no resampling of the random effects).
sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("treatment", "WT", "CRCL"),
  omega  = NA,
  returnType = "data.frame"
)

stopifnot(!anyNA(sim$Cc), all(sim$Cc >= 0))
```

# Replicate published figures and findings

## Steady-state concentration-time profiles

Jung 2024 does not publish a typical-value concentration-time figure (its
eFigures 4-7 are diagnostics on observed data), so the panel below is
illustrative rather than a figure replication. It shows the biphasic decline the
paper describes, "with the second phase beginning at approximately 8 hours after
the end of infusion".

```{r profiles}
# rxSolve returns observation records only (no evid column), so no filter is
# needed here.
sim |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05), Q50 = median(Cc), Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~treatment) +
  scale_y_log10() +
  labs(
    x = "Time within the steady-state dosing interval window (h)",
    y = "Serum vancomycin (ug/mL)",
    title = "Steady-state profiles, median and 5th-95th percentile",
    caption = "Illustrative; Jung 2024 publishes no typical-value profile figure."
  )
```

## Clearance is a less-than-proportional function of eGFR

This is the paper's headline structural finding (Key Points: "As the eGFR
increased, the clearance of vancomycin showed a less than proportional
increase"). The exponent 0.5259 is what encodes it.

```{r cl-vs-egfr}
egfr_grid <- seq(20, 250, by = 1)
tibble::tibble(
  CRCL = egfr_grid,
  `Jung 2024 (exponent 0.5259)` = 2.763 * (egfr_grid / 141)^0.5259,
  `Proportional (exponent 1)`   = 2.763 * (egfr_grid / 141)
) |>
  tidyr::pivot_longer(-CRCL, names_to = "Scaling", values_to = "cl") |>
  ggplot(aes(CRCL, cl, linetype = Scaling)) +
  geom_line() +
  geom_vline(xintercept = 141, colour = "grey50") +
  labs(
    x = "Bedside-Schwartz eGFR (mL/min/1.73 m^2)",
    y = "Clearance at WT 20 kg (L/h)",
    title = "Less-than-proportional eGFR effect on vancomycin clearance",
    caption = "Vertical line: the 141 mL/min/1.73 m^2 reference. Jung 2024 Table 2."
  ) +
  theme(legend.position = "bottom")
```

# PKNCA validation

AUC24 is computed with PKNCA over the steady-state 0-24 h window. The dose
object declares `duration = 1` so PKNCA treats the records as 1-hour infusions
rather than bolus doses.

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# ss = 1 means the time-zero record already carries the steady-state trough,
# so a time-zero concentration exists for every subject; assert it rather than
# patching one in.
stopifnot(sum(sim_nca$time == 0) == 4L * n_per_arm)

dose_df <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | treatment + id,
  concu = "ug/mL", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | treatment + id,
  doseu = "mg", duration = 1
)

intervals <- data.frame(
  start = 0, end = 24,
  auclast = TRUE, cmax = TRUE, cmin = TRUE, cav = TRUE
)

nca_res <- suppressWarnings(
  PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
)

auc24 <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "auclast", start == 0, end == 24) |>
  dplyr::select(id, treatment, AUC24 = PPORRES)
stopifnot(nrow(auc24) == 4L * n_per_arm)
```

## Gate 1 - PKNCA AUC24 equals the paper's own analytic AUC24

Jung 2024 never integrates a simulated profile: both its target-attainment
analysis and its AUC24-versus-trough regression define
`AUC24 = Dose / CL * (24 / tau)` (Methods; supplement eMethods "Target
Attainment Analysis"). For a linear model at steady state that identity is
exact, so reproducing it from the ODE solution end-to-end validates the
structural model, the steady-state setup, and the NCA all at once.

```{r gate-analytic}
analytic <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::group_by(id, treatment) |>
  dplyr::summarise(dose_24h = sum(amt), .groups = "drop") |>
  dplyr::left_join(
    sim |> dplyr::group_by(id) |> dplyr::summarise(cl = dplyr::first(cl), .groups = "drop"),
    by = "id"
  ) |>
  dplyr::mutate(AUC24_analytic = dose_24h / cl)

gate1 <- auc24 |>
  dplyr::left_join(analytic, by = c("id", "treatment")) |>
  dplyr::mutate(rel_err = AUC24 / AUC24_analytic - 1)

max_rel_err <- max(abs(gate1$rel_err))
max_rel_err
stopifnot(max_rel_err < 0.005)
```

The worst per-subject discrepancy across all `r nrow(gate1)` subjects is
`r sprintf("%.3f%%", 100 * max_rel_err)`, which is the residual
lin-up/log-down interpolation error of the observation grid, not a model
discrepancy.

## Gate 2 - AUC24 is exactly proportional to the total daily dose

The four regimens deliver 60, 45, 80, and 90 mg/kg/day. Because the model is
linear and every arm shares one cohort, the AUC24 distributions must be exact
rescalings of one another - and Jung 2024's own eTable 4 medians are
(452, 339, 603, 678), whose ratios are (1, 0.750, 1.334, 1.500) against a daily
dose ratio of (1, 0.750, 1.333, 1.500). Making that a hard assertion turns a
published-table property into a regression guard.

```{r gate-linearity}
daily_dose <- regimens$mg_per_kg * 24 / regimens$tau
names(daily_dose) <- regimens$treatment

sim_median <- vapply(
  split(auc24$AUC24, auc24$treatment), median, numeric(1)
)[regimens$treatment]

gate2 <- max(abs(
  sim_median / sim_median[[1]] - daily_dose / daily_dose[[1]]
))
gate2
stopifnot(gate2 < 1e-4)
```

## Comparison against the published AUC24 distribution

Jung 2024 eTable 4 reports the median [IQR] of AUC24/MIC by regimen and sampling
strategy. At MIC = 1 ug/mL the ratio equals AUC24 itself, so the
infinite-sampling rows are directly comparable. `ncaComparisonTable()` pools the
simulated side by median, which is the statistic eTable 4 reports.

```{r nca_compare}
published <- tibble::tribble(
  ~treatment,      ~auclast,
  "15 mg/kg q6h",  452,
  "15 mg/kg q8h",  339,
  "20 mg/kg q6h",  603,
  "30 mg/kg q8h",  678
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res$result |> dplyr::filter(PPTESTCD == "auclast"),
  reference = published,
  by        = "treatment",
  params    = "auclast",
  units     = c(auclast = "ug*h/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Steady-state AUC24 versus Jung 2024 eTable 4",
                  "(MIC = 1 ug/mL, infinite-sampling rows).",
                  "* differs from reference by >20%."),
  align = c("l", "l", "r", "r", "r")
)
```

The full distribution, not just the median:

```{r auc-distribution}
auc24 |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    Simulated = sprintf("%.0f [%.0f, %.0f]", median(AUC24),
                        quantile(AUC24, 0.25), quantile(AUC24, 0.75)),
    .groups = "drop"
  ) |>
  dplyr::left_join(
    tibble::tibble(
      treatment = regimens$treatment,
      Published = c("452 [331, 630]", "339 [248, 473]",
                    "603 [442, 840]", "678 [497, 945]")
    ),
    by = "treatment"
  ) |>
  dplyr::rename(
    "Dosing regimen"                = treatment,
    "Published AUC24, median [IQR]" = Published,
    "Simulated AUC24, median [IQR]" = Simulated
  ) |>
  dplyr::relocate(`Published AUC24, median [IQR]`,
                  .before = `Simulated AUC24, median [IQR]`) |>
  knitr::kable(
    caption = paste("AUC24 (ug*h/mL) distribution versus Jung 2024 eTable 4,",
                    "MIC = 1, infinite sampling."),
    align = c("l", "r", "r")
  )
```

## Target attainment of AUC24/MIC >= 400 (Table 4)

```{r target-attainment}
attainment <- auc24 |>
  tidyr::expand_grid(MIC = c(0.5, 1, 2)) |>
  dplyr::group_by(treatment, MIC) |>
  dplyr::summarise(Simulated = 100 * mean(AUC24 / MIC >= 400), .groups = "drop")

published_ta <- tibble::tribble(
  ~treatment,     ~MIC, ~Published,
  "15 mg/kg q6h", 0.5,  96.2,
  "15 mg/kg q6h", 1.0,  60.4,
  "15 mg/kg q6h", 2.0,  13.2,
  "15 mg/kg q8h", 0.5,  86.8,
  "15 mg/kg q8h", 1.0,  36.5,
  "15 mg/kg q8h", 2.0,   5.1,
  "20 mg/kg q6h", 0.5,  99.4,
  "20 mg/kg q6h", 1.0,  80.8,
  "20 mg/kg q6h", 2.0,  27.8,
  "30 mg/kg q8h", 0.5,  99.8,
  "30 mg/kg q8h", 1.0,  86.8,
  "30 mg/kg q8h", 2.0,  36.5
)

ta_cmp <- published_ta |>
  dplyr::left_join(attainment, by = c("treatment", "MIC")) |>
  dplyr::mutate(
    Difference = round(Simulated - Published, 1),
    Simulated  = round(Simulated, 1)
  )

ta_cmp |>
  dplyr::rename(
    "Dosing regimen"          = treatment,
    "MIC (ug/mL)"             = MIC,
    "Published % >= 400"      = Published,
    "Simulated % >= 400"      = Simulated,
    "Difference (percentage points)" = Difference
  ) |>
  knitr::kable(
    caption = paste("Target attainment of AUC24/MIC >= 400 versus Jung 2024",
                    "Table 4 (infinite-sampling rows)."),
    align = c("l", "r", "r", "r", "r")
  )

max_ta_gap <- max(abs(ta_cmp$Difference))
max_ta_gap
```

Every simulated attainment is below the published value by
`r sprintf("%.1f", max_ta_gap)` percentage points or less, and the sign is
consistently negative. That is the expected consequence of the covariate
approximation: the virtual cohort's AUC24 IQR is slightly wider than the
paper's (see the distribution table above), so slightly more subjects fall below
the 400 threshold. The rank ordering of the regimens, the shape of the
MIC dependence, and the paper's clinical conclusion - that 15 mg/kg q6h and q8h
leave a large minority of critically ill children under target, while
20 mg/kg q6h and 30 mg/kg q8h exceed 80% at MIC = 1 - all reproduce.

```{r attainment-plot}
auc24 |>
  ggplot(aes(x = treatment, y = AUC24)) +
  geom_boxplot(outlier.size = 0.4) +
  geom_hline(yintercept = 400, linetype = "dashed") +
  scale_y_log10() +
  labs(
    x = NULL, y = "Steady-state AUC24 (ug*h/mL)",
    title = "AUC24 by regimen against the AUC24/MIC >= 400 target at MIC = 1",
    caption = "Dashed line: AUC24 = 400. Compare Jung 2024 Table 4 and eTable 4."
  ) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
```

# Independent cross-check: the AUC24-from-trough regression

Jung 2024's second product is a regression predicting AUC24 from a single
observed trough. The main text prints the equation as an undecoded image; its
coefficients are in eTable 3:

| eTable 3 parameter | Estimate | %RSE |
|---|---|---|
| Coefficient of trough concentration = 7 ug/mL | 366 | 1.98 |
| Additive shift for dose frequency < 8 hours | 54.9 | 24.2 |
| Slope of (DOSE [mg] - 500) | 0.0894 | 17.9 |
| Exponent of trough concentration | 0.495 | 3.92 |
| Residual variability (CV%) | 20.2 | 18.8 |

Those four coefficients reconstruct the equation as

```
AUC24 = (366 + 54.9 * DFRQ8 + 0.0894 * (Dose_mg - 500)) * (Ctrough / 7)^0.495
```

where `DFRQ8` is 1 when the dosing interval is under 8 hours. The reconstruction
is confirmed against the paper's own Table 3, which tabulates, for 12
regimen-by-weight combinations, the trough concentration that yields
AUC24/MIC = 400. Feeding those published troughs back through the equation must
return 400 for every row.

```{r regression-check}
auc24_from_trough <- function(dose_mg, interval_h, ctrough) {
  (366 + 54.9 * as.integer(interval_h < 8) + 0.0894 * (dose_mg - 500)) *
    (ctrough / 7)^0.495
}

table3 <- tibble::tribble(
  ~treatment,     ~tau, ~WT, ~ctrough,
  "15 mg/kg q6h",  6,   10,   7.4,
  "15 mg/kg q6h",  6,   20,   6.9,
  "15 mg/kg q6h",  6,   40,   6.1,
  "15 mg/kg q8h",  8,   10,  10.0,
  "15 mg/kg q8h",  8,   20,   9.3,
  "15 mg/kg q8h",  8,   40,   8.0,
  "20 mg/kg q6h",  6,   10,   7.2,
  "20 mg/kg q6h",  6,   20,   6.6,
  "20 mg/kg q6h",  6,   40,   5.6,
  "30 mg/kg q8h",  8,   10,   9.3,
  "30 mg/kg q8h",  8,   20,   8.0,
  "30 mg/kg q8h",  8,   40,   6.1
) |>
  dplyr::mutate(
    mg_per_kg = as.numeric(sub(" .*", "", treatment)),
    dose_mg   = mg_per_kg * WT,
    AUC24     = round(auc24_from_trough(dose_mg, tau, ctrough), 1)
  )

# Table 3 troughs are printed to one decimal, so the round trip lands near 400
# rather than exactly on it; 2 ug*h/mL is the rounding budget.
stopifnot(max(abs(table3$AUC24 - 400)) < 2)

table3 |>
  dplyr::select(treatment, WT, ctrough, AUC24) |>
  dplyr::rename(
    "Dosing regimen"               = treatment,
    "Body weight (kg)"             = WT,
    "Published trough (ug/mL)"     = ctrough,
    "Regression AUC24 (ug*h/mL)"   = AUC24
  ) |>
  knitr::kable(
    caption = paste("Jung 2024 Table 3 round trip: the published trough",
                    "achieving AUC24/MIC = 400 at MIC = 1, fed back through the",
                    "eTable 3 regression, returns 400."),
    align = c("l", "r", "r", "r")
  )
```

The two products can now be checked against each other over the whole virtual
cohort, which is how the regression was actually built: each subject's
steady-state trough (PKNCA `cmin` over the dosing window) is fed through the
regression, and the result compared with that subject's own model AUC24. No
extra simulation is needed - both quantities come from the PKNCA run above.

```{r cross-check}
cross <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "cmin", start == 0, end == 24) |>
  dplyr::select(id, treatment, ctrough = PPORRES) |>
  dplyr::left_join(auc24, by = c("id", "treatment")) |>
  dplyr::left_join(
    events |> dplyr::filter(evid == 1) |>
      dplyr::group_by(id) |>
      dplyr::summarise(dose_mg = dplyr::first(amt),
                       tau = 24 / dplyr::n(), .groups = "drop"),
    by = "id"
  ) |>
  dplyr::mutate(
    AUC24_regression = auc24_from_trough(dose_mg, tau, ctrough),
    pct_diff = 100 * (AUC24_regression / AUC24 - 1)
  )

pooled_bias <- median(cross$pct_diff)
pooled_cv <- 100 * sd(cross$AUC24_regression / cross$AUC24) /
  mean(cross$AUC24_regression / cross$AUC24)

cross |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    `Median % difference` = round(median(pct_diff), 1),
    `IQR of % difference` = sprintf("%.1f to %.1f", quantile(pct_diff, 0.25),
                                    quantile(pct_diff, 0.75)),
    .groups = "drop"
  ) |>
  dplyr::rename("Dosing regimen" = treatment) |>
  knitr::kable(
    caption = paste("Regression AUC24 (from each subject's simulated",
                    "steady-state trough) versus that subject's model AUC24,",
                    "by regimen."),
    align = c("l", "r", "r")
  )
```

Three things fall out, and all three are properties of the paper rather than of
this extraction:

* **The reconstruction is unbiased in aggregate.** Pooled across the cohort the
  median difference is `r sprintf("%+.1f%%", pooled_bias)`, matching Jung 2024's
  own claim that the equation "provided an unbiased prediction of AUC24 across
  the observed range of trough vancomycin concentrations."
* **Its scatter matches the published residual error.** The ratio of the two
  AUC24 routes has a CV of `r sprintf("%.1f%%", pooled_cv)`, against the 20.2%
  residual variability reported for the regression in eTable 3. Recovering that
  number independently is a third confirmation - alongside the Table 3 round
  trip and the unbiasedness - that the eTable 3 coefficients were reassembled
  into the right functional form.
* **But the bias is regimen-dependent, and worst where the paper extrapolates.**
  The per-regimen medians run from about +12% for 15 mg/kg q6h to about -21% for
  30 mg/kg q8h. The regression carries dose as a *linear additive* term centred
  at 500 mg and was fit to the observed regimens only (5-31 mg/kg q6h, most
  commonly 15 mg/kg every 6-9 hours). 30 mg/kg q8h is a purely simulated
  regimen, reaching 1200 mg per dose in the heaviest children - well outside
  that calibration range - so the linear dose term is being extrapolated. A user
  who wants AUC24 for one of the higher simulated regimens should take it from
  the PK model (`Dose * (24/tau) / CL`), not from the trough regression.

# Assumptions and deviations

* **Covariate joint distribution.** Only marginal medians and IQRs are published
  (Table 1), so the virtual cohort uses log-normal marginals matched to those
  and treats body weight and eGFR as **independent**. The real cohort almost
  certainly has some weight-eGFR association. Because AUC24 scales as
  `WT^0.25 * eGFR^-0.5259`, a positive association would partly cancel and
  narrow the AUC24 distribution - which is the direction of the residual
  disagreement observed above (simulated IQR slightly wider, target attainment
  4-6 percentage points lower than published). No parameter was adjusted to
  close that gap.
* **Cohort size.** 200 subjects per arm rather than the paper's 2250, per the
  repository cap. The comparison statistics are medians, IQRs, and proportions,
  all of which are stable at n = 200; the cohort is deterministic (exact
  marginal quantiles, seeded coupling) so the reported numbers do not move
  under reseeding.
* **Time-varying covariates held constant.** Jung 2024 treats both body weight
  and eGFR as time-varying (eTable 1) and fit them that way. The target-
  attainment simulation resampled *baseline* covariate vectors, so this vignette
  holds both constant over the 24-hour steady-state window, matching the paper's
  own simulation rather than its estimation dataset.
* **Steady state imposed, not simulated.** `ss = 1` replaces a multi-day burn-in.
  This matches the paper's analytic `Dose * (24/tau) / CL` definition of AUC24
  exactly (Gate 1) and avoids the long-half-life tail that the 102.6% CV on Vc
  would otherwise produce in a finite burn-in.
* **Body weight is not the dosing basis in the model file.** The model takes
  `amt` in mg; the mg/kg regimens are converted to mg in the vignette's event
  table using each subject's `WT`.
* **eGFR equation.** The packaged model uses the bedside Schwartz eGFR that
  Jung 2024 Table 2 reports. The paper also refit the model with the CKiD
  cystatin-C equation in the 248-subject subset having both serum creatinine and
  cystatin C, and found equivalent predictive performance (eFigures 8-9), but no
  separate parameter table is published for that fit, so it is not packaged as a
  second model.
* **The AUC24-from-trough regression is not packaged as a model file.** It is a
  nonlinear regression on model-predicted AUC24, not an ODE model, so it does
  not belong in `inst/modeldb/`. Its coefficients are reproduced in this
  vignette from eTable 3 and validated three independent ways (Table 3 round
  trip, aggregate unbiasedness, and recovery of the published 20.2% residual CV).
* **The trough regression should not be used for the higher simulated
  regimens.** Its dose term is linear and centred at 500 mg and it was fit only
  to the observed regimens, so it disagrees with the PK model by a median of
  about -21% at 30 mg/kg q8h (see the cross-check section). This is a property
  of Jung 2024, not of the packaged model, and no parameter was altered because
  of it; it is recorded here so a downstream user does not mistake the
  regression for a general-purpose AUC24 estimator.
* **The extracorporeal-therapy cohort is out of scope.** Jung 2024 recruited 33
  additional subjects on ECMO / CRRT / extracorporeal liver support and
  explicitly excluded them from this analysis (eMethods); no model was published
  for them.
* **No non-paper-derived parameter values.** Every value in `ini()` comes from
  Jung 2024 Table 2 or the supplement's eMethods; nothing was digitised from a
  figure, obtained by correspondence, or carried from another model. The
  supplement was retrieved from the publisher's own permalink
  (<http://links.lww.com/CCX/B406>) and is archived alongside the lead PDF.
