MI-210 Essentials of Population PKPD M&S — Pumas Edition
Model qualification (or evaluation) methods help assess how well a model describes the observed data and whether it is suitable for its intended purpose:
Standard goodness-of-fit diagnostics each target a specific aspect of model adequacy:
If the model is correctly specified, CWRES should be approximately \(N(0,1)\) with no systematic trends. Deviations from this expectation guide targeted model refinement.
Asymptotic standard errors derived from the Fisher information matrix assume that the likelihood surface is approximately quadratic (i.e., the parameter estimates are normally distributed) near the minimum. This assumption can fail for small datasets, when variance parameters are estimated near their boundary (e.g., \(\omega^2 \approx 0\)), or when the likelihood surface is flat or irregular. In such cases, asymptotic confidence intervals may be misleading. Bootstrap resampling makes no distributional assumptions: by repeatedly resampling subjects with replacement and re-estimating the model, it constructs empirical confidence intervals that reflect the actual sampling distribution of the estimator. However, bootstrap is not without limitations — if the convergence rate across resamples falls below approximately 80%, this signals model instability or overparameterization, and the bootstrap results themselves become unreliable.
Bootstrap provides non-parametric confidence intervals for parameter estimates by resampling subjects with replacement.
The book uses a 1-compartment IV model with combined error (ADVAN1 TRANS2):
pk_df = CSV.read(joinpath(data_dir, "prob_5_1.csv"), DataFrame)
rename!(pk_df,
:ID => :id, :DV => :dv, :AMT => :amt, :TIME => :time,
:II => :ii, :ADDL => :addl, :RATE => :rate,
)
# Prepare for Pumas
pk_df.evid = ifelse.(pk_df.amt .> 0, 1, 0)
pk_df.cmt = fill(1, nrow(pk_df))
pk_df.dv = ifelse.(pk_df.evid .== 1, missing, pk_df.dv)
pop = read_pumas(pk_df; covariates = [:HT, :WT, :CLCR, :SEX, :AGE])
println("Subjects: ", length(pop))Subjects: 40# 1-cpt IV model with BLOCK(2) OMEGA and combined error
# From 501.ctl: ADVAN1 TRANS2, Y=F + F*ERR(1) + ERR(2)
pk_model = @model begin
@metadata begin
desc = "1-Cpt IV PopPK (Bootstrap Base Model)"
end
@param begin
tvCL ∈ RealDomain(lower = 0.0, init = 13.0)
tvV ∈ RealDomain(lower = 0.0, init = 75.0)
Ω ∈ PSDDomain(2)
σ_prop ∈ RealDomain(lower = 0.0, init = 0.2)
σ_add ∈ RealDomain(lower = 0.0, init = 1.0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates WT CLCR SEX AGE HT
@pre begin
CL = tvCL * exp(η[1])
Vc = tvV * exp(η[2])
end
@dynamics Central1
@derived begin
cp := @. Central / Vc
dv ~ @. Normal(cp, sqrt((cp * σ_prop)^2 + σ_add^2))
end
end┌ Warning: Covariates WT, CLCR, SEX, AGE and HT are not used in the model. └ @ Pumas ~/.julia/packages/Pumas/GZeMg/src/dsl/model_macro.jl:3399
PumasModel
Parameters: tvCL, tvV, Ω, σ_prop, σ_add
Random effects: η
Covariates: WT, CLCR, SEX, AGE, HT
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: dv
Observed: dvinit_params = (
tvCL = 13.0, tvV = 75.0,
Ω = [0.04 0.02; 0.02 0.04],
σ_prop = 0.2, σ_add = 1.0,
)
base_fit = fit(pk_model, pop, init_params, Pumas.FOCEI())
base_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.140471e+02 1.997472e+02 * time: 0.01035308837890625 1 8.501166e+02 1.437327e+02 * time: 0.9784431457519531 2 8.413936e+02 3.418607e+02 * time: 0.9819841384887695 3 7.401033e+02 1.272967e+02 * time: 0.9856510162353516 4 7.248790e+02 4.166194e+01 * time: 1.3897311687469482 5 7.178590e+02 2.541013e+01 * time: 1.3927490711212158 6 7.106983e+02 7.097696e+01 * time: 1.3959240913391113 7 6.951798e+02 4.877067e+01 * time: 1.399273157119751 8 6.881324e+02 3.937486e+01 * time: 1.402212142944336 9 6.810644e+02 1.182857e+01 * time: 1.4055390357971191 10 6.795778e+02 6.428324e+00 * time: 1.4088010787963867 11 6.792989e+02 2.558174e+00 * time: 1.41218900680542 12 6.792423e+02 1.421276e+00 * time: 1.4154059886932373 13 6.792198e+02 1.683910e-01 * time: 1.4188671112060547 14 6.792196e+02 7.468799e-02 * time: 1.4218711853027344 15 6.792196e+02 3.711636e-02 * time: 1.4250400066375732 16 6.792196e+02 5.935271e-03 * time: 1.4278860092163086 17 6.792196e+02 1.131421e-03 * time: 1.430298089981079 18 6.792196e+02 1.236725e-04 * time: 1.4328479766845703
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 40
Observation records: Active Missing
dv: 560 0
Total: 560 0
Number of parameters: Constant Optimized
0 7
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -679.2196
--------------------
Estimate
--------------------
tvCL 11.806
tvV 84.841
Ω₁,₁ 0.05785
Ω₂,₁ 0.011617
Ω₂,₂ 0.061883
σ_prop 0.20242
σ_add 0.0028996
--------------------base_infer = infer(base_fit)
base_infer[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: Closed form
Number of subjects: 40
Observation records: Active Missing
dv: 560 0
Total: 560 0
Number of parameters: Constant Optimized
0 7
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -679.2196
-------------------------------------------------------------
Estimate SE 95.0% C.I.
-------------------------------------------------------------
tvCL 11.806 0.46518 [ 10.895 ; 12.718 ]
tvV 84.841 3.5631 [ 77.858 ; 91.825 ]
Ω₁,₁ 0.05785 0.01398 [ 0.030449 ; 0.08525 ]
Ω₂,₁ 0.011617 0.0088497 [ -0.0057282; 0.028962 ]
Ω₂,₂ 0.061883 0.013041 [ 0.036322 ; 0.087443 ]
σ_prop 0.20242 0.0074533 [ 0.18781 ; 0.21703 ]
σ_add 0.0028996 0.00059726 [ 0.001729 ; 0.0040702]
-------------------------------------------------------------In Pumas, bootstrap is performed via infer(fit, Bootstrap()):
# Bootstrap with 200 replicates (use more for production — e.g., 1000)
boot_result = infer(base_fit, Pumas.Bootstrap(samples = 200))
boot_result[ Info: Bootstrap inference finished.
Bootstrap inference results
Dynamical system type: Closed form
Number of subjects: 40
Observation records: Active Missing
dv: 560 0
Total: 560 0
Number of parameters: Constant Optimized
0 7
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -679.2196
------------------------------------------------------------
Estimate SE 95.0% C.I.
------------------------------------------------------------
tvCL 11.806 0.48494 [ 10.889 ; 12.684 ]
tvV 84.841 3.5461 [ 78.503 ; 92.273 ]
Ω₁,₁ 0.05785 0.013835 [ 0.032432 ; 0.0842 ]
Ω₂,₁ 0.011617 0.0090141 [ -0.0059066; 0.026267 ]
Ω₂,₂ 0.061883 0.012881 [ 0.035084 ; 0.082688 ]
σ_prop 0.20242 0.0072407 [ 0.1903 ; 0.21753 ]
σ_add 0.0028996 0.0010607 [ 1.1109e-5; 0.0049531]
------------------------------------------------------------
Unique resampled populations: 200 out of 200
No stratification.In NONMEM, bootstrap requires external tools (PsN, metrumrg) to resample datasets and re-run the model hundreds of times. In Pumas, infer(fit, Bootstrap()) handles everything in a single call.
VPC compares the distribution of observed data to simulated data from the model. It assesses whether the model can reproduce the key features of the observed data.
In NONMEM, VPC requires $SIMULATION and post-processing with R/PsN. In Pumas, use vpc() and vpc_plot():
vpc_result = vpc(base_fit, 200) # 200 simulations[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 200
Subjects in data: 40
Stratification variable(s): None
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])vpc_plot(vpc_result)
Log-likelihood profiling provides confidence intervals for individual parameters by fixing one parameter at a range of values and re-estimating the rest:
# Profile tvCL over a range around the estimate
# (computationally intensive — commented for rendering speed)
# cl_range = range(coef(base_fit).tvCL * 0.8, coef(base_fit).tvCL * 1.2, length = 20)
# profiles = [fit(pk_model, pop, merge(coef(base_fit), (tvCL = cl,)), FOCEI();
# constantcoef = (tvCL = cl,)) for cl in cl_range]
# ll_values = loglikelihood.(profiles)insp = inspect(base_fit)
insp_df = DataFrame(insp)
obs = @rsubset(insp_df, !ismissing(:dv))
fig = Figure(size = (800, 600))
ax1 = Axis(fig[1, 1], xlabel = "PRED", ylabel = "DV", title = "DV vs PRED")
scatter!(ax1, obs.dv_pred, obs.dv, color = :navy, markersize = 4, alpha = 0.5)
ablines!(ax1, 0, 1, color = :black)
ax2 = Axis(fig[1, 2], xlabel = "IPRED", ylabel = "DV", title = "DV vs IPRED")
scatter!(ax2, obs.dv_ipred, obs.dv, color = :navy, markersize = 4, alpha = 0.5)
ablines!(ax2, 0, 1, color = :black)
ax3 = Axis(fig[2, 1], xlabel = "PRED", ylabel = "CWRES", title = "CWRES vs PRED")
scatter!(ax3, obs.dv_pred, obs.dv_wres, color = :navy, markersize = 4, alpha = 0.5)
hlines!(ax3, 0, color = :black)
ax4 = Axis(fig[2, 2], xlabel = "Time (hr)", ylabel = "CWRES", title = "CWRES vs TIME")
scatter!(ax4, obs.time, obs.dv_wres, color = :navy, markersize = 4, alpha = 0.5)
hlines!(ax4, 0, color = :black)
fig[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
Before evaluating a model’s adequacy, the right model structure must be selected. Exploratory data analysis provides the first line of evidence:
A concentration-effect (CE) plot reveals whether the response is in equilibrium with plasma concentrations:
| CE Plot Pattern | Interpretation | Model Implication |
|---|---|---|
| No loop (single curve) | Equilibrium — effect tracks concentration | Direct response (Emax, sigmoid Emax) |
| Counter-clockwise loop | Effect lags behind concentration | Effect compartment (ke0 delay) |
| Clockwise loop | Tolerance, sensitization, or indirect mechanism | Indirect response (IDR Types I–IV) or tolerance model |
When comparing more than two candidate models:
# Compare multiple models systematically
fits = [fit1, fit2, fit3]
for (i, f) in enumerate(fits)
println("Model $i — AIC: ", round(aic(f), digits=1),
" BIC: ", round(bic(f), digits=1),
" -2LL: ", round(-2*loglikelihood(f), digits=1))
end
# Or use metrics_table for side-by-side comparison
metrics_table(fit1, fit2, fit3)
# LRT for nested pairs only
lrtest(fit1, fit2) # only valid when fit1 is nested within fit2| Topic | Tutorial | Key Additions |
|---|---|---|
| Lecture: Graphical Diagnostics | PKPD W6 | DV vs PRED/IPRED interpretation, CWRES patterns, time-dependent misspecification |
| Lecture: Residual Diagnostics & Error Models | PKPD W7 | Additive vs proportional vs combined error, residual analysis methodology |
| Lecture: Base Model & Likelihood Comparison | PKPD W5 | -2LL, LRT, AIC/BIC theory, model comparison workflow |
| PD model selection framework | PD Model Selection | EDA → mechanism → statistics workflow, hysteresis detection, warfarin case study |
| Bayesian model comparison | Bayesian Model Comparison | WAIC, LOO-CV, Bayes factors for non-nested models |
| VPC customization | Continuous VPC | Stratified VPC, prediction-corrected VPC, log-scale |