using Pumas
using PharmaDatasets
using DataFramesMeta
using CairoMakie
using AlgebraOfGraphics
using PumasUtilities

Visual Predictive Checks for Continuous Data Models
1 Learning Objectives
By the end of this tutorial, you will be able to generate and interpret VPC plots in Pumas, explore misspecification via simulation, and apply stratification and prediction correction.
2 Introduction
The original Visual Predictive Check was proposed to compare the observed data to prediction intervals from an estimated model. If the model is not in agreement with the data we will observe data curves outside of the prediction intervals and we may reject the model. In a Bayesian context, a similar plot is used to verify the posterior predictive distribution against the data.
In practice, the interpretation is not so easy. Covariates including different treatments, biologically significant differences between patients, and different sample times complicates matters significantly. This means that the clear visual check quickly becomes muddied. In this tutorial, we will see how to produce a VPC for continuous data in Pumas and we will try to adjust for some of the complications involved with using the VPC in practice.
First, we set up a Warfarin PK model.
= @model begin
warfarin_pk_model
@metadata begin
= "Warfarin 1-compartment PK model (PD removed)"
desc = u"hr"
timeu end
@param begin
# PK parameters
"""
Clearance (L/hr)
"""
∈ RealDomain(lower = 0.0, init = 0.134)
pop_CL """
Central volume (L)
"""
∈ RealDomain(lower = 0.0, init = 8.11)
pop_Vc """
Absorption lag time (hr)
"""
∈ RealDomain(lower = 0.0, init = 0.523)
pop_tabs
# Inter-individual variability
"""
- ΩCL: Clearance
- ΩVc: Central volume
"""
∈ PDiagDomain([0.01, 0.01])
pk_Ω # Residual variability
"""
σ_prop: Proportional error
"""
∈ RealDomain(lower = 0.0, init = 0.00752)
σ_prop end
@random begin
~ MvNormal(pk_Ω) # mean = 0, covariance = pk_Ω
pk_η end
@covariates begin
"""
FSZCL: Clearance scaling factor
"""
FSZCL"""
FSZV: Volume scaling factor
"""
FSZVend
@pre begin
= FSZCL * pop_CL * exp(pk_η[1])
CL = FSZV * pop_Vc * exp(pk_η[2])
Vc
= pop_tabs
tabs = log(2) / tabs
Ka end
@vars begin
:= Central / Vc
cp end
@dynamics begin
' = -Ka * Depot
Depot' = Ka * Depot - (CL / Vc) * Central
Centralend
@derived begin
"""
Concentration (ng/mL)
"""
~ @. Normal(cp, sqrt((σ_prop * cp)^2))
conc end
end
PumasModel
Parameters: pop_CL, pop_Vc, pop_tabs, pk_Ω, σ_prop
Random effects: pk_η
Covariates: FSZCL, FSZV
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived: conc
Observed: conc
# Transform the data in a single chain of operations
= @chain dataset("pumas/warfarin_pumas") begin
warfarin_data_wide @rtransform begin
# Scaling factors
:FSZV = :wtbl / 70 # volume scaling
:FSZCL = (:wtbl / 70)^0.75 # clearance scaling (allometric)
end
end
Row | id | time | evid | amt | cmt | conc | pca | wtbl | age | sex | FSZV | FSZCL |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Int64 | Float64 | Int64 | Float64? | Int64? | Float64? | Float64? | Float64 | Int64 | String1 | Float64 | Float64 | |
1 | 1 | 0.0 | 1 | 100.0 | 1 | missing | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
2 | 1 | 0.5 | 0 | missing | missing | 0.0 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
3 | 1 | 1.0 | 0 | missing | missing | 1.9 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
4 | 1 | 2.0 | 0 | missing | missing | 3.3 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
5 | 1 | 3.0 | 0 | missing | missing | 6.6 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
6 | 1 | 6.0 | 0 | missing | missing | 9.1 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
7 | 1 | 9.0 | 0 | missing | missing | 10.8 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
8 | 1 | 12.0 | 0 | missing | missing | 8.6 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
9 | 1 | 24.0 | 0 | missing | missing | 5.6 | 44.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
10 | 1 | 36.0 | 0 | missing | missing | 4.0 | 27.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
11 | 1 | 48.0 | 0 | missing | missing | 2.7 | 28.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
12 | 1 | 72.0 | 0 | missing | missing | 0.8 | 31.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
13 | 1 | 96.0 | 0 | missing | missing | missing | 60.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
319 | 32 | 48.0 | 0 | missing | missing | 6.9 | 24.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
320 | 32 | 72.0 | 0 | missing | missing | 4.4 | 23.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
321 | 32 | 96.0 | 0 | missing | missing | 3.5 | 20.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
322 | 32 | 120.0 | 0 | missing | missing | 2.5 | 22.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
323 | 33 | 0.0 | 1 | 100.0 | 1 | missing | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
324 | 33 | 0.0 | 0 | missing | missing | missing | 100.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
325 | 33 | 24.0 | 0 | missing | missing | 9.2 | 49.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
326 | 33 | 36.0 | 0 | missing | missing | 8.5 | 32.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
327 | 33 | 48.0 | 0 | missing | missing | 6.4 | 26.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
328 | 33 | 72.0 | 0 | missing | missing | 4.8 | 22.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
329 | 33 | 96.0 | 0 | missing | missing | 3.1 | 28.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
330 | 33 | 120.0 | 0 | missing | missing | 2.5 | 33.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
= read_pumas(
pop
warfarin_data_wide;= :id,
id = :time,
time = :amt,
amt = :cmt,
cmt = :evid,
evid = [:sex, :wtbl, :FSZV, :FSZCL],
covariates = [:conc],
observations )
Population
Subjects: 32
Covariates: sex, wtbl, FSZV, FSZCL
Observations: conc
= fit(warfarin_pk_model, pop, init_params(warfarin_pk_model), FOCE()) foce_fit
[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 9.744619e+04 1.950953e+05 * time: 0.04359602928161621 1 1.313047e+04 2.609679e+04 * time: 1.8302209377288818 2 9.658030e+03 1.909689e+04 * time: 1.8672621250152588 3 4.233250e+03 8.101602e+03 * time: 2.1364200115203857 4 2.367624e+03 4.268593e+03 * time: 2.172765016555786 5 1.313943e+03 2.054844e+03 * time: 2.2109858989715576 6 8.374399e+02 1.010773e+03 * time: 2.253291130065918 7 6.115489e+02 4.863561e+02 * time: 2.296687126159668 8 5.146179e+02 2.031565e+02 * time: 2.338041067123413 9 4.860731e+02 7.536873e+01 * time: 2.3803749084472656 10 4.801210e+02 5.238009e+01 * time: 2.421694040298462 11 4.795642e+02 5.405505e+01 * time: 2.4608869552612305 12 4.795291e+02 5.414517e+01 * time: 2.49833607673645 13 4.794853e+02 5.383276e+01 * time: 2.5354840755462646 14 4.793492e+02 5.211648e+01 * time: 2.5733790397644043 15 4.790578e+02 4.736685e+01 * time: 2.6153299808502197 16 4.784683e+02 3.589484e+01 * time: 2.6899759769439697 17 4.776996e+02 1.930725e+01 * time: 2.761275053024292 18 4.772243e+02 1.316223e+01 * time: 2.7976839542388916 19 4.770920e+02 1.346641e+01 * time: 2.8353099822998047 20 4.770711e+02 1.365789e+01 * time: 2.872580051422119 21 4.770686e+02 1.368452e+01 * time: 2.908911943435669 22 4.770612e+02 1.372786e+01 * time: 2.94673490524292 23 4.770441e+02 1.378138e+01 * time: 2.9840610027313232 24 4.770002e+02 1.384496e+01 * time: 3.0224690437316895 25 4.769031e+02 1.387480e+01 * time: 3.0620369911193848 26 4.767248e+02 2.003031e+01 * time: 3.1024129390716553 27 4.765128e+02 3.331465e+01 * time: 3.1424739360809326 28 4.763904e+02 3.753509e+01 * time: 3.1818110942840576 29 4.763566e+02 3.512805e+01 * time: 3.2200989723205566 30 4.763399e+02 3.256639e+01 * time: 3.2580010890960693 31 4.763051e+02 2.857696e+01 * time: 3.2957959175109863 32 4.762186e+02 2.633532e+01 * time: 3.334536075592041 33 4.759910e+02 2.972803e+01 * time: 3.371752977371216 34 4.754039e+02 3.513095e+01 * time: 3.4115920066833496 35 4.739900e+02 4.343181e+01 * time: 3.541343927383423 36 4.718772e+02 5.172641e+01 * time: 3.5806899070739746 37 4.699605e+02 5.109410e+01 * time: 3.618601083755493 38 4.679301e+02 1.837720e+01 * time: 3.656856060028076 39 4.677728e+02 9.841366e+00 * time: 3.6952130794525146 40 4.677644e+02 9.868356e+00 * time: 3.730889081954956 41 4.677619e+02 9.877285e+00 * time: 3.765291929244995 42 4.677607e+02 9.881556e+00 * time: 3.7983651161193848 43 4.677536e+02 9.897820e+00 * time: 3.8326659202575684 44 4.677390e+02 9.912120e+00 * time: 3.867448091506958 45 4.676965e+02 9.918954e+00 * time: 3.9025161266326904 46 4.675898e+02 1.060357e+01 * time: 3.9380240440368652 47 4.673112e+02 1.524992e+01 * time: 3.9743709564208984 48 4.666264e+02 2.323583e+01 * time: 4.013211965560913 49 4.649869e+02 2.870754e+01 * time: 4.053476095199585 50 4.638128e+02 3.688505e+01 * time: 4.0930869579315186 51 4.617598e+02 3.172279e+01 * time: 4.132256031036377 52 4.609261e+02 2.389881e+01 * time: 4.169701099395752 53 4.601257e+02 7.743696e+00 * time: 4.2052929401397705 54 4.600344e+02 2.839021e+00 * time: 4.240345001220703 55 4.600202e+02 6.791380e-01 * time: 4.274462938308716 56 4.600185e+02 1.392758e-01 * time: 4.307496070861816 57 4.600183e+02 3.003773e-02 * time: 4.339316129684448 58 4.600183e+02 4.113068e-03 * time: 4.372266054153442 59 4.600183e+02 3.717389e-03 * time: 4.406907081604004 60 4.600183e+02 3.347682e-03 * time: 4.441642999649048 61 4.600183e+02 3.347291e-03 * time: 4.482463121414185 62 4.600183e+02 3.347291e-03 * time: 4.534295082092285 63 4.600183e+02 2.864074e-03 * time: 4.5990989208221436 64 4.600183e+02 2.864074e-03 * time: 4.645575046539307 65 4.600183e+02 2.864074e-03 * time: 4.692521095275879 66 4.600183e+02 2.864074e-03 * time: 4.740606069564819 67 4.600183e+02 2.864074e-03 * time: 4.788588047027588 68 4.600183e+02 2.864074e-03 * time: 4.838531017303467 69 4.600183e+02 2.864074e-03 * time: 4.890316009521484 70 4.600183e+02 2.864074e-03 * time: 4.940639972686768
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
0 6
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: NoObjectiveChange
Log-likelihood value: -460.01832
--------------------
Estimate
--------------------
pop_CL 0.13782
pop_Vc 8.4196
pop_tabs 1.0986
pk_Ω₁,₁ 0.056591
pk_Ω₂,₂ 0.013582
σ_prop 0.23567
--------------------
3 What exactly is a VPC?
A Visual Predictive Check (VPC) is a plot constructed to identify misspecification of the model. It attempts to compare the distribution of the data with the distribution of simulations from the estimated model. Essentially constructed in parts:
flowchart TD Start(["Start"]) --> Q1["For pcVPC calculate the prediction median function<br>and apply it to the observed independent variable."] Q1 --> Q2["Calculate the sequence of quantiles (0.1, 0.5, 0.9)<br>for the (corrected) dependent variable observations<br>across the independent variable nodes"] Q2 --> Q3 Q2 --> Q6 Q3["Simulate all datasets using estimated parameters with<br>random effects and error model sampling (and apply<br>correction for pcVPC)"] --> Q4["Calculate quantiles for each simulated dataset"] Q4 --> Q5["Determine prediction intervals across<br>simulations per bin/location"] Q5 --> Q6["Plot results (observed vs simulated intervals)"] Q6 --> Q7["Check if observed medians lie within intervals"] Q7 --> End(["End"]) %% Labels and styles subgraph n1 ["VPC Workflow"] end classDef step fill:#9ecfdf,stroke:#000,stroke-width:1px,color:#000 classDef Peach stroke-width:1px,stroke-dasharray:none,stroke:#FBB35A,fill:#FFEFDB,color:#8F632D class Q1,Q2,Q3,Q4,Q5,Q6,Q7 step class n1 Peach
This differs from other diagnostics such as weighted residuals by considering not only data compared to the modeled mean, but tries to make sure that more or less extreme subjects are also well-predicted by the model.
4 Applying the VPC workflow in Pumas
In Pumas we have a predefined workflow for simulating the required data as well as summarizing the data using the binless VPC method. In short, the binless method avoids binning data across the independent variable. Instead, it applies a smoothed local quantile regression to combine information from data that is close together according to the kernel used. Binless VPCs are advantageous over traditional VPCs as they allow smoother representation of model-data agreement rather than across few discrete intervals.
The workflow consists of just two functions:
vpc(...)
to simulate the required data and calculate the quantities used for lines and ribbons in the final plotsvpc_plot(...)
fromPumasUtilities
that takes the data and presents it in a plot
simobs
It is also possible to simulate the data using simobs
but we will not cover that approach in this tutorial.
4.1 How to simulate?
To simulate, we use the vpc
function along with input on how the simulations and quantile calculations should be performed. The most simple invocation of the function is:
= vpc(foce_fit) myvpc
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 499
Subjects in data: 32
Stratification variable(s): None
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
This creates an object that has all the information we’ll need. Here’s a brief explanation for each line of the VPC output:
Visual Predictive Check:
This is the name of the diagnostic tool used to assess how well your model predicts observed data by comparing simulated data distributions.Type of VPC: Continuous VPC:
The VPC was performed for a dependent variable (e.g., concentration) that changes continuously over time.Simulated populations: 499:
The model generated 499 synthetic datasets to create prediction intervals for comparison.Subjects in data: 31:
Your original dataset contains 31 individual subjects.Stratification variable(s): None:
The VPC was performed on the entire dataset without splitting it into subgroups (e.g., by sex or age).Confidence level: 0.95:
The shaded prediction intervals on the plot are calculated to contain 95% of the simulated data.VPC lines: quantiles ([0.1, 0.5, 0.9]):
The VPC plot will display the 10th percentile (0.1), median (0.5), and 90th percentile (0.9) of the simulated data.
While the above function call will produce data for a VPC we often need to make adjustments in practice. Consider the (partial) signature where we omitted some options for clarity of this tutorial:
vpc(
::AbstractFittedPumasModel;
fpm::Integer = 499,
samples::Vector{Symbol} = [keys(fpm.data[1].observations)[1]],
observations::Vector{Symbol} = [:time],
covariates::Vector{Symbol} = Symbol[],
stratify_by= nothing,
bandwidth ::Bool = false,
prediction_correction )
Some important potential keyword arguments in the vpc
function are:
samples
(default is 500) to control the number times the population is re-simulated according to the fitted modelobservations
is a vector of symbols (e.g.[:conc]
) to control what endpoint to perform the VPC for if there are multiple present in the modelcovariates
(default is:time
) is a vector of symbols (e.g.[:tad]
) to control the first axis variablestratify_by
to control stratification of VPC calculations and plottingbandwidth
is a number (default is 2.0 for continuous VPCs) to control the level of smoothingprediction_correction
(default is false) to apply a type of correction for differences in subjects. If prediction correction is applied it is currently not possible to also use stratification.
Please see the docstring (?vpc
) for all options and defaults.
With the output of vpc
in hand, we can continue to plotting the VPC.
4.2 How to plot?
Plotting is done with the vpc_plot
function. This function accepts the output of vpc()
and can be customized in the same way as the rest of the built-in plotting functions covered earlier.
vpc_plot(myvpc, figurelegend = (; position = :b, nbanks = 2, tellheight = true))
The VPC looks good. The data medians are generally well within the ribbons.
4.3 A misspecified model
To produce a VPC that does not look good we can intentionally constrain the model to not produce the estimates found above. Let us test whether forcing a higher pop_CL
can be detected as misspecification in our setup:
= fit(
foce_fit_CL
warfarin_pk_model,
pop,init_params(warfarin_pk_model)..., pop_CL = 0.15),
(FOCE();
= (:pop_CL,),
constantcoef = (; show_trace = false),
optim_options )
[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
1 5
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -461.70032
----------------------
Estimate
----------------------
† pop_CL 0.15
pop_Vc 8.4758
pop_tabs 1.0994
pk_Ω₁,₁ 0.063764
pk_Ω₂,₂ 0.013783
σ_prop 0.23778
----------------------
† indicates constant parameters
Then we perform the required calculations
= vpc(foce_fit_CL) myvpc_CL
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 499
Subjects in data: 32
Stratification variable(s): None
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
And produce the plot
vpc_plot(myvpc_CL, figurelegend = (; position = :b, nbanks = 2, tellheight = true))
In this case, the model does not look well-specified, and we should consider not fixing the population level clearance parameter. Note, that the interpretation is up to the person looking at the plot. In this case, it appears as if the fit is decent for early time points and bad for later time points. This would tend to rule out a problem in the absorption phase. The question is then if it’s misspecified because of clearance or maybe we didn’t correctly account for a distribution phase with many peripheral compartments? In this case, we introduced the problem ourselves so we have a good idea, but in general it’s not always easy.
4.4 Stratification
If a covariate with a limited number of values is thought to be significant in the model or if we want to make sure that the fit is equally good across strata we should produce a stratified plot.
= vpc(foce_fit; stratify_by = [:sex]) my_vpc_sex
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 499
Subjects in data: 32
Stratification variable(s): [:sex]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
vpc_plot(my_vpc_sex; figurelegend = (; position = :b, nbanks = 2, tellheight = true))
We see that this stratified plot looks good as well.
4.5 Bandwidth selection
Let us try to experiment with the bandwidth. The bandwidth controls how wide of a window we will use to construct the local regression lines. If we consider a time point on the first axis, the quantile lines are constructed in a way where local data points carry more weight in the construction of the line than data points further away. The bandwidth
keyword argument controls how far away to put significant weight. In the dataset we observe that there are some very large observations at :time == 2.0
. Let us increase the bandwidth such that when we construct the value for the very early time points we end up putting too much weight to these later time points.
= vpc(foce_fit; stratify_by = [:sex], bandwidth = 5) my_vpc_sex
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 499
Subjects in data: 32
Stratification variable(s): [:sex]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
vpc_plot(my_vpc_sex; figurelegend = (; position = :b, nbanks = 2, tellheight = true))
According to this VPC, the data median starts at a large positive number at time point 0! This leads to a model misspecification according to the red outlier band. When outliers are observed we need to consider at least three possible scenarios:
- Is the data wrong? There could be unphysiological values that we need to explain and that our model could not be expected to explain. Hopefully, exploratory data analysis (EDA) caught such values.
- Is the model wrong? If the data looks good, maybe the model simply cannot fit these data points. We may need to take a closer look at the model fit for specific subjects that contribute to this upper part of the observed distribution.
- Is the VPC misspecified? Just as binned VPCs require appropriate binning, the binless VPC has aggregation of similar data points by use of local regression. We can try to vary the
bandwidth
parameter to see if deviating from the default value of2
However, in the original plot we saw that the black curve went down to a small value at :time == 0.0
as expected. The problem highlights something that is true for all VPCs and should always be kept in mind. We are aggregating data across time and covariate differences, so it is not always easy to know what the VPC curve “should” look like when the scenarios become complex. Here, we felt comfortable concluding that the first VPC did not invalidate our fit because we don’t think any subjects should start with positive values of the drug in plasma, but it may not always be easy to say if the VPC or the fit is the reason for an apparent indication of misspecification.
4.6 Changing to Time-After-Dose (tad
)
When producing a VPC it is often suggested to plot time-after-dose as the first axis variable and stratify by the number of the dose unless the dosage regimen is more complicated and several doses are administered at the same time or very close together. To change the variable on the first axis, we use the covariates
keyword. Stratification is done as before.
= vpc(foce_fit; covariates = [:tad], stratify_by = [:dosenum]) my_vpc_sex
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 499
Subjects in data: 32
Stratification variable(s): [:dosenum]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
vpc_plot(my_vpc_sex; figurelegend = (; position = :b, nbanks = 2, tellheight = true))
In our example, there is not real change because we only have a single dose per subject.
4.7 Prediction correction
Prediction-corrected Visual Predictive Checks (pcVPCs) were developed to account for variability in independent variables (e.g., dose, covariates, or time-varying predictors) that may influence predictions but cannot easily be stratified. In such cases, raw observed and simulated values may differ simply due to differences in design rather than model misspecification. The pcVPC corrects these design-related differences by normalizing both observed and simulated outcomes relative to the population-level prediction at each time point (typically using the model-predicted median). This allows for a clearer comparison of the model’s predictive performance across heterogeneous conditions.
To perform a prediction corrected VPC in Pumas, simply set the prediction_correction
keyword to true.
= vpc(foce_fit; prediction_correction = true)
my_pcvpc vpc_plot(my_pcvpc; figurelegend = (; position = :t, nbanks = 2, tellheight = true))
[ Info: Continuous VPC
In this instance, we see that the VPC is maybe leaning a bit more towards model misspecification if we think prediction correction is appropriate. The model bands seem to clear (got to zero) faster than the observed data.
5 Conclusion
In this tutorial, we saw how to produce a VPC from a fitted model. We showed how to look for misspecification and some of the levers and controls a user has to customize the plot.
6 Further Reading
If you want to read more on the application and properties of visual predictive checks, please have a look at the following papers.
1. Karlsson, M. O., & Holford, N. H. (2008). A tutorial on visual predictive checks. PAGE 17, Abstract 1434. http://www.page-meeting.org/?abstract=1434
2. Bergstrand, M., Hooker, A. C., Wallin, J. E., & Karlsson, M. O. (2011). Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models. CPT: Pharmacometrics & Systems Pharmacology, 89(2), 204–217. https://doi.org/10.1038/clpt.2010.268
3. Ette, E. I., & Williams, P. J. (2004). Population pharmacokinetics II: estimation methods. Annals of Pharmacotherapy, 38(11), 1907–1915. https://doi.org/10.1345/aph.1E225
4. Bonate, P. L. (2011). Pharmacokinetic-Pharmacodynamic Modeling and Simulation (2nd ed.). Springer. Chapters on model evaluation and simulation diagnostics.