Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › Interactions Between Covariates In PD Models

Population PK/PD Model Building · Tutorial 439 of 598

Interactions Between Covariates In PD Models

Specify and interpret covariate interactions on PD parameters, calculate their conditional effects, and recognize when the data cannot reliably estimate them.

Advanced 11 min read

What You'll Learn

  • Write an interaction on a transformed PD-parameter scale and identify its reference values.
  • Calculate how a covariate effect changes across levels of a second covariate.
  • Interpret interactions for continuous and categorical covariates.
  • Translate interaction coefficients into group-specific parameter values and ratios.
  • Check whether the observed covariate combinations contain enough variation to estimate an interaction.
  • Recognize scale dependence, sparse-data problems, and common misinterpretations.

Why Covariate Interactions Matter

An additive covariate model assumes that the effect of one covariate on a PD parameter is the same at every value of another covariate. That may be inadequate when, for example, the relationship between age and a PD parameter differs by genotype, or when the contrast between two treatment groups varies with a continuous covariate. An interaction represents that change in one covariate effect across values of another.

An interaction is a feature of the specified model scale. A zero interaction on the log-parameter scale means that covariate effects combine additively on that scale, or multiplicatively on the parameter scale. It does not mean that the covariates have no effects, nor does it establish a biological mechanism. Interpretation therefore starts by stating both the PD parameter and the scale on which its covariates enter.

Definition: A covariate interaction is a model term that makes the effect of one covariate on a PD parameter depend on the value of another covariate. For a positive typical PD parameter \(\theta_{\mathrm{typ}}\), a log-linear interaction between covariates \(x\) and \(z\) can be written as \(\log(\theta_{\mathrm{typ}})=\beta_0+\beta_x x+\beta_z z+\beta_{xz}xz\). The coefficient \(\beta_{xz}\) quantifies how the conditional slope for \(x\), on the log-parameter scale, changes per unit increase in \(z\).

Model Form and Conditional Effects

For a positive parameter, the log link ensures that the typical parameter remains positive. Covariates should be coded and scaled so that their reference values and units are explicit. If \(x\) and \(z\) are centered, zero denotes their chosen centering values; if a categorical covariate is coded as zero and one, zero denotes its reference category. The intercept \(\beta_0\) is the log typical parameter at \(x=0,z=0\).

$$ \log(\theta_{\mathrm{typ}})=\beta_0+\beta_x x+\beta_z z+\beta_{xz}xz $$

Holding \(z\) fixed, the slope for \(x\) on the log scale is \(\beta_x+\beta_{xz}z\). Thus \(\beta_x\) is the slope only when \(z=0\), and \(\beta_{xz}\) is the change in that slope per unit increase in \(z\). Similarly, the conditional slope for \(z\) is \(\beta_z+\beta_{xz}x\). On the parameter scale, a change \(\Delta x\) at fixed \(z\) multiplies \(\theta_{\mathrm{typ}}\) by \(\exp\{(\beta_x+\beta_{xz}z)\Delta x\}\).

The product term must be interpreted together with its corresponding main effects. Retaining \(\beta_{xz}xz\) while omitting \(\beta_xx\) or \(\beta_zz\) generally changes the model's reference interpretation and imposes an additional restriction. Unless there is a compelling parameterization reason, include both main effects when fitting their interaction.

Worked Example: Age Effect That Differs by Renal-Function Group

Suppose a hypothetical population PD model describes a positive typical \(EC_{50}\), in mg/L. Let \(x\) be age in decades above 60 years and let \(z=1\) indicate a specified renal-function group, with \(z=0\) the reference group. Suppose the fitted covariate component is \(\log(\theta_{\mathrm{typ}})=\log(4)+0.12x+0.20z-0.15xz\). Calculate the age effect in each group and the typical parameter at \(x=2\) in each group.

At the reference value \(z=0\), the age slope on the log scale is \(0.12\) per decade. A one-decade increase multiplies \(EC_{50}\) by \(\exp(0.12)\approx1.1275\), an increase of about \(12.75\%\).

At \(z=1\), the age slope is \(0.12-0.15=-0.03\) per decade. A one-decade increase multiplies the parameter by \(\exp(-0.03)\approx0.9704\), or a decrease of about \(2.96\%\). The interaction does not say that the age effect is \(-0.15\): that is the difference in age slopes between the groups on the log scale.

For \(x=2,z=0\), the typical value is \(\exp\{\log(4)+0.12(2)\}=4\exp(0.24)=5.085\ \mathrm{mg/L}\). For \(x=2,z=1\), it is \(\exp\{\log(4)+0.12(2)+0.20-0.15(2)\}=4\exp(0.14)=4.601\ \mathrm{mg/L}\).

At \(x=2\), the ratio of the group-specific typical parameters is \(\exp(0.20-0.15(2))=\exp(-0.10)\approx0.9048\). Thus the group coded one has a typical \(EC_{50}\) about \(9.52\%\) lower than the reference group at that age, according to this model. The group contrast changes with age; it is not the same as the coefficient \(0.20\), which is the log-scale contrast only at \(x=0\).

Categorical Covariates and the Difference in Differences

With two binary covariates \(A\) and \(B\), coded zero for their reference categories and one for the alternative categories, the same model gives four group-specific log parameters. The interaction coefficient is the difference between the \(A\) contrast when \(B=1\) and the \(A\) contrast when \(B=0\). Equivalently, it is the difference between the \(B\) contrasts across levels of \(A\). This is often called a difference in differences on the model's linear-predictor scale.

The identity follows directly by evaluating the model in each of its four cells. Write \(\eta_{ab}=\log(\theta_{\mathrm{typ}})\) at \(A=a,B=b\). Then the difference between the two \(A\) contrasts is:

$$ (\eta_{11}-\eta_{01})-(\eta_{10}-\eta_{00})=\beta_{AB} $$

To verify it, substitute \(\eta_{ab}=\beta_0+\beta_Aa+\beta_Bb+\beta_{AB}ab\). The first contrast is \(\beta_A+\beta_{AB}\), the second is \(\beta_A\), and their difference is \(\beta_{AB}\). On the original parameter scale, exponentiating \(\beta_{AB}\) gives a ratio of ratios, not a difference of parameter values.

Worked Example: Two Binary Covariates on a PD Parameter

Consider a hypothetical model for typical \(EC_{50}\) in mg/L, with \(A\) indicating a concomitant medication and \(B\) indicating a genotype category. Suppose \(\log(\theta_{\mathrm{typ}})=\log(10)+0.10A+0.30B+0.25AB\). Calculate the typical parameter in each covariate group and interpret the interaction.

For \(A=0,B=0\), \(\theta_{\mathrm{typ}}=\exp\{\log(10)\}=10.00\ \mathrm{mg/L}\). For \(A=1,B=0\), it is \(10\exp(0.10)=11.05\ \mathrm{mg/L}\). For \(A=0,B=1\), it is \(10\exp(0.30)=13.50\ \mathrm{mg/L}\). For \(A=1,B=1\), it is \(10\exp(0.10+0.30+0.25)=10\exp(0.65)=19.16\ \mathrm{mg/L}\).

When \(B=0\), the \(A=1\) versus \(A=0\) ratio is \(\exp(0.10)=1.105\). When \(B=1\), that ratio is \(\exp(0.10+0.25)=\exp(0.35)\approx1.419\). The ratio of these two ratios is \(\exp(0.25)\approx1.284\), consistent with the interaction coefficient: the \(A\) contrast is 1.284 times as large on the multiplicative parameter scale in the \(B=1\) category.

The coefficient \(0.25\) is not a 25% increase in \(EC_{50}\) for every person with both covariates. It is a log-scale difference in contrasts; the corresponding ratio of ratios is about 1.284. The resulting group-specific values also depend on both main-effect coefficients and the intercept.

Continuous-by-Continuous Interactions

When both covariates are continuous, the product term still changes conditional slopes. Its units are the log-parameter scale divided by the product of the covariate units. Centering can make the main effects more interpretable: with \(x\) and \(z\) centered, \(\beta_x\) is the slope for \(x\) at the centering value of \(z\), rather than at a potentially irrelevant value of zero. Centering does not remove the interaction or create information absent from the data.

Worked Example: Continuous Covariates and Conditional Slopes

Suppose, hypothetically, that a positive typical PD parameter in arbitrary units has model \(\log(\theta_{\mathrm{typ}})=\log(2)+0.08x+0.05z+0.04xz\), where \(x\) is age in decades above its reference value and \(z\) is a centered laboratory covariate measured in its stated units. Calculate the age slope at \(z=-1\) and \(z=1\), and predict the typical parameter at \(x=2,z=1\).

The conditional age slope is \(0.08+0.04z\). At \(z=-1\), it is \(0.04\) per decade; at \(z=1\), it is \(0.12\) per decade. Thus a one-decade increase in age multiplies the parameter by \(\exp(0.04)\approx1.041\) at \(z=-1\), and by \(\exp(0.12)\approx1.1275\) at \(z=1\).

At \(x=2,z=1\), the linear predictor beyond the intercept is \(0.08(2)+0.05(1)+0.04(2)(1)=0.29\). The prediction is \(2\exp(0.29)=2.673\) parameter units. These are conditional model predictions for the stated covariate values, not evidence that changing age or the laboratory covariate would cause the parameter to change.

What Data Support Is Needed?

An interaction is estimable only if the covariate data contain information that distinguishes the product term from the intercept and main effects. For two binary covariates, observations are needed in all four combinations to estimate the full intercept-plus-main-effects-plus-interaction parameterization without relying on additional restrictions. If one cell is empty, the corresponding design matrix has fewer than four distinct rows and cannot have rank four; the four coefficients are not separately identifiable from those data alone.

For continuous covariates, observations should cover a meaningful range of their combinations. If \(x\) and \(z\) are strongly associated, the product column \(xz\) may be nearly explained by the intercept and main-effect columns. The model may then return an interaction estimate with large uncertainty, strong parameter correlations, or sensitivity to a few observations. A large total sample does not solve a lack of support in particular combinations.

Before interpreting an interaction, inspect covariate ranges and cross-tabulations, the number of subjects contributing to each relevant region or category, and the uncertainty and stability of the fitted coefficients. Confirm that the PD data are informative about the parameter itself; a well-populated covariate table cannot compensate for an inadequately identified PD parameter. Predictions at combinations with little or no observed support are extrapolations and should not be treated as equally reliable.

1
State the parameter and scale.
Specify which PD parameter is modeled and whether the interaction is on the parameter, log-parameter, or another link scale.
2
Define covariate coding.
Give units, centering values, and reference categories so the intercept and main effects have clear meanings.
3
Check support and estimability.
Look for empty or sparse categories, narrow joint ranges, collinearity, and weak information about the PD parameter.
4
Calculate conditional effects.
Evaluate the slope or group contrast at relevant covariate values, then transform to the parameter scale if needed.
5
Assess uncertainty and scientific relevance.
Use uncertainty estimates and model comparison alongside biological plausibility; do not equate a fitted interaction with a causal mechanism.

Diagnostics, Alternatives, and Limitations

An interaction may be going wrong if its estimate is very imprecise, changes greatly when a few subjects are removed, has a large correlation with other parameter estimates, or produces extreme predictions in sparsely observed covariate combinations. For a nested comparison of a model with versus without an interaction, the likelihood-ratio test is applicable under its regularity conditions; the models must use the same data and differ by the specified interaction terms. Statistical evidence should be considered with the uncertainty in the group-specific effects and the scientific rationale.

The main alternative is an additive covariate model without the product term. It is more parsimonious and is appropriate when a common covariate effect is scientifically reasonable and supported by the data. Separate subgroup models allow effects to differ freely, but use parameters less efficiently and may be unstable when subgroups are small. A prespecified interaction offers a structured way to represent effect modification while retaining common information; it does not guarantee that the data can estimate it precisely.

Interactions depend on scale. A model additive on the log-parameter scale implies multiplicative effects on the parameter scale, so an interaction can be present on one scale and absent on another. Centering changes the interpretation of main-effect coefficients but not the fitted values when the model is reparameterized consistently. Finally, estimates from observational covariates describe modeled associations conditional on the included covariates; they do not by themselves establish causal effect modification.

Common Mistake

A frequent error is to report \(\beta_{xz}\) as the effect of \(x\) for every value of \(z\), or as a percentage change in the PD parameter. In the log-linear model, the conditional slope for \(x\) is \(\beta_x+\beta_{xz}z\), and a one-unit change in \(x\) multiplies the parameter by \(\exp(\beta_x+\beta_{xz}z)\). The interaction coefficient describes how that conditional slope changes with \(z\); it is not, by itself, the subgroup-specific effect or a direct percent change.

Key takeaway: Interpret an interaction by evaluating the conditional effect at specified covariate values, then translating it to the PD-parameter scale. Its estimate is meaningful only when the observed data support the relevant covariate combinations and the PD parameter is itself sufficiently informed.

Check Your Understanding

Answer each question using the stated model scale and covariate coding.

  1. In \(\log(\theta_{\mathrm{typ}})=\beta_0+\beta_xx+\beta_zz+\beta_{xz}xz\), what is the slope for \(x\) when \(z=2\)?
  2. For a binary-by-binary interaction, what does exponentiating the interaction coefficient represent on the parameter scale?
  3. In a log-linear model, \(\beta_x=0.06\) and \(\beta_{xz}=-0.02\). Calculate the one-unit multiplier for \(x\) when \(z=1\).
  4. Why can an interaction be unidentifiable when one cell of a two-by-two categorical covariate table is empty?
  5. Does a statistically supported interaction establish that one covariate causally modifies the effect of another? Explain.

Solutions. (1) The conditional slope is \(\beta_x+2\beta_{xz}\). (2) It is the ratio of the two group-specific covariate ratios, also called a ratio of ratios. (3) The slope is \(0.06-0.02=0.04\), so the multiplier is \(\exp(0.04)\approx1.041\), a 4.08% increase per unit of \(x\) at \(z=1\). (4) With fewer than four observed covariate combinations, the intercept, two main-effect columns, and product column cannot all be separately identified from the observed design. (5) No. The interaction describes a scale-specific conditional association in the fitted model; causal interpretation requires additional design and assumptions.

Clinical Biostats