Time-varying covariates, overadjustment bias and mediation in pharmacokinetic/pharmacodynamic modeling

Overadjustment bias is a term used in epidemiology to refer to the situation where bias is introduced when controlling for a variable in the analysis, for example, by including the variable as a covariate in a model.^{1, 2} There are several situations in which inclusion of covariates can cause overadjustment bias, two of which are graphically illustrated in Figure 1. In the first example, the drug exposure has a certain causal effect on the outcome and part of this causal effect is mediated by an intermediate. Because part of the causal effect of exposure on outcome is explained by the effect of exposure on the intermediate, including the intermediate as a covariate in the model will result in a biased (i.e., under-estimated) estimate of the total effect of exposure on the outcome. In the second example, in Figure 1, the potential covariate is a descendant or consequence of the outcome.

The risk of overadjustment bias lies particularly in those covariates that are included in the model in a time-varying fashion or when using a time-constant, but post-baseline value as a covariate (e.g., the value observed at 3-months after treatment start in each subject). In contrast, the covariate values at baseline are not affected by treatment (assuming treatment starts at baseline) and one would therefore not expect overadjustment bias when including baseline covariates in a model.

Although overadjustment bias is rarely discussed within the context of PK/PD modeling, it is by no means less important for this setting. In a PK/PD model, the scenarios above would result in an under-estimated treatment effect. This bias could negatively impact real-life decisions. It may lead one to incorrectly conclude that the treatment is not effective (enough) or that a higher concentration may be needed to reach the target efficacy. Overadjustment bias could also result in overly conservative sample size calculations for a future clinical study.

To demonstrate the importance of overadjustment bias, we performed several illustrative simulations of a hypothetical placebo-controlled study (R code included in the Supplemental Material S1). Simulations were conducted for three scenarios involving a time-varying covariate influencing a numerical outcome variable (Figure 2a). Each scenario was simulated 1000 times, with each dataset containing 100 subjects that each have six observations for both covariate and outcome variables.

Scenario 1 assumes a direct treatment effect on the outcome, which is not mediated by the time-varying covariate (i.e., treatment does not have any effect on the covariate in this scenario). In Scenario 2, the treatment effect on the outcome is entirely mediated by the time-varying covariate: drug treatment lowers the covariate, which lowers the outcome value. Scenario 3 represents a hybrid scenario with 50% direct and 50% mediated treatment effects. The total treatment effect for a specific drug concentration was kept the same across the scenarios. In each scenario, linear relationships were used for exposure–response (E-R) of the treatment effect(s) on outcome and covariate, and for the effect of the covariate on outcome.

Two E-R models were fitted to each dataset: Model A used the baseline value of the covariate as a predictor of the outcome, while Model B used the time-varying covariate instead. The bias was estimated by comparing the estimated E-R slope with the true total exposure effect on the outcome (Figure 2b).

In Scenario 1, where the time-varying covariate is not an intermediate for the treatment effect, both models yield unbiased E-R slope estimates on average. However, model B exhibits a 2.2-fold reduction in the average standard error of the E-R slope estimate compared with Model A and smaller variance in the estimated bias (Figure 2b). In Scenarios 2 and 3, Model A produces unbiased E–R slope estimates (average bias of 0.4% and − 0.3%, respectively), while Model B, incorporating the time-varying covariate, results in biased treatment effects (average bias of −100.1% and −50.0%, respectively).

Scenario 1 highlights the potential benefits of incorporating time-varying covariates,^3-5 not only for enhancing statistical fit but also for improving the statistical power to estimate treatment effects, reflected in the reduced standard error of the estimated E–R slope. In Scenario 1, the covariate is not affected by the treatment (so there is no overadjustment bias) but is strongly predictive for the outcome. But when the time-varying covariate acts as an intermediate on the causal pathway between drug treatment and outcome (Scenarios 2 and 3), its inclusion in the model can lead to substantial underestimation of the treatment effect, reaching up to 100%.

Interestingly, the difference we see in the estimated slope between model A and model B is directly associated with the percentage of the total effect that is mediated by the time-varying covariate: the mediation of 0%, 50%, and 100% in Scenarios 1–3 gave an average estimated bias in model B of ~ 0, −50% and −100%, respectively. This phenomenon is leveraged in epidemiology in so called mediation analyses, to estimate the level of mediation of the treatment effect by an intermediate, and these techniques can also be applied within the context of PKPD modeling.

In this perspective, we will briefly discuss the traditional approach to mediation analysis (i.e., the difference method). In mediation analysis, the total treatment effect represents the sum of the direct effect (the effect not mediated by the intermediate) and the indirect effect (the effect mediated via the intermediate).⁶ For the difference method, two models are fitted to the data. One regression model includes the exposure and baseline covariates to predict the outcome (model A in previous example), while the second model also includes the potential postbaseline mediator as a predictor of the outcome (model B in previous example). In the model without the potential mediator, the estimated effect of exposure is interpreted as the total effect on the outcome. In the model with the potential mediator, the estimated effect of exposure represents the direct effect (i.e., the effect remaining after controlling for the mediator, which “explains away” the indirect part of the total effect of exposure on outcome). The indirect effect is obtained by subtracting the direct effect from the total effect. The level of mediation is than calculated as indirect effect/total effect.⁷

The traditional mediation analysis method described above, remains the most commonly used method.⁸ However, it relies on several assumptions, such as the lack of interaction effects or uncontrolled confounders.⁷ If one or more of these assumptions are violated, the traditional mediation analysis may provide flawed results. More advanced mediation analysis techniques do exist, such as those that employ counterfactual techniques from causal inference. For readers interested in such techniques, we refer them to more in-depth articles.^{6, 7}

When considering the use of post-baseline covariates (time-constant or time-varying) in PK/PD models, it is important to evaluate the risk for overadjustment bias. For example, there will likely be overadjustment bias when the treatment has a non-negligible effect on a time-varying covariate that is also included in the model. This evaluation should ideally be based on both exploratory data analysis (example shown in Figure S1) and expert opinion. One could also perform a mediation analysis to assess the overadjustment bias, as the estimated level of mediation is equal to the estimated overadjustment bias of the total treatment effect when the time-varying covariate is included in the model.

When the risk of overadjustment bias is high, one might restrict the covariate analysis to baseline values to avoid underestimation of the total treatment effect on the outcome. This will in general be the simplest solution. However, a limitation of this approach is that it doesn't characterize the underlying relationships of treatment, intermediate and outcome. If disentangling these relationships is of interest, one might consider more elaborate sequential or joint PK/PD modeling approaches in which the intermediate variable is modeled as an additional outcome variable, rather than as a covariate, which would also remove overadjustment bias.

We also discussed the use of mediation analysis to quantify the level in which an intermediate variable may mediate the total effect of a treatment on an outcome. High levels of mediation—while not sufficient to demonstrate causality—may support the use of the intermediate as a surrogate endpoint or bridging biomarker.⁹ This has particular relevance in drug development, especially when clinical end points are rare or necessitate prolonged follow-up. An example is a PK/PD analysis of finerenone that showed the drug effect of slowing down kidney function decline (late outcome, requiring long follow-up) could be characterized via its effect on albuminuria (early outcome), supporting albuminuria as a potential surrogate end point for future studies.¹⁰ Additionally, using mediation analysis to support biomarker development may also contribute to personalized treatment.

In conclusion, the use of postbaseline or time-varying covariates can lead to overadjustment bias of the estimated treatment effect. When the risk of overadjustment bias from time-varying covariates is high, it would be better to use baseline covariates or joint modeling approaches instead. Additionally, the mediation analysis may be applied within the context of PK/PD modeling to support biomarker development.

SCG received a personal grant from the Van Wersch Springboard fund.

The authors declared no competing interests for this work.