Covariate Modeling in Population Pharmacokinetics: Intuition, Diagnostics, and Methodology

Lecture Notes

1 1. The Role of Covariates: Explaining Variability

1.1 1.1 Intuitive Foundation

The fundamental premise of covariate modeling is that covariates explain between-subject variability. To build this intuition, consider two scenarios:

Scenario A — Undifferentiated population: When concentration-time data from all subjects appears as a single cloud of observations, the between-subject variability in clearance is captured entirely by the random effect term:

\[CL_i = \theta_{CL} \cdot e^{\eta_{CL,i}}\]

where $\eta_{CL,i} \sim \mathcal{N}(0, \omega_{CL}^2)$. If the coefficient of variation of clearance across subjects is approximately 40%, the standard deviation $\omega_{CL}$ encodes all of that unexplained variability.

Scenario B — Two distinct subpopulations: When the data visibly separates into two groups with different pharmacokinetic profiles (e.g., faster vs. slower clearance), a single population parameter is no longer adequate. We can instead parameterize each group explicitly:

\[CL_i = \theta_{CL,\text{group1}} \cdot e^{\eta_{CL,i}} \quad \text{(Group 1)}\] \[CL_i = \theta_{CL,\text{group2}} \cdot e^{\eta_{CL,i}} \quad \text{(Group 2)}\]

Two key observations from comparing these equations to the single-group model:

The population clearance differs between groups.
The variability is shared — the distribution from which $\eta$ is drawn is the same for both groups (same $\omega_{CL}$).

1.2 1.2 The Fundamental Property of Individual Parameters

A critical principle: $\eta$ will be whatever it needs to be to fit the data. The individual parameter estimate is fixed by the observed data. Given the individual’s true clearance, the $\eta$ simply adjusts to the difference between that individual value and whatever the current population parameter is.

Numerical illustration:

Setting	$\theta_{CL}$	$CL_{\text{subj 1}}$	$CL_{\text{subj 11}}$	$\eta_{\text{subj 1}}$	$\eta_{\text{subj 11}}$
One group	10 L/h	9	11	−1	+1
Two groups	7 (Grp 1), 12 (Grp 2)	9	11	+2	−1

The individual clearances do not change; $\eta$ simply shifts to accommodate the new population parameter. As a consequence, the variability per group decreases — the $\omega_{CL}$ that was 40% for the single model may drop to approximately 20% after the group split. This reduction represents the variability that the covariate has explained.

1.3 1.3 Equivalent Parameterizations for a Binary Covariate

For a binary covariate (e.g., two groups), there are three equivalent ways to write the clearance model:

Single population (no covariate): \[CL_i = \theta_{CL} \cdot e^{\eta_{CL,i}}\]
Separate parameters per group: \[CL_i = \theta_{CL,\text{group1}} \cdot e^{\eta_{CL,i}} \quad \text{or} \quad \theta_{CL,\text{group2}} \cdot e^{\eta_{CL,i}}\]
Proportional covariate effect (single base + offset): \[CL_i = \theta_{CL} \cdot \left(1 + \theta_{\Delta CL} \cdot \mathbf{1}[\text{group}=1]\right) \cdot e^{\eta_{CL,i}}\]

Parameterizations 2 and 3 each introduce one additional parameter relative to parameterization 1 (i.e., one additional degree of freedom). The choice between them is a communication decision: whether the audience is better served by absolute group clearances or by a proportional difference from a reference group.

The indicator notation $\mathbf{1}[\text{group}=1]$ evaluates to 1 when the condition is true and 0 otherwise — standard behavior in any programming language. When the condition is false, the term collapses to zero and the equation reduces to $\theta_{CL}$, the reference group clearance.

2 2. Empirical Bayes Estimates and Covariate Diagnostics

2.1 2.1 The Linear Interpretation of $\eta$

Taking the log of the base model gives:

\[\ln CL_i = \ln \theta_{CL} + \eta_{CL,i}\]

This is a linear equation of the form $Y = mX + b$, where: - $Y = \eta_{CL,i}$ (or equivalently, $\ln CL_i$) - The slope $m = 0$ (no systematic trend against any covariate) - The intercept $b = 0$ (mean of $\eta$ is zero by assumption)

The distribution assumption $\eta_{CL} \sim \mathcal{N}(0, \omega_{CL}^2)$ must hold across all conditions.

2.2 2.2 Empirical Bayes Estimates vs. Covariate Plots

After fitting the base model, the empirical Bayes estimates (EBEs) — the individual $\hat{\eta}$ values — are extracted and regressed against candidate covariates. The diagnostic expectation is that a plot of $\hat{\eta}_{CL}$ versus any covariate should show:

Mean near zero (no bias)
Slope near zero (no trend)
Even spread (homogeneous variability)

If a covariate $X$ produces a systematic trend (e.g., a positive slope), it indicates that $\eta_{CL}$ contains unexplained structure attributable to $X$, and the covariate should be incorporated into the structural model:

\[CL_i = \theta_{CL} \cdot \left(\frac{X_i}{X_{\text{ref}}}\right)^{\theta_X} \cdot e^{\eta_{CL,i}}\]

After including the covariate, the $\eta_{CL}$ vs. $X$ plot should return to a slope of zero, confirming that the relationship has been captured.

2.3 2.3 $\eta$-Shrinkage and Its Consequences

Shrinkage occurs when subjects have sparse or informative-poor data. In such cases, the individual estimates are “pulled” toward the population mean (zero), resulting in a compressed distribution of $\hat{\eta}$ values.

When shrinkage is substantial (generally greater than 30%), the EBE-versus-covariate plots become unreliable in two ways:

True correlations may be hidden — a covariate relationship that genuinely exists in the data will be obscured because the $\hat{\eta}$ values have collapsed toward zero.
Spurious correlations may appear — artificial structure can emerge from the shrinkage pattern rather than from a true biological relationship.

Therefore, before interpreting EBE-covariate diagnostic plots, $\eta$-shrinkage must be assessed, and plots should be interpreted with caution when shrinkage exceeds approximately 30%.

3 3. Covariate Functional Forms

3.1 3.1 Binary Covariates

As described above, a binary covariate (two categories) introduces one additional parameter. Using indicator variable notation:

\[CL_i = \theta_{CL} \cdot \left(1 + \theta_{\text{sex}} \cdot \mathbf{1}[\text{sex} = \text{male}]\right) \cdot e^{\eta_{CL,i}}\]

3.2 3.2 Categorical Covariates with More Than Two Levels

For an $N$-level categorical covariate, $N - 1$ coefficients are required. This is directly analogous to the one-hot encoding scheme used in machine learning and survival analysis.

Example — three-level covariate (low, medium, high):

Create three binary indicator columns and include $N - 1 = 2$ of them:

\[CL_i = \theta_{CL} \cdot \left(1 + \theta_{\text{low}} \cdot \mathbf{1}[\text{level} = \text{low}] + \theta_{\text{med}} \cdot \mathbf{1}[\text{level} = \text{medium}]\right) \cdot e^{\eta_{CL,i}}\]

The omitted level (here, “high”) serves as the reference category.

3.3 3.3 Continuous Covariates

Continuous covariates should never be discretized for modeling purposes. Categorizing a continuous variable introduces substantial bias and reduces predictive power. The model can then only predict within the observed range of each category, and extrapolation beyond those bins is invalid.

If a continuous covariate appears to exhibit a non-linear or segmented relationship with $\hat{\eta}$, the correct approach is to develop a continuous mathematical equation (e.g., a power function, polynomial, or piecewise function) that captures the shape across the full covariate range — provided data exist across that range. Categorization is acceptable only as a post-hoc communication tool after continuous modeling is complete.

4 4. Covariate Selection Methodology

Three principal methodologies are used for covariate selection in population pharmacokinetic modeling.

4.1 4.1 Generalized Additive Modeling (GAM)

GAM operates outside the full pharmacokinetic model, using the extracted EBEs as the dependent variable.

Procedure:

For each candidate covariate $X_k$, fit a simple linear regression: \[\hat{\eta}_{CL} = m \cdot X_k + b\]
Select the covariate with the best fit metric (e.g., lowest AIC) — the “winner” of round one.
Condition on the winner and repeat: test all remaining covariates in combination with the winner.
Continue iterating until no additional covariate further reduces the metric.

This process is run separately for each random-effect parameter ($\eta_{CL}$, $\eta_V$, etc.). The result is a ranked set of covariate-parameter relationships that informs the subsequent pharmacokinetic modeling step.

GAM is computationally efficient because it uses simple linear regressions rather than full nonlinear mixed-effects model fits.

4.2 4.2 Stepwise Covariate Modeling (SCM)

SCM performs covariate selection within the full pharmacokinetic model and is the most widely used approach in practice.

4.2.1 Forward Addition

Starting from the base model (no covariates):

For each candidate covariate $X_k$, fit the full PK model with $X_k$ added to the parameter of interest.
Compare each candidate model against the base using a selection criterion (AIC, BIC, or likelihood ratio test).
Select the best-performing covariate — the “winner.”
Condition on the winner, repeat step 1 for remaining covariates.
Continue until no further improvement is achieved.

4.2.2 Backward Elimination

Starting from a model containing all covariates under consideration:

Remove one covariate at a time and re-fit.
If removing a covariate causes a large increase in AIC (i.e., the model worsens substantially), that covariate is important and retained.
The covariate whose removal causes the smallest AIC increase (or no increase) is eliminated.
Repeat until all remaining covariates are confirmed as important.

4.2.3 Combined Forward-Backward Approach

The recommended practice is to perform forward addition followed by backward elimination. The rationale is:

Forward addition can be somewhat liberal (less strict significance threshold, e.g., $p < 0.05$).
Backward elimination should be more conservative (e.g., $p < 0.001$) to confirm that retained covariates are truly necessary.

Each step in SCM constitutes a hypothesis test (null hypothesis: the covariate has no effect on the parameter), which raises concerns about alpha spending from multiple comparisons. In practice, Bonferroni or similar corrections are rarely applied in SCM, which is a known limitation of the methodology.

4.3 4.3 Full Fixed Effects Modeling (FFM)

In FFM, all biologically plausible covariates are included in the model simultaneously from the outset, and no hypothesis testing is performed.

Procedure:

Identify all covariates that have physiologically plausible relationships with each PK parameter, based on subject-matter expertise.
Remove redundant covariates before model fitting (e.g., body weight and body surface area are strongly correlated — include only one).
Fit the model with all selected covariates and obtain parameter estimates with confidence intervals.
Assess clinical (or operational) significance of each covariate coefficient based on the magnitude of its effect, not its p-value.

Key distinction from SCM/GAM: In FFM, statistical significance (i.e., whether the confidence interval excludes zero) is a necessary but not sufficient criterion for calling a covariate “important.” A coefficient may be estimated precisely (CI does not include zero) while being clinically irrelevant (e.g., a 10% effect on clearance). Conversely, a covariate may be clinically meaningful but imprecisely estimated due to limited sample size — which reflects study design, not biology.

FFM is technically more rigorous because: - It avoids the multiple-testing problem inherent in SCM. - It allows the modeler to exercise physiological judgment in covariate selection. - It requires fitting only one model (assuming covariate cleaning is done beforehand) rather than dozens or hundreds.

The primary challenge is that covariate cleaning — identifying and removing correlated covariates — requires careful prior analysis and domain knowledge.

5 5. Model Comparison for Covariate Models

The same likelihood-based framework used for structural model selection applies to covariate model comparison. When adding one covariate (one degree of freedom):

The likelihood ratio test statistic follows a $\chi^2$ distribution with 1 degree of freedom.
The critical value at $\alpha = 0.05$ is 3.84.

A covariate model is considered superior to the base model if it improves the objective function by more than this threshold. Graphical diagnostics (EBE plots, goodness-of-fit plots) must also confirm the improvement.

The iterative process is: base model → add one covariate → compare → retain or reject → add next covariate → repeat.

6 6. Summary

The key concepts covered in this lecture are:

Covariates explain variability — adding covariates shifts unexplained between-subject variability from $\omega$ into the structural model.
$\eta$ adjusts automatically — individual ETAs absorb whatever residual difference exists between an individual’s true parameter and the population prediction, regardless of parameterization.
EBE-covariate plots are diagnostic tools — they test whether the assumption $\mathbb{E}[\eta] = 0$ holds across covariate strata. Systematic trends indicate a missing covariate relationship.
Shrinkage invalidates EBE diagnostics — when $\eta$-shrinkage exceeds approximately 30%, EBE-covariate plots are unreliable.
Never discretize continuous covariates during modeling — this introduces bias and destroys predictive power.
Three covariate selection methods — GAM (regression-based, outside the PK model), SCM (stepwise hypothesis testing within the full PK model), and FFM (all-inclusive, no hypothesis testing, judged by clinical significance).

--- title: "Covariate Modeling in Population Pharmacokinetics: Intuition, Diagnostics, and Methodology" subtitle: "Lecture Notes" format: html: toc: true toc-depth: 3 --- --- ## 1. The Role of Covariates: Explaining Variability ### 1.1 Intuitive Foundation The fundamental premise of covariate modeling is that covariates explain between-subject variability. To build this intuition, consider two scenarios: **Scenario A — Undifferentiated population:** When concentration-time data from all subjects appears as a single cloud of observations, the between-subject variability in clearance is captured entirely by the random effect term: $$CL_i = \theta_{CL} \cdot e^{\eta_{CL,i}}$$ where $\eta_{CL,i} \sim \mathcal{N}(0, \omega_{CL}^2)$. If the coefficient of variation of clearance across subjects is approximately 40%, the standard deviation $\omega_{CL}$ encodes all of that unexplained variability. **Scenario B — Two distinct subpopulations:** When the data visibly separates into two groups with different pharmacokinetic profiles (e.g., faster vs. slower clearance), a single population parameter is no longer adequate. We can instead parameterize each group explicitly: $$CL_i = \theta_{CL,\text{group1}} \cdot e^{\eta_{CL,i}} \quad \text{(Group 1)}$$ $$CL_i = \theta_{CL,\text{group2}} \cdot e^{\eta_{CL,i}} \quad \text{(Group 2)}$$ Two key observations from comparing these equations to the single-group model: 1. The **population clearance differs** between groups. 2. The **variability is shared** — the distribution from which $\eta$ is drawn is the same for both groups (same $\omega_{CL}$). ### 1.2 The Fundamental Property of Individual Parameters A critical principle: **$\eta$ will be whatever it needs to be to fit the data.** The individual parameter estimate is fixed by the observed data. Given the individual's true clearance, the $\eta$ simply adjusts to the difference between that individual value and whatever the current population parameter is. **Numerical illustration:** | Setting | $\theta_{CL}$ | $CL_{\text{subj 1}}$ | $CL_{\text{subj 11}}$ | $\eta_{\text{subj 1}}$ | $\eta_{\text{subj 11}}$ | |---|---|---|---|---|---| | One group | 10 L/h | 9 | 11 | −1 | +1 | | Two groups | 7 (Grp 1), 12 (Grp 2) | 9 | 11 | +2 | −1 | The individual clearances do not change; $\eta$ simply shifts to accommodate the new population parameter. As a consequence, the **variability per group decreases** — the $\omega_{CL}$ that was 40% for the single model may drop to approximately 20% after the group split. This reduction represents the variability that the covariate has explained. ### 1.3 Equivalent Parameterizations for a Binary Covariate For a binary covariate (e.g., two groups), there are three equivalent ways to write the clearance model: 1. **Single population (no covariate):** $$CL_i = \theta_{CL} \cdot e^{\eta_{CL,i}}$$ 2. **Separate parameters per group:** $$CL_i = \theta_{CL,\text{group1}} \cdot e^{\eta_{CL,i}} \quad \text{or} \quad \theta_{CL,\text{group2}} \cdot e^{\eta_{CL,i}}$$ 3. **Proportional covariate effect (single base + offset):** $$CL_i = \theta_{CL} \cdot \left(1 + \theta_{\Delta CL} \cdot \mathbf{1}[\text{group}=1]\right) \cdot e^{\eta_{CL,i}}$$ Parameterizations 2 and 3 each introduce **one additional parameter** relative to parameterization 1 (i.e., one additional degree of freedom). The choice between them is a communication decision: whether the audience is better served by absolute group clearances or by a proportional difference from a reference group. The indicator notation $\mathbf{1}[\text{group}=1]$ evaluates to 1 when the condition is true and 0 otherwise — standard behavior in any programming language. When the condition is false, the term collapses to zero and the equation reduces to $\theta_{CL}$, the reference group clearance. --- ## 2. Empirical Bayes Estimates and Covariate Diagnostics ### 2.1 The Linear Interpretation of $\eta$ Taking the log of the base model gives: $$\ln CL_i = \ln \theta_{CL} + \eta_{CL,i}$$ This is a linear equation of the form $Y = mX + b$, where: - $Y = \eta_{CL,i}$ (or equivalently, $\ln CL_i$) - The slope $m = 0$ (no systematic trend against any covariate) - The intercept $b = 0$ (mean of $\eta$ is zero by assumption) The distribution assumption $\eta_{CL} \sim \mathcal{N}(0, \omega_{CL}^2)$ must hold across all conditions. ### 2.2 Empirical Bayes Estimates vs. Covariate Plots After fitting the base model, the empirical Bayes estimates (EBEs) — the individual $\hat{\eta}$ values — are extracted and regressed against candidate covariates. The diagnostic expectation is that a plot of $\hat{\eta}_{CL}$ versus any covariate should show: - **Mean near zero** (no bias) - **Slope near zero** (no trend) - **Even spread** (homogeneous variability) If a covariate $X$ produces a systematic trend (e.g., a positive slope), it indicates that $\eta_{CL}$ contains unexplained structure attributable to $X$, and the covariate should be incorporated into the structural model: $$CL_i = \theta_{CL} \cdot \left(\frac{X_i}{X_{\text{ref}}}\right)^{\theta_X} \cdot e^{\eta_{CL,i}}$$ After including the covariate, the $\eta_{CL}$ vs. $X$ plot should return to a slope of zero, confirming that the relationship has been captured. ### 2.3 $\eta$-Shrinkage and Its Consequences **Shrinkage** occurs when subjects have sparse or informative-poor data. In such cases, the individual estimates are "pulled" toward the population mean (zero), resulting in a compressed distribution of $\hat{\eta}$ values. When shrinkage is substantial (generally greater than 30%), the EBE-versus-covariate plots become unreliable in two ways: 1. **True correlations may be hidden** — a covariate relationship that genuinely exists in the data will be obscured because the $\hat{\eta}$ values have collapsed toward zero. 2. **Spurious correlations may appear** — artificial structure can emerge from the shrinkage pattern rather than from a true biological relationship. Therefore, before interpreting EBE-covariate diagnostic plots, $\eta$-shrinkage must be assessed, and plots should be interpreted with caution when shrinkage exceeds approximately 30%. --- ## 3. Covariate Functional Forms ### 3.1 Binary Covariates As described above, a binary covariate (two categories) introduces one additional parameter. Using indicator variable notation: $$CL_i = \theta_{CL} \cdot \left(1 + \theta_{\text{sex}} \cdot \mathbf{1}[\text{sex} = \text{male}]\right) \cdot e^{\eta_{CL,i}}$$ ### 3.2 Categorical Covariates with More Than Two Levels For an $N$-level categorical covariate, $N - 1$ coefficients are required. This is directly analogous to the one-hot encoding scheme used in machine learning and survival analysis. **Example — three-level covariate (low, medium, high):** Create three binary indicator columns and include $N - 1 = 2$ of them: $$CL_i = \theta_{CL} \cdot \left(1 + \theta_{\text{low}} \cdot \mathbf{1}[\text{level} = \text{low}] + \theta_{\text{med}} \cdot \mathbf{1}[\text{level} = \text{medium}]\right) \cdot e^{\eta_{CL,i}}$$ The omitted level (here, "high") serves as the reference category. ### 3.3 Continuous Covariates Continuous covariates should **never be discretized** for modeling purposes. Categorizing a continuous variable introduces substantial bias and reduces predictive power. The model can then only predict within the observed range of each category, and extrapolation beyond those bins is invalid. If a continuous covariate appears to exhibit a non-linear or segmented relationship with $\hat{\eta}$, the correct approach is to develop a continuous mathematical equation (e.g., a power function, polynomial, or piecewise function) that captures the shape across the full covariate range — provided data exist across that range. Categorization is acceptable **only as a post-hoc communication tool** after continuous modeling is complete. --- ## 4. Covariate Selection Methodology Three principal methodologies are used for covariate selection in population pharmacokinetic modeling. ### 4.1 Generalized Additive Modeling (GAM) GAM operates outside the full pharmacokinetic model, using the extracted EBEs as the dependent variable. **Procedure:** 1. For each candidate covariate $X_k$, fit a simple linear regression: $$\hat{\eta}_{CL} = m \cdot X_k + b$$ 2. Select the covariate with the best fit metric (e.g., lowest AIC) — the "winner" of round one. 3. Condition on the winner and repeat: test all remaining covariates in combination with the winner. 4. Continue iterating until no additional covariate further reduces the metric. This process is run separately for each random-effect parameter ($\eta_{CL}$, $\eta_V$, etc.). The result is a ranked set of covariate-parameter relationships that informs the subsequent pharmacokinetic modeling step. GAM is computationally efficient because it uses simple linear regressions rather than full nonlinear mixed-effects model fits. ### 4.2 Stepwise Covariate Modeling (SCM) SCM performs covariate selection within the full pharmacokinetic model and is the most widely used approach in practice. #### Forward Addition Starting from the base model (no covariates): 1. For each candidate covariate $X_k$, fit the full PK model with $X_k$ added to the parameter of interest. 2. Compare each candidate model against the base using a selection criterion (AIC, BIC, or likelihood ratio test). 3. Select the best-performing covariate — the "winner." 4. Condition on the winner, repeat step 1 for remaining covariates. 5. Continue until no further improvement is achieved. #### Backward Elimination Starting from a model containing all covariates under consideration: 1. Remove one covariate at a time and re-fit. 2. If removing a covariate causes a large increase in AIC (i.e., the model worsens substantially), that covariate is important and retained. 3. The covariate whose removal causes the smallest AIC increase (or no increase) is eliminated. 4. Repeat until all remaining covariates are confirmed as important. #### Combined Forward-Backward Approach The recommended practice is to perform **forward addition followed by backward elimination**. The rationale is: - Forward addition can be somewhat liberal (less strict significance threshold, e.g., $p < 0.05$). - Backward elimination should be more conservative (e.g., $p < 0.001$) to confirm that retained covariates are truly necessary. Each step in SCM constitutes a hypothesis test (null hypothesis: the covariate has no effect on the parameter), which raises concerns about **alpha spending** from multiple comparisons. In practice, Bonferroni or similar corrections are rarely applied in SCM, which is a known limitation of the methodology. ### 4.3 Full Fixed Effects Modeling (FFM) In FFM, all biologically plausible covariates are included in the model simultaneously from the outset, and no hypothesis testing is performed. **Procedure:** 1. Identify all covariates that have physiologically plausible relationships with each PK parameter, based on subject-matter expertise. 2. Remove redundant covariates before model fitting (e.g., body weight and body surface area are strongly correlated — include only one). 3. Fit the model with all selected covariates and obtain parameter estimates with confidence intervals. 4. Assess **clinical (or operational) significance** of each covariate coefficient based on the magnitude of its effect, not its p-value. **Key distinction from SCM/GAM:** In FFM, statistical significance (i.e., whether the confidence interval excludes zero) is a necessary but not sufficient criterion for calling a covariate "important." A coefficient may be estimated precisely (CI does not include zero) while being clinically irrelevant (e.g., a 10% effect on clearance). Conversely, a covariate may be clinically meaningful but imprecisely estimated due to limited sample size — which reflects study design, not biology. FFM is technically more rigorous because: - It avoids the multiple-testing problem inherent in SCM. - It allows the modeler to exercise physiological judgment in covariate selection. - It requires fitting only one model (assuming covariate cleaning is done beforehand) rather than dozens or hundreds. The primary challenge is that covariate cleaning — identifying and removing correlated covariates — requires careful prior analysis and domain knowledge. --- ## 5. Model Comparison for Covariate Models The same likelihood-based framework used for structural model selection applies to covariate model comparison. When adding one covariate (one degree of freedom): - The likelihood ratio test statistic follows a $\chi^2$ distribution with 1 degree of freedom. - The critical value at $\alpha = 0.05$ is 3.84. A covariate model is considered superior to the base model if it improves the objective function by more than this threshold. Graphical diagnostics (EBE plots, goodness-of-fit plots) must also confirm the improvement. The iterative process is: base model → add one covariate → compare → retain or reject → add next covariate → repeat. --- ## 6. Summary The key concepts covered in this lecture are: 1. **Covariates explain variability** — adding covariates shifts unexplained between-subject variability from $\omega$ into the structural model. 2. **$\eta$ adjusts automatically** — individual ETAs absorb whatever residual difference exists between an individual's true parameter and the population prediction, regardless of parameterization. 3. **EBE-covariate plots are diagnostic tools** — they test whether the assumption $\mathbb{E}[\eta] = 0$ holds across covariate strata. Systematic trends indicate a missing covariate relationship. 4. **Shrinkage invalidates EBE diagnostics** — when $\eta$-shrinkage exceeds approximately 30%, EBE-covariate plots are unreliable. 5. **Never discretize continuous covariates** during modeling — this introduces bias and destroys predictive power. 6. **Three covariate selection methods** — GAM (regression-based, outside the PK model), SCM (stepwise hypothesis testing within the full PK model), and FFM (all-inclusive, no hypothesis testing, judged by clinical significance).