nlmixr2 and tidyvpc
By Matthew Fidler in nlmixr2
February 27, 2026
VPCs – a new method for VPCs in nlmixr2
I am excited to announce a new VPC method available in nlmixr2 –
tidyvpc.
This method can be built on-top of previous code and provides slightly
different visualizations than the vpc package.
Below is an example that shows the two different methods:
library(nlmixr2)## ── Attaching packages ───────────────────────────────────────────────────────── nlmixr2 5.0.0 ──## ✔ lotri 1.0.2 ✔ monolix2rx 0.0.6
## ✔ nlmixr2data 2.0.9 ✔ nlmixr2lib 0.3.2
## ✔ nlmixr2est 5.0.2 ✔ nlmixr2rpt 0.2.2
## ✔ nlmixr2extra 5.0.0 ✔ nonmem2rx 0.1.9
## ✔ nlmixr2plot 5.0.1 ✔ posologyr 1.2.8
## ✔ rxode2 5.0.1 ✔ shinyMixR 0.5.3
## ✔ babelmixr2 0.1.11 ✔ xpose.nlmixr2 0.4.1
## ✔ ggPMX 1.3.2## ── Optional Packages Loaded/Ignored ─────────────────────────────────────────── nlmixr2 5.0.0 ──## ✔ babelmixr2 ✔ nonmem2rx
## ✔ ggPMX ✔ posologyr
## ✔ monolix2rx ✔ shinyMixR
## ✔ nlmixr2lib ✔ xpose.nlmixr2
## ✔ nlmixr2rpt## ── Conflicts ───────────────────────────────────────────────────────────── nlmixr2conflicts() ──
## ✖ rxode2::boxCox() masks nlmixr2est::boxCox()
## ✖ rxode2::yeoJohnson() masks nlmixr2est::yeoJohnson()library(patchwork)
pheno <- function() {
ini({
tcl <- log(0.008) # typical value of clearance
tv <- log(0.6) # typical value of volume
## var(eta.cl)
eta.cl + eta.v ~ c(1,
0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
# interindividual variability on clearance and volume
add.err <- 0.1 # residual variability
})
model({
cl <- exp(tcl + eta.cl) # individual value of clearance
v <- exp(tv + eta.v) # individual value of volume
ke <- cl / v # elimination rate constant
d/dt(A1) = - ke * A1 # model differential equation
cp = A1 / v # concentration in plasma
cp ~ add(add.err) # define error model
})
}
fit <- nlmixr(pheno, pheno_sd, "saem",
control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))## ℹ parameter labels from comments are typically ignored in non-interactive mode## ℹ Need to run with the source intact to parse comments## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## ℹ calculate uninformed etas## ℹ done## Calculating covariance matrix## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem predOnly model 0...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 1...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 2...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## → Calculating residuals/tables## → loading into symengine environment...## → pruning branches (`if`/`else`) of full model...## ✔ done## → calculate jacobian## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate sensitivities## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate ∂(f)/∂(η)## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate ∂(R²)/∂(η)## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in inner model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in inner model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00
##
## [====|====|====|====|====|====|====|====|====|====] 0:00:00
##
## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in EBE model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in EBE model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → compiling inner model...## ✔ done## → finding duplicate expressions in FD model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in FD model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → compiling EBE model...## ✔ done## → compiling events FD model...## ✔ done## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## → loading into symengine environment...## → pruning branches (`if`/`else`) of full model...## ✔ done## → calculate jacobian## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate sensitivities## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate ∂(f)/∂(η)## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → calculate ∂(R²)/∂(η)## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in inner model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in inner model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00
##
## [====|====|====|====|====|====|====|====|====|====] 0:00:00
##
## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in EBE model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in EBE model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → compiling inner model...## ✔ done## → finding duplicate expressions in FD model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in FD model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → compiling EBE model...## ✔ done## → compiling events FD model...## ✔ done# Using the VPC package, a regular VPC
p1a <- vpcPlot(fit, show=list(obs_dv=TRUE))## [====|====|====|====|====|====|====|====|====|====] 0:00:00## Warning in filter_dv(obs, verbose): No software packages matched for filtering values, not filtering.
## Object class: other, data.frame
## Available filters: phoenix, nonmem## Warning in filter_dv(sim, verbose): No software packages matched for filtering values, not filtering.
## Object class: other, nlmixr2vpcSim, data.frame
## Available filters: phoenix, nonmem## Registered S3 method overwritten by 'classInt':
## method from
## print.classIntervals ggPMX## A prediction-corrected VPC
p2a <- vpcPlot(fit, pred_corr = TRUE, show=list(obs_dv=TRUE))## [====|====|====|====|====|====|====|====|====|====] 0:00:00## Warning in filter_dv(obs, verbose): No software packages matched for filtering values, not filtering.
## Object class: other, data.frame
## Available filters: phoenix, nonmem
## Warning in filter_dv(obs, verbose): No software packages matched for filtering values, not filtering.
## Object class: other, nlmixr2vpcSim, data.frame
## Available filters: phoenix, nonmem# using the tidyvpc package with method="tidyvpc"
# Using the VPC package, a regular VPC
p1b <- vpcPlot(fit, method="tidyvpc")## [====|====|====|====|====|====|====|====|====|====] 0:00:00## A prediction-corrected VPC
p2b <- vpcPlot(fit, pred_corr = TRUE, method="tidyvpc")## [====|====|====|====|====|====|====|====|====|====] 0:00:00(p1a * p1b) / (p2a * p2b)
On the left is the vpc package and on the right is the tidyvpc package.
What inspired adding tidyvpc package now?
The addition of tidyvpc was driven by a CRAN request to fix
nlmixr2plot. These errors started when the vpc package was
updated. For people experiencing similar errors, using the new tidyvpc
method may be helpful.
Still, I have submitted a vpc bug-fix which also fixes all the
issues we had and will hopefully be integrated in the next CRAN
release of vpc.
If you want to try out the bug-fix I submitted before it is released use:
remotes::install_github("mattfidler/vpc")