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.053817033767700195
     1     1.313047e+04     2.609679e+04
 * time: 2.3301382064819336
     2     9.658030e+03     1.909689e+04
 * time: 2.3730661869049072
     3     4.233250e+03     8.101602e+03
 * time: 2.416604995727539
     4     2.367624e+03     4.268593e+03
 * time: 2.460447072982788
     5     1.313943e+03     2.054844e+03
 * time: 2.503981113433838
     6     8.374399e+02     1.010773e+03
 * time: 2.549386978149414
     7     6.115489e+02     4.863561e+02
 * time: 2.5960099697113037
     8     5.146179e+02     2.031565e+02
 * time: 2.641162157058716
     9     4.860731e+02     7.536873e+01
 * time: 2.684295177459717
    10     4.801210e+02     5.238009e+01
 * time: 2.7266781330108643
    11     4.795642e+02     5.405505e+01
 * time: 2.768522024154663
    12     4.795291e+02     5.414517e+01
 * time: 2.80948805809021
    13     4.794853e+02     5.383276e+01
 * time: 2.8515281677246094
    14     4.793492e+02     5.211648e+01
 * time: 3.1045820713043213
    15     4.790578e+02     4.736685e+01
 * time: 3.1456902027130127
    16     4.784683e+02     3.589484e+01
 * time: 3.1865429878234863
    17     4.776996e+02     1.930725e+01
 * time: 3.2278270721435547
    18     4.772243e+02     1.316223e+01
 * time: 3.2704110145568848
    19     4.770920e+02     1.346641e+01
 * time: 3.312412977218628
    20     4.770711e+02     1.365789e+01
 * time: 3.3549320697784424
    21     4.770686e+02     1.368452e+01
 * time: 3.391733169555664
    22     4.770612e+02     1.372786e+01
 * time: 3.4304020404815674
    23     4.770441e+02     1.378138e+01
 * time: 3.4701991081237793
    24     4.770002e+02     1.384496e+01
 * time: 3.5111000537872314
    25     4.769031e+02     1.387480e+01
 * time: 3.5530591011047363
    26     4.767248e+02     2.003031e+01
 * time: 3.5961031913757324
    27     4.765128e+02     3.331465e+01
 * time: 3.637253999710083
    28     4.763904e+02     3.753509e+01
 * time: 3.6781251430511475
    29     4.763566e+02     3.512805e+01
 * time: 3.7175650596618652
    30     4.763399e+02     3.256639e+01
 * time: 3.75640606880188
    31     4.763051e+02     2.857696e+01
 * time: 3.7967259883880615
    32     4.762186e+02     2.633532e+01
 * time: 3.8376049995422363
    33     4.759910e+02     2.972803e+01
 * time: 3.878458023071289
    34     4.754039e+02     3.513095e+01
 * time: 3.923084020614624
    35     4.739900e+02     4.343181e+01
 * time: 3.9696691036224365
    36     4.718772e+02     5.172641e+01
 * time: 4.014244079589844
    37     4.699605e+02     5.109410e+01
 * time: 4.056824207305908
    38     4.679301e+02     1.837720e+01
 * time: 4.1008219718933105
    39     4.677728e+02     9.841366e+00
 * time: 4.144431114196777
    40     4.677644e+02     9.868356e+00
 * time: 4.185774087905884
    41     4.677619e+02     9.877285e+00
 * time: 4.223919153213501
    42     4.677607e+02     9.881556e+00
 * time: 4.2623021602630615
    43     4.677536e+02     9.897820e+00
 * time: 4.301727056503296
    44     4.677390e+02     9.912120e+00
 * time: 4.419183969497681
    45     4.676965e+02     9.918954e+00
 * time: 4.462076187133789
    46     4.675898e+02     1.060357e+01
 * time: 4.502827167510986
    47     4.673112e+02     1.524992e+01
 * time: 4.542540073394775
    48     4.666264e+02     2.323583e+01
 * time: 4.582597017288208
    49     4.649869e+02     2.870754e+01
 * time: 4.624516010284424
    50     4.638128e+02     3.688505e+01
 * time: 4.664412021636963
    51     4.617598e+02     3.172279e+01
 * time: 4.704107999801636
    52     4.609261e+02     2.389881e+01
 * time: 4.743338108062744
    53     4.601257e+02     7.743696e+00
 * time: 4.781212091445923
    54     4.600344e+02     2.839021e+00
 * time: 4.818940162658691
    55     4.600202e+02     6.791380e-01
 * time: 4.855629205703735
    56     4.600185e+02     1.392758e-01
 * time: 4.890884160995483
    57     4.600183e+02     3.003773e-02
 * time: 4.924974203109741
    58     4.600183e+02     4.113068e-03
 * time: 4.95789098739624
    59     4.600183e+02     3.717389e-03
 * time: 4.992350101470947
    60     4.600183e+02     3.347682e-03
 * time: 5.02669620513916
    61     4.600183e+02     3.347291e-03
 * time: 5.066483974456787
    62     4.600183e+02     3.347291e-03
 * time: 5.117146015167236
    63     4.600183e+02     2.864074e-03
 * time: 5.158914089202881
    64     4.600183e+02     2.864074e-03
 * time: 5.206866979598999
    65     4.600183e+02     2.864074e-03
 * time: 5.258052110671997
    66     4.600183e+02     2.864074e-03
 * time: 5.312542200088501
    67     4.600183e+02     2.864074e-03
 * time: 5.392261028289795
    68     4.600183e+02     2.864074e-03
 * time: 5.4431211948394775
    69     4.600183e+02     2.864074e-03
 * time: 5.497143983840942
    70     4.600183e+02     2.864074e-03
 * time: 5.548158168792725

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.