Visual Predictive Check (VPC)

Pumas allows you to generate VPC quantiles for your data and model simulations and utilize the Julia plotting capabilities to generate the relevant VPC plots. This is allowed with the vpc function discussed below that returns a VPC object, this can then be passed to the vpc_plot function to generate the plots.

Overview

vpc(fpm::AbstractFittedPumasModel;
  samples::Integer=499,
  qreg_method=IP(),
  observations::Array{Symbol} = [keys(extract_subject(fpm.data[1]).observations)[1]],
  stratify_by::Array{Symbol} = Symbol[],
  quantiles::NTuple{3,Float64}=(0.1, 0.5, 0.9),
  level::Real=0.95,
  ensemblealg::DiffEqBase.EnsembleAlgorithm=EnsembleSerial(),
  bandwidth::Real=2,
  maxnumstrats::Array{<:Integer} = [6 for i in 1:length(stratify_by)],
  covariates::Array{Symbol} = [:time],
  smooth::Bool = true,
  rng::AbstractRNG=default_rng(),
  obstimes::AbstractVector = [])

The vpc computes the quantiles for VPC for an AbstractFittedPumasModel with simulated confidence intervals around the empirical quantiles based on reps simulated populations.

The number of repetitions of the simulations are passed in with the reps positional argument.

The following keyword arguments are supported:

  • samples: The number of simulated populations to generate, defaults to 499.
  • quantiles: A three-tuple of the quantiles for which the quantiles will be computed. The default is (0.1, 0.5, 0.9) which computes the 10th, 50th and 90th percentile.
  • level: Probability level to use for the simulated confidence intervals around each of the quantiles. The default is 0.95.
  • observations: The name of the dependent variable to use for the VPCs, passed in as a Symbol wrapped in an array. The default is the first dependent variable in the population. Even though the argument takes an array note that only single derived variable at a time is supported.
  • stratify_by: The covariates to be used for stratification. Takes an array of the Symbols of the stratification covariates.
  • ensemblealg: This is passed to the simobs call while the reps simulations. For more description check the docs for simobs.
  • bandwidth: The kernel bandwidth in the quantile regression. If you are seeing NaNs or an error, increasing the bandwidth should help in most cases. With higher values of the bandwidth you will get more smoothened plots of the quantiles so it's a good idea to check with your data the right bandwidth.
  • maxnumstrats: The maximum number of strata for each covariate. Takes an array with the number of strata for the corresponding covariate, passed in stratify_by. It defaults to 6 for each of the covariates.
  • covariates: The independent variable for VPC, defaults to [:time]. Similar to observations it requires a single element array.
  • obstimes: The times for simulation in case of continuous VPC, same as obstimes in simobs. Defaults to union of all subject's unique times in the data where observation is not missing.
  • qreg_method: Defaults to IP(). For most users the method used in quantile regression is not going to be of concern, but if you see large run times switching qreg_method to IP(true) should help in improving the performance with a tradeoff in the accuracy of the fitting.

Hence a vpc call with 100 repetitions of simulation would be:

res = fit(model, data, param,Pumas.FOCE())
vpc_fpm = vpc(res, samples = 100)

Stratification on wt covariate

vpc_fpm_stratwt = vpc(res, stratify_by=[:wt])
# run only 100 simulations 
vpc_fpm_stratispm = vpc(res, samples = 100, stratify_by=[:isPM])

The VPC object obtained as the result contains the following fields:

  • simulated_quantiles::DataFrame: DataFrame of the simulated quantiles result.
  • popvpc::PopVPC: A PopVPC object with the observation quantiles (data_quantiles), data (data), stratification covariate if specified (stratify_by) and the dependent variable (dv).
  • level::Float64: The simulation CI's level.

Since the quantiles are stored as DataFrames it is very easy to use the result of vpc and also export it using CSV.jl to your disk.

vpc_fpm.simulated_quantiles
vpc_fpm.popvpc.data_quantiles

using CSV
CSV.write("sim_quantiles.csv", vpc_fpm.simulated_quantiles)

In Pumas, for models with continuous derived variables, the non-parametric quantile regression approach discussed in Jamsen et al. is used. Discrete and Time to Event models are also supported through the vpc function and can be used without any change in syntax.

Full examples:

Let's take a look at a complete example, first we define the model and generate round-trip data.

julia> model = @model begin
         @param begin
           tvcl ∈ RealDomain(lower=0)
           tvv ∈ RealDomain(lower=0)
           pmoncl ∈ RealDomain(lower = -0.99)
           Ω ∈ PDiagDomain(2)
           σ_prop ∈ RealDomain(lower=0)
         end
       
         @random begin
           η ~ MvNormal(Ω)
         end
       
         @covariates wt isPM
       
         @pre begin
           CL = tvcl * (1 + pmoncl*isPM) * (wt/70)^0.75 * exp(η[1])
           Vc    = tvv * (wt/70) * exp(η[2])
         end
       
         @dynamics Central1
       
         @derived begin
           cp = @. 1000*(Central / Vc)
           dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
         end
       endPumasModel
  Parameters: tvcl, tvv, pmoncl, Ω, σ_prop
  Random effects: η
  Covariates: wt, isPM
  Dynamical variables: Central
  Derived: cp, dv
  Observed: cp, dv
julia> ev = DosageRegimen(100, time=0, addl=2, ii=24)DosageRegimen Row │ time cmt amt evid ii addl rate duration ss route │ Float64 Int64 Float64 Int8 Float64 Int64 Float64 Float64 Int8 NCA.Route ─────┼─────────────────────────────────────────────────────────────────────────────────── 1 │ 0.0 1 100.0 1 24.0 2 0.0 0.0 0 NullRoute
julia> s1 = Subject(id=1, events=ev, covariates=(isPM=1, wt=70))Subject ID: 1 Events: 3 Covariates: isPM, wt
julia> param = ( tvcl = 4.0, tvv = 70, pmoncl = -0.7, Ω = Diagonal([0.09,0.09]), σ_prop = 0.04 )(tvcl = 4.0, tvv = 70, pmoncl = -0.7, Ω = [0.09 0.0; 0.0 0.09], σ_prop = 0.04,)
julia> choose_covariates() = (isPM = rand([1, 0]), wt = rand(55:80))choose_covariates (generic function with 1 method)
julia> pop_with_covariates = Population(map(i -> Subject(id=i, events=ev, covariates=choose_covariates()),1:10))Population Subjects: 10 Covariates: isPM, wt Observations:
julia> obs = simobs(model, pop_with_covariates, param, obstimes=0:1:60)Simulated population (Vector{<:Subject}) Simulated subjects: 10 Simulated variables: cp, dv
julia> simdf = DataFrame(obs)640×10 DataFrame Row │ id time cp dv evid amt cmt rate isPM wt │ String Float64 Float64? Float64? Int8 Float64? Symbol? Float64? Int64? Int64? ─────┼─────────────────────────────────────────────────────────────────────────────────────────────── 1 │ 1 0.0 missing missing 1 100.0 Central 0.0 1 77 2 │ 1 0.0 1132.47 1570.88 0 missing missing 0.0 1 77 3 │ 1 1.0 1118.51 1134.38 0 missing missing 0.0 1 77 4 │ 1 2.0 1104.71 1316.7 0 missing missing 0.0 1 77 5 │ 1 3.0 1091.09 1189.15 0 missing missing 0.0 1 77 6 │ 1 4.0 1077.64 1042.32 0 missing missing 0.0 1 77 7 │ 1 5.0 1064.35 1292.43 0 missing missing 0.0 1 77 8 │ 1 6.0 1051.22 875.673 0 missing missing 0.0 1 77 9 │ 1 7.0 1038.26 1167.83 0 missing missing 0.0 1 77 10 │ 1 8.0 1025.46 986.493 0 missing missing 0.0 1 77 11 │ 1 9.0 1012.81 832.218 0 missing missing 0.0 1 77 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 631 │ 10 51.0 2362.4 2133.28 0 missing missing 0.0 1 74 632 │ 10 52.0 2287.82 1918.42 0 missing missing 0.0 1 74 633 │ 10 53.0 2215.6 2218.34 0 missing missing 0.0 1 74 634 │ 10 54.0 2145.65 2242.76 0 missing missing 0.0 1 74 635 │ 10 55.0 2077.91 2199.11 0 missing missing 0.0 1 74 636 │ 10 56.0 2012.31 2039.51 0 missing missing 0.0 1 74 637 │ 10 57.0 1948.78 1610.14 0 missing missing 0.0 1 74 638 │ 10 58.0 1887.26 1395.52 0 missing missing 0.0 1 74 639 │ 10 59.0 1827.68 1689.23 0 missing missing 0.0 1 74 640 │ 10 60.0 1769.98 1790.8 0 missing missing 0.0 1 74 619 rows omitted
julia> data = read_pumas(simdf, time=:time, covariates=[:isPM, :wt])Population Subjects: 10 Covariates: isPM, wt Observations: dv
julia> res = fit(model, data, param, Pumas.FOCE())Iter Function value Gradient norm 0 4.335786e+03 2.271171e+01 * time: 0.033969879150390625 1 4.334556e+03 1.197261e+01 * time: 1.4230308532714844 2 4.333938e+03 5.276184e+00 * time: 1.43110990524292 3 4.332686e+03 5.925609e+00 * time: 1.4363529682159424 4 4.332413e+03 5.449928e+00 * time: 1.44455885887146 5 4.332202e+03 9.754982e-01 * time: 1.449753999710083 6 4.332059e+03 2.036235e+00 * time: 1.4547648429870605 7 4.331952e+03 2.479859e+00 * time: 1.4599719047546387 8 4.331927e+03 8.509045e-01 * time: 1.465043067932129 9 4.331922e+03 1.404262e-01 * time: 1.4699969291687012 10 4.331920e+03 4.233475e-01 * time: 1.474808931350708 11 4.331918e+03 4.145562e-01 * time: 1.4796299934387207 12 4.331917e+03 1.369603e-01 * time: 1.484264850616455 13 4.331917e+03 1.745298e-02 * time: 1.488813877105713 14 4.331917e+03 7.403672e-03 * time: 1.4932708740234375 15 4.331917e+03 1.438446e-04 * time: 1.4976270198822021 FittedPumasModel Successful minimization: true Likelihood approximation: Pumas.FOCE Log-likelihood value: -4331.9173 Number of subjects: 10 Number of parameters: Fixed Optimized 0 6 Observation records: Active Missing dv: 610 0 Total: 610 0 --------------------- Estimate --------------------- tvcl 5.865 tvv 77.285 pmoncl -0.80491 Ω₁,₁ 0.061998 Ω₂,₂ 0.054158 σ_prop 0.043397 ---------------------
julia> vpc_fpm = vpc(res);

julia> vpc_plot(vpc_fpm)

FPM VPC

Alternatively, using a PumasEMModel

emmodel = @emmodel begin
  @random begin
    CLdw ~ 1 + isPM | LogNormal # CLdw = μ_CL * exp(μ_CL_isPM * isPM) * exp(η_CL)
    Vcdw ~ 1        | LogNormal # Vcdw = μ_Vc * exp(η_Vc)
  end
  @covariates wt
  @pre begin
    CL = CLdw * (wt/70)^0.75
    Vc = Vcdw * (wt/70)
  end
  @dynamics Central1
  @post begin
    cp = 1000 * Central / Vc
  end
  @error begin
    dv ~ ProportionalNormal(cp)
  end
end

emparam = (
    CLdw = (4.0, -1.2), # -1.2 ≈ log(1 - 0.7) = log(1 - pmoncl)
    Vcdw    = 70,
    Ω = (0.09,0.09),
    σ = (0.04,)
)

emres = fit(emmodel, data, emparam, Pumas.LaplaceI())

vpc_emfpm = vpc(emres)

Note that because isPM equals either 0 or 1, pmoncl = exp(μ_CL_isPM) - 1, allowing us to perform a change of variables and express the earlier PumasModel as a PumasEMModel.

Tip

In PumasModels and PumasEMModels, you can increase performance by preprocessing covariates. Instead of wt, you can try storing (wt/70) and (wt/70)^0.75 as appropriately named covariates.

vpc_plot(vpc_emfpm)

EMMODEL-FPM VPC

Let's generate a VPC stratified on isPM covariate with stratify_by. Note that multiple covariates can be passed to stratify_by to get cross stratified VPC plots, in that case the maxnumstrats will have to be adjusted accordingly.

resem_saem = fit(emmodel, data, (CLdw = (5.0,0.0), Vcdw = 50.0), Pumas.SAEM())

vpc_fpmem_saem = vpc(resem_saem, samples = 100)
vpc_plot(vpc_fpmem_saem)

SAEM-FPM VPC

Tip

When fitting with SAEM, it is strongly recomended to use values for the Ω parameters much larger than your estiamtes of the true values, or to avoid providing them so they default to 1.0 along the diagonal. Using small values such as emparams will result in convergence failure.

vpc_fpm = vpc(res, stratify_by = [:isPM])
vpc_plot(vpc_fpm)

vpc_fpm_strat

vpc_saem_fpm = vpc(resem_saem, stratify_by = [:isPM]);
vpc_plot(vpc_saem_fpm)

vpc_saemfpm_strat

Let's pass in the time grid with the obstimes keyword argument to vpc for the simulation time points.

vpc_fpm = vpc(res, 100, obstimes = 0.0:5:60.0)
vpc_plot(vpc_fpm)

vpc_fpm_obstimes

By specifying the plot and the path the VPC plot can be saved to the disk as below:

p = vpc_plot(vpc_fpm)
CairoMakie.save("obstimesvpc.png", p)

More examples of VPC plots can be found in the plotting page's vpc_plot section.

Tip

If your vpc run gives you an error message similar to ERROR: LinearAlgebra.LAPACKException(2), increasing the bandwidth should fix it in most cases.

Note

For most users, the method used in quantile regression is not going to be of concern, but if you see large run times switching qreg_method to IP(true) after loading QuantileRegressions.jl should help in improving the performance with a tradeoff in the accuracy of the regression fitting.

Prediction corrected VPCs (pcVPCs) are work in progress and will become available soon.