Visual Predictive Checks for Continuous Data Models

Author

Patrick Kofod Mogensen 

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

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.

warfarin_pk_model = @model begin

    @metadata begin
        desc = "Warfarin 1-compartment PK model (PD removed)"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Clearance (L/hr)
        """
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        """
        Central volume (L)
        """
        pop_Vc ∈ RealDomain(lower = 0.0, init = 8.11)
        """
        Absorption lag time (hr)
        """
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)

        # Inter-individual variability
        """
          - ΩCL: Clearance
          - ΩVc: Central volume
        """
        pk_Ω ∈ PDiagDomain([0.01, 0.01])
        # Residual variability
        """
        σ_prop: Proportional error
        """
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)    # mean = 0, covariance = pk_Ω
    end

    @covariates begin
        """
        FSZCL: Clearance scaling factor
        """
        FSZCL
        """
        FSZV: Volume scaling factor
        """
        FSZV
    end

    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])

        tabs = pop_tabs
        Ka = log(2) / tabs
    end

    @vars begin
        cp := Central / Vc
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end

    @derived begin
        """
        Concentration (ng/mL)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2))
    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
warfarin_data_wide = @chain dataset("pumas/warfarin_pumas") begin
    @rtransform begin
        # Scaling factors
        :FSZV = :wtbl / 70            # volume scaling
        :FSZCL = (:wtbl / 70)^0.75     # clearance scaling (allometric)
    end
end
330×12 DataFrame
305 rows omitted
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 
pop = read_pumas(
    warfarin_data_wide;
    id = :id,
    time = :time,
    amt = :amt,
    cmt = :cmt,
    evid = :evid,
    covariates = [:sex, :wtbl, :FSZV, :FSZCL],
    observations = [:conc],
)
Population
  Subjects: 32
  Covariates: sex, wtbl, FSZV, FSZCL
  Observations: conc
foce_fit = fit(warfarin_pk_model, pop, init_params(warfarin_pk_model), FOCE())
[ 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 plots
  • vpc_plot(...) from PumasUtilities that takes the data and presents it in a plot
VPC with 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:

myvpc = vpc(foce_fit)
[ 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(
    fpm::AbstractFittedPumasModel;
    samples::Integer = 499,
    observations::Vector{Symbol} = [keys(fpm.data[1].observations)[1]],
    covariates::Vector{Symbol} = [:time],
    stratify_by::Vector{Symbol} = Symbol[],
    bandwidth = nothing,
    prediction_correction::Bool = false,
)

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 model
  • observations is a vector of symbols (e.g. [:conc]) to control what endpoint to perform the VPC for if there are multiple present in the model
  • covariates (default is :time) is a vector of symbols (e.g. [:tad]) to control the first axis variable
  • stratify_by to control stratification of VPC calculations and plotting
  • bandwidth is a number (default is 2.0 for continuous VPCs) to control the level of smoothing
  • prediction_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:

foce_fit_CL = fit(
    warfarin_pk_model,
    pop,
    (init_params(warfarin_pk_model)..., pop_CL = 0.15),
    FOCE();
    constantcoef = (:pop_CL,),
    optim_options = (; show_trace = false),
)
[ 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

myvpc_CL = vpc(foce_fit_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.

my_vpc_sex = vpc(foce_fit; stratify_by = [: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.

my_vpc_sex = vpc(foce_fit; stratify_by = [:sex], bandwidth = 5)
[ 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 of 2

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.

my_vpc_sex = vpc(foce_fit; covariates = [:tad], stratify_by = [:dosenum])
[ 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.

my_pcvpc = vpc(foce_fit; prediction_correction = true)
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.