skills/pkpd-modeling/references/population-pk.md
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%.
| Method | Character | When |
|---|---|---|
| FO | First-order linearisation about eta = 0 | Obsolete for final models; biased with large IIV. Still useful for initial estimates |
| FOCE | Linearises about the individual eta | The workhorse |
FOCE-I (INTERACTION) | FOCE with eta-epsilon interaction | Required whenever the residual error is proportional or combined. Omitting INTERACTION with a proportional error model is a common and consequential mistake |
| Laplace | Second-order expansion | Needed for non-continuous data (categorical, count, time-to-event) and for LIKELIHOOD/-2LL models |
| SAEM | Stochastic approximation EM | Robust 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 integration | Accurate objective function; usually run after SAEM |
| MCMC / NUTS | Full Bayesian | NONMEM 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.
Beal's methods, and what they actually do:
| Method | Treatment | Verdict |
|---|---|---|
| M1 | Discard all BLQ observations | Biased upward whenever a meaningful fraction is BLQ; acceptable only when that fraction is small |
| M2 | Discard BLQ, condition the likelihood on being above LLOQ | Better than M1, rarely used |
| M3 | Maximum likelihood: BLQ contributes the probability that the observation is below LLOQ | The reference method. Needs F_FLAG=1 and the Laplace method |
| M4 | M3 conditioned on the concentration being positive | Marginal improvement over M3 |
| M5 | Substitute LLOQ/2 | Biased; still common; at least state it |
| M6 | Substitute LLOQ/2 for the first BLQ in a run, discard the rest | Ad hoc |
| M7 | Substitute 0 | Worst; 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.
Order of operations that avoids the usual traps:
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.
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:
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.
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.