Back to Claude Scientific Skills

Population PK: estimation, covariates, BLQ, and diagnostics

skills/pkpd-modeling/references/population-pk.md

2.61.08.0 KB
Original Source

Population PK: estimation, covariates, BLQ, and diagnostics

The model has three layers, and they are diagnosed separately

Structural:   C_ij = f(theta_i, dose_i, t_ij)
Individual:   theta_i = theta_pop * exp(eta_i),   eta_i ~ N(0, Omega)
Residual:     y_ij = C_ij * (1 + eps_prop) + eps_add,   eps ~ N(0, Sigma)

Almost every difficult population model problem is a layer confusion: an extra compartment added to absorb between-occasion variability, a covariate added to fix a misspecified absorption model, a proportional error inflated to cover a structural bias at low concentrations. Identify the layer before changing anything.

Exponential (log-normal) inter-individual variability is the default because clearances and volumes are positive and right-skewed. omega is reported as an approximate CV: CV ≈ sqrt(exp(omega²)−1), which is close to omega itself below about 30%.

Estimation methods

MethodCharacterWhen
FOFirst-order linearisation about eta = 0Obsolete for final models; biased with large IIV. Still useful for initial estimates
FOCELinearises about the individual etaThe workhorse
FOCE-I (INTERACTION)FOCE with eta-epsilon interactionRequired whenever the residual error is proportional or combined. Omitting INTERACTION with a proportional error model is a common and consequential mistake
LaplaceSecond-order expansionNeeded for non-continuous data (categorical, count, time-to-event) and for LIKELIHOOD/-2LL models
SAEMStochastic approximation EMRobust to poor initial estimates and to complex models; less prone to local minima; does not itself give an objective function for comparison — follow with an IMP evaluation
Importance sampling (IMP)Monte Carlo integrationAccurate objective function; usually run after SAEM
MCMC / NUTSFull BayesianNONMEM 7.6 adds NUTS; also Stan, Torsten, nlmixr2

Compare objective functions only between models fitted with the same method on the same records. An OFV drop obtained by switching from FOCE to IMP is not evidence about the model.

Below the limit of quantification

Beal's methods, and what they actually do:

MethodTreatmentVerdict
M1Discard all BLQ observationsBiased upward whenever a meaningful fraction is BLQ; acceptable only when that fraction is small
M2Discard BLQ, condition the likelihood on being above LLOQBetter than M1, rarely used
M3Maximum likelihood: BLQ contributes the probability that the observation is below LLOQThe reference method. Needs F_FLAG=1 and the Laplace method
M4M3 conditioned on the concentration being positiveMarginal improvement over M3
M5Substitute LLOQ/2Biased; still common; at least state it
M6Substitute LLOQ/2 for the first BLQ in a run, discard the restAd hoc
M7Substitute 0Worst; badly biases the terminal phase

Rule of thumb: below roughly 10% BLQ, M1 is defensible; above that, use M3. The bias from M1 falls on the terminal phase, so it propagates directly into half-life, Vz and accumulation predictions.

Covariate model building

Order of operations that avoids the usual traps:

  1. Get the structural and variability models right first. A covariate added to a misspecified structural model can absorb the misspecification and look significant.
  2. Include size and maturation on mechanistic grounds, not on statistical ones. Allometric weight scaling on clearance and volume is a prior, not a hypothesis to test. Fixing the exponents at 0.75/1.0 is standard and usually preferable to estimating them.
  3. Screen on eta-covariate plots, not on the objective function, to generate candidates.
  4. Forward inclusion at p < 0.05 (ΔOFV > 3.84, 1 df), backward elimination at p < 0.001 (ΔOFV > 10.83). The asymmetry exists because forward selection on the same data inflates the false-positive rate.
  5. Prefer full-model or full random-effects approaches when the goal is to quantify a covariate effect and its uncertainty rather than to select covariates. Stepwise selection produces biased effect sizes (selection bias towards large effects) and confidence intervals that are too narrow.

The likelihood-ratio test relies on the ΔOFV being chi-square distributed, which is only approximately true, and is not true on the boundary — testing whether a variance component is zero puts the null on the edge of the parameter space, and the nominal p-value is conservative.

Correlated covariates (weight and BMI; age and renal function) cannot both be included informatively. Collinearity shows up as an inflated condition number and unstable estimates rather than as a failure.

Diagnostics

Goodness-of-fit plots. DV vs PRED and DV vs IPRED; CWRES vs time and vs PRED; |IWRES| vs IPRED. Use CWRES, not WRES — WRES is computed under the FO approximation and is misleading for a FOCE-estimated model. Trends in CWRES against time indicate structural misspecification; a fan shape against PRED indicates the residual error model.

Shrinkage. eta shrinkage = 1 − SD(eta_i)/omega. Above roughly 20-30%, individual estimates have collapsed towards the population mean, and:

  • eta-covariate plots become uninformative and can show spurious relationships;
  • IPRED-based diagnostics look artificially good, because IPRED is being pulled towards the data;
  • individual parameter estimates should not be used for secondary analysis.

Epsilon shrinkage above ~30% makes IWRES-based diagnostics unreliable for the same reason. FDA's 2022 population pharmacokinetics guidance states that model selection based on shrinkage is not necessary — shrinkage is a caveat on how you may interpret diagnostics, not a selection criterion.

Visual predictive check. Simulate many replicates of the study design; plot observed percentiles against the simulated prediction intervals for those percentiles. Use prediction-corrected VPC whenever doses or covariates differ across subjects, otherwise the between-subject spread in predictions swamps the comparison. A VPC evaluates the whole model — structural, variability and residual — and is the single most informative diagnostic.

NPDE are the recommended numerical counterpart: decorrelated, and should be N(0,1) under a correct model. Test mean, variance, and normality, and plot against time and predictions.

Parameter uncertainty. The $COVARIANCE sandwich estimator is standard but fails on overparameterised models. Alternatives: nonparametric bootstrap (expensive; also fails when the model is unstable, which is informative in itself), log-likelihood profiling (best for a single poorly determined parameter), and sampling importance resampling (SIR), which is much cheaper than bootstrap and works when the covariance step fails.

Condition number — the ratio of largest to smallest eigenvalue of the correlation matrix of the estimates. Above about 1000 indicates ill-conditioning and unreliable standard errors. fit_compartmental.py reports this quantity using the same definition, so the familiar threshold applies.

Model evaluation checklist

  • Minimisation successful, and the covariance step completed
  • Parameter RSEs: fixed effects below ~30%, variance components below ~50%
  • No parameter at a bound
  • Condition number below 1000
  • Eta shrinkage below 30% for any eta used in a covariate plot
  • CWRES without trend against time or PRED
  • pcVPC with observed percentiles inside the simulated intervals
  • The model reproduces the quantity the analysis exists to predict — not just the observations

What a population analysis has to report

FDA's 2022 guidance expects the analysis plan to be prospective and the report to state: the data (including exclusions and how BLQ was handled), the structural and statistical models with justification, the covariate strategy defined before analysis, the estimation method and software version, the diagnostics, and the model's intended use. Deviations from the plan are documented, not silently absorbed.