V2 ACHER: Visualization of complex trial data in pharmacometric analyses with covariates.

Pharmacometric models can enhance clinical decision making, with covariates exposing potential contributions to variability of subpopulation characteristics, for example, demographics or disease status. Intuitive visualization of models with multiple covariates is needed because sparsity of data in visualizations trellised by covariate values can raise concerns about the credibility of the underlying model. V2 ACHER, introduced here, is a stepwise transformation of data that can be applied to a variety of static (non-ordinary-differential-equation-based) pharmacometric analyses. This work uses four examples of increasing complexity to show how the transformation elucidates the relationship between observations and model results and how it can also be used in visual predictive checks to confirm the quality of a model. V2 ACHER facilitates consistent, intuitive, single-plot visualization of a multicovariate model with a complex data set, thereby enabling easier model communication for modelers and for cross-functional development teams and facilitating confident use in support of decisions.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Efficient communication of pharmacometric models through visualization is essential for drug development decision making, and trellis and visual predictive check (VPC) plots are standard tools.

WHAT QUESTION DID THIS STUDY ADDRESS?

It presents an alternative visualization method to help understand, evaluate, interpret, and communicate pharmacometric models.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

V2ACHER facilitates model communication by enhancing visualization of pharmacometric analyses with covariate effects, including population pharmacokinetics, population pharmacokinetics–pharmacodynamics, and model‐based meta‐analyses models. The method accounts for covariate effects on the dependent and independent variables (and, optionally, for random effects, a.k.a. between‐subject variability), retaining the response curve of a selected reference population. Applying V2ACHER to VPCs (V3PC) also simplifies model qualification, enhancing the properties of prediction‐correction VPCs, and potentially facilitating the communication of modeling results to nonmodelers.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

V2ACHER and V3PC are intuitively interpretable visualizations able to reveal model strengths and weaknesses. They can help engage discovery and development teams in modeling and, thus, increase the likelihood of appropriate model impact on decisions.

INTRODUCTION

Pharmacometric (PMx) models, including pharmacokinetics (PK), pharmacodynamics (PD), and exposure–response models, are used to quantify and interpret observations, support optimal experimental nonclinical and clinical design, and, thus, to accelerate drug discovery and development. There is increasing endorsement of PMx modeling from regulatory agencies (US Food and Drug Administration, European Medicines Agency, International Conference on Harmonisation), especially with the publication of the model‐informed drug discovery and development white paper., , Despite widespread regulatory support of, and even requirement for, model‐informed drug discovery and development, the impact of models on decisions can be reduced if nonmodeling team members, key stakeholders, or regulators have less experience in interpreting and evaluating the properties of PMx models.

Nonlinear mixed‐effects (nlme) PK‐PD models have been used for more than 40 years to describe the relationship between dosing parameters and physiological outcomes in a clinical population that includes multiple subgroups as well as across multiple clinical trials and nonclinical experiments., This approach differs from standard group statistical analysis in that the data from all individuals are considered together, but the effects of individual characteristics, such as age, sex, and disease status, can be used to account for between‐subject variability (BSV). Incorporation of covariates into the model reduces unexplained variability in the results by reducing the confounding effects of these variables and permitting meaningful evaluation of efficacy and safety from sparser data sets on fewer subjects. These model results can also be used to inform dose selection and dose modification during the drug discovery and development process, thereby minimizing the costs of testing ineffective dosing regimens. Modeling with covariates is essential to consider when there are different types of observations, observation sites, formulations, and (for translational modeling from nonclinical data) species. Clinically, it is especially important for understanding the impact of drugs on subpopulations, such as children or patients with renal or hepatic impairment, whose responses are typically represented by small populations, and who might require adjustment to dosing regimens based on difference in exposure to standard populations.

One challenge in communicating results of population‐based nlme PMx models or of model‐based meta‐analyses (MBMAs) is a lack of visualization tools that account for the impact of covariates on BSV and allow the original data (often sparse) to be shown in the context of the model (often with complex structure, random effects, and covariates). The objective of this work is to share a PMx tool/method that (1) fills this need by providing intuitive visualization of variability‐aligned, covariate‐harmonized effects with a reference (V2ACHER), (2) results in easier and even more productive communication with nonmodelers about model properties and results, and (3) can also substantially increase the interpretability of visual predictive checks (VPCs; often used to ensure the quality of a model and its covariate structure).

An overview of V2ACHER will be followed by a description of the algorithm details. Three examples with increasing complexity will show how the transformation steps work together to achieve the desired visualization in the presence of random effects (BSV) and/or one or more covariates. A fourth example will illustrate how using V2ACHER‐transformed data to create a VPC (V3PC, for V2ACHER‐VPC) can substantially simplify the interpretation of a VPC and enhance the ability to detect model misspecification.

METHODS

Overview of the V2ACHER data transformation algorithm

The V2ACHER method transforms the raw data (original observations) in several steps, resulting in more effective visualization of the data and the model in a single plot. The method assumes that the data have already been modeled using nlme, generalized nonlinear least squares (gnls), or related methods. The goal of the transformation steps is to match the difference between a transformed data point and the reference curve to the distance between that point and the model curve corresponding to that point's covariate values. Here, “difference” typically means vertical distance on the visualization appropriate to the model structure: for log‐normal variability, for example, the difference preserved would likely be on a logarithmic axis and, thus, correspond to preserving the ratio between the value and the (covariate‐adjusted) prediction (which is the same as preserving arithmetic difference between the logarithms of those two values).

At a summary level, V2ACHER creates a consistent, intuitive, single plot showing the model predictions (including confidence intervals) along with the raw data that underlie the model. The transformation of the data relative to the reference curve adjusts for covariate effects (CEs) for the independent and/or dependent variables, depending on the model. For the independent variable, the units are (if required by CEs) effectively matched to those of the reference data by a multiplicative factor (or the equivalent shift in logarithmic units). In a PK‐PD Emax (sigmoidal) model, for example, if male subjects have fourfold higher EC50 values than females (and no other CEs), this step scales male concentrations by ¼ so that the male and female curves (and, thus, data fit well by the model) would coincide since concentrations would now be consistent relative to the reference EC50 value (i.e., in EC50 units). (EC50 is the concentration producing 50% of the maximal effect.)

The dependent variable data values may also be scaled (linearly) to account for the CEs on them. Scaling in this step can be illustrated by a compound with linear PK and with bioavailability of 2f in males and only f in females. To properly visualize concentration data on the same plot using the female dose‐concentration model as a reference curve, the data for male subjects should be multiplied by ½. Confidence intervals (if any) associated with the observations are scaled so that they visually represent, in an appropriate manner, the variation in the data relative to the position of the reference curve. (Typically the end points of the confidence interval would be transformed in the same way as the measurement, but the intervals could be unchanged or a different technique used in cases where that was considered to yield an overly stringent or flexible impression for the purpose at hand.)

Details of the V2ACHER data transformation algorithm

The four steps of the method are shown in Figure 1.

FIGURE 1

Schematic of the visualization of variability‐aligned, covariate‐harmonized effects with a reference (V2ACHER) method. Symbols are defined in the main text. Illustration plots are of response (R) versus independent variable (IDV). Both axes are here assumed to be on a logarithmic scale for consistency with the scaling used to illustrate Steps 3 and 4. The purple points correspond to data transformed using all steps through the current step (the number of which appears above the point), the grayed‐out points correspond to data from previous steps. The red and solid blue curves are the reference and covariate curves defined in Step 3, respectively. The dashed purple curve has been transformed by Step 3, and the dashed purple curve on top of the red curve by all steps. Gray arrows indicate the general direction of transformation in each step. Step 1 is model creation, and no transformation is performed. In Step 2 (which is optional), only data are transformed, and the shift will generally be toward the typical value R. Step 3 transforms both data and the corresponding covariate model curve along the IDV axis. Step 4 transforms them on the R axis so that the vertical distance from point 4 to the reference curve is the same as the distance from point 1 to the solid blue curve. When the blue curve has been transformed using all steps, it will coincide with the red curve, as shown. BSV, between‐subject variability

In Step 1, a PMx model with covariate(s) and/or random effect(s) is identified: this model can be described using gnls, nlme, empirical, or mechanistic models.

In Step 2, if the data set is derived from individuals or multiple groups, such as separate clinical trials, Equation (1) is used to shift the data to account for BSV (with trials treated as subjects in MBMA):

Rcov represents the response value (dependent physiological or clinical value, the datum ordinate) and R ipred and R pred represent (respectively) the corresponding individual and typical responses predicted by the model. This step removes the BSV and transforms the observations to the values they would have if the data could be completely explained using the typical prediction and the residual error. Depending on the intended use of the visualization, this step might be omitted from the scaling process (see Examples 1 and 2): as explained in Example 4, this step might need to be omitted when V2ACHER is used in the context of model qualification using V3PC.

Step 3 accounts for CEs impacting the independent variable (IDV); for example, dose for a PK model, or exposure for a PK‐PD model. This could include attributes such as types of observations, observations sites, formulations, and (for translational modeling from nonclinical data) species. The CEs are relative to those for a “reference curve,” the model prediction for typical response in the “reference population,” where a “population” is effectively a particular combination of covariate values such as low birth weight, male, born at full term. (By analogy, the “covariate,” or “nonreference” curve refers to one for nonreference covariate values.) The reference population is defined as the population to which other populations’ model curves will be aligned (by Steps 3 and 4). The selection of this population is completely arbitrary (as long as its response is predicted by the model) and is typically based on the interest of the modelers and/or development team, for example, the target population for a compound that also has response data in other populations. Equations ((2), (3), (4), (5a), (5b), (6a), (6b), (7), (7a), (7b), (7c), (8), (8a), (8b), (9)) (next) show how the IDV values should be transformed for the nonreference population and will be used both for its model prediction (covariate curve) and its associated data points.

The effect of a covariate (cov) on a model parameter is here assumed (without loss of generality) to have the form in Equation (2), where θ cov is the covariate‐adjusted value for the model parameter θ. The CE can, for example, be defined as a fraction using Equation ((3), (5a)) for a categorical covariate and Equation ((3), (5a)) for a continuous covariate: where cov is the covariate value and θ cov is the model CE impacting the model parameter θ (such as EC50). It is assumed that θ, and θ cov, , were determined during model building. Typically the definition of the “reference population” (defined in the text between Equations (1) and (2)) and the parameterization would be arranged so that that CE =1 for the reference curve (i.e., there are no CEs because the effects are defined relative to the model curve for the reference population).

CE is then used to adjust the value of the independent variable, IDV, transforming it to IDVtr to account for CEs using Equation ((3), (5a)),

The response curve, R nonref(x), of typical predicted values for the nonreference population at the transformed IDV = x, is now plotted at the abscissa value IDVtrans(x) (instead of at x) creating the curve (x, R step3(x)) so that R step3 is defined using Equation ((4), (5b)),

This means that, if the original point was plotted at (x, y) = (IDVcov, R BSV), that point is now plotted (as shown in Figure 1) at (IDVtr, R BSV), where IDVtr is the transformed value of IDV corresponding to the original (untransformed) value IDVcov. Analogously, before Step 3 the nonreference curve is plotted as R(x) versus x with R(x) = R nonref(x) (x is the value on the IDV axis). In Step 3, the curve goes through the same transformation as an observation (IDVcov, R BSV), and so after Step 3 the curve would be plotted as R(x) versus x with R(x) = R step3(x) (as shown in Figure 1). This means that the value of the typical prediction curve at IDV = x is now used (plotted) at the value of x transformed to the equivalent reference coordinate IDVtrans(x). So if the EC50 of my reference curve is 1, and that of my covariate curve is 2, the values of the covariate curve at x = 1, 2, and 3 would be used to form the transformed curve at x = 0.5, 1, 1.5 (respectively), as illustrated by the left shift in the dashed purple curve relative to the blue curve in Figure 1.

Step 3 can be repeated multiple times (once per covariate) if covariates act independently on the model parameters. (Generalization to interacting covariates can often be achieved by reparameterization of the model in terms of independent covariates.) The result of Step 3 is that observations are aligned along the IDV axis in a way that results in the ability to more easily interpret the IDV independent of any CE. The reason for this can be illustrated by considering a dose–response relationship that depends on the population, such as when the pediatric dose needs to be half of that for adults to obtain the same response, CE = 0.5 (e.g., the EC50 for a typical child is half that of an adult). Thus, scaling the dose (IDV) from pediatric to adults (i.e., using adults as the reference population) would use IDVtr =IDV/0.5. As the result, all pediatric doses “move” to equivalent adult doses. The covariate(s) can include discrete and/or continuous variables.

Covariates such as level of disease severity and demographics can, of course, also impact a response, for example, a sex effect on exposure–response in which females experience twice the effect as males at the same concentration. Step 4 will transform the typical predicted response curve for the nonreference population, R step3(IDV), so that its predicted (BSV‐corrected, IDV‐transformed) response curve, R tr(IDV), will (with R tr defined using the transformation defined in Equation (6a)) appear to be identical to that of the reference population. Furthermore, transformed versions of the observations preserve the (plotted vertical) distances from original observations to their respective prediction curves. This principle ensures that the transformed curve and data will transparently represent the degree of fit of data, conveying the same attributes of model fit (or lack thereof) as the original plots comparing observations to model curves. Stated more rigorously, the transformation must be chosen so that, when the observed (BSV‐corrected) data, R BSV, are transformed the same way as their (typical) prediction curve, the (plotted vertical) distance from point to prediction (curve) is the same as when the original (untransformed) point is compared with its corresponding curve in the original (untransformed) plot (cf. Figure 1).

To obey these principles, Step 4 plots data and predictions as R (ordinate) versus IDV (abscissa) after transformation. Step 4 scales an observed data point that is at (IDVtr, R BSV) after Step 3 so that it is now plotted at (IDVtr, R tr), with IDVtr defined in Equation (5a) and the transformed value of the observation, R tr, given by Equation (6a) (next). Step 4 scales the response curve in the new visualization so it is plotted as R trans(x) versus x (instead of R step3(x) versus x) with the transformed curve, R trans(x), given in Equation (6b): with the transformed value, IDVtr, of IDV and definition of R step  given in Equation (5a,b). The cancellation of terms in Equation (6b) (resulting in ) is because of the definition of in Equation ((4), (5b)) and shows that the transformation of points and curves has the desired property of relating everything to the reference curve: if all the observed points were exactly on their respective predicted covariate curves, all the transformed points would be on the reference curve. For an additive error model, multiplication and division in Equation (6a,b) are replaced by addition and subtraction, respectively. For a mixed error model, both shifting and scaling are needed in proportion to the coefficients of the respective terms (in this case the preservation of plotted vertical distance could be approximate).

Software

For Examples 1 and 3, R version 4.0.2 (2020–06–22) was used for all modeling, simulation, data processing, and plotting. For Examples 2 and 4, NONMEM (version 7.4.3) was used for modeling and simulation using PsN version 4.8.1, and R (version 4.0.2 [2020-06-22]) was used for data processing and plotting. The Certara tidyvpc package (version 1.1.0 [2020-09-28]) provided VPC percentiles and prediction intervals for Example 4.

Research Involving Human Participants

The studies involving human participants were approved by the appropriate institutional and/or national research ethics committee and were performed in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards.

RESULTS

Example 1: model with single discrete covariate and without BSV

To evaluate the methodology with relatively little uncertainty and noise, a simulation data set was taken from a published example where a different method was used to visualize the data and for which only Step 4 of Figure 1 is needed. The simulation data provided exposure–response in which the males responded differently than the females. The V2ACHER method was then applied to this example to adjust for sex effects on the dependent variable. There were two simulated arms in the study: one with 500 males and the other with 500 females. A gnls model was fit to the simulated data (Step 1) using the gnls function in the R package “nlme,” accounting for the sex effect. To enable relatively direct comparison with results of the previous method, the male and female data sets were each divided into seven binned exposure quantiles representing 7 different “dose” levels, and the mean responses and 95% confidence intervals were calculated. Since there was no BSV and the covariate affected the dependent variable only, neither Step 2 (BSV correction) nor Step 3 (transformation of exposure) were required, leaving only Step 4 (response transformation). The predicted values for a typical male individual were selected as a reference in this example. Applying Step 4 (Figure 1) scaled the female response data to the male reference curve using Equation (6a).

Figure 2 compares four approaches to visualizing this simulated exposure–response data set, illustrating the impact of scaling the response for sex as the population covariate. Figure 2a shows the original data. Figure 2b shows how each sex grouping fits its respective (covariate‐adjusted) model; it does not allow immediate, intuitive understanding of the relationship between the two sets of data. The V2ACHER‐transformed data (Figure 2c) reveal the relationship between the data set and a single, coherent model. Although the difference between Figure 2b,c is likely to be less impactful in a relatively simple example with parallel model curves (as seen here), it is often an impediment when there are more (combinations of) covariate values and in the presence of covariate‐driven changes to curve shapes. If the data are combined using the method originally proposed (Figure 2d), the predominance of lower exposure values among the males biases the apparent model fit, potentially complicating assessment of the model's goodness of fit to the data.

FIGURE 2

Visualization of raw and transformed exposure–response data for Example 1. Data are from a simulation of the model in Overgaard et al., showing the exposure (area under the curve [AUC]) and response for 500 male and 500 female subjects binned into seven exposure (AUC) quantile groups separately for each sex. Shown are the mean ± 95% confidence interval for the response in each group. The red symbols and curves are for simulated males, the blue symbols and curves are for simulated females in (a), (b), and (c), and the axes scales are the same for all four panels (a)‐(d). (a) Unclear exposure–response relationship in the data set with sex information. The data show high apparent (sex‐effect confounded) variability prior to adjustment for the effects of sex. (b) Superimposition of data and the separate covariate‐adjusted prediction curves for males and females. (c) Data from panel (a) are shown after scaling to adjust for the sex effect (using the male curve as the reference) elucidating how data from both sexes are fit by a single, coherent model. (d) Visualization of raw data and model curve as proposed in the original reference. The dark blue points show the medians of the “observed” responses (and 95% confidence intervals) from their respective AUC bins, where the dark blue curve is generated similarly using the predicted responses. The light blue curve is the mean of the male and female model curves. The data used are from a pseudo‐random number generator, and therefore the figure is not expected to be identical to that originally published

Example 2: model with two discrete covariates

This example is intended to illustrate the impact of either including or excluding correction for BSV (as described in Figure 1) and to demonstrate the use of V2ACHER when there is more than one covariate, and the covariates were assumed to impact only the dependent variable. (The example does not include any IDV scaling.) A simulation analysis was constructed based on a published PK analysis of neonatal phenobarbital clearance with multiple demographic covariates (additional details in Supplementary Material S5). Birth weight and age were found to impact clearance, and weight the volume of distribution. The simulation parameters were the following: body weight at birth (BW) in kg = 0.5, 1.5, and 3.5; postnatal age in days at first dose administration (AGE) = 1, 11, and 21; single intravenous dose (DOSE) = 4 to 40 mg/kg; and infusion rate (INF) = 36 to 72 mg/h. The weight at first dose administration was highly correlated with age and birth weight, thus the model used Equation (7) to simulate the weight:where “noise” represents multiplicative variability distributed normally with a mean of 1 and a standard deviation (SD) of 0.1.

The CEs were modeled as follows:

Where: where CLx is the effect of covariates on clearance for covariate x = AGE, BW, or (for COV) both.

A sparse, randomly sampled data set containing 50 samples from 20 neonates was generated based on the published model simulation. Model parameters for the two‐compartment model with first‐order absorption and elimination, including BSV on clearance and central volume of distribution, were then re‐estimated using the smaller data set (original and re‐estimated parameter values are in Supplementary Material S1). The 50 data points were then adjusted either:

Without correction for BSV, scaling data points with Step 4 of Figure 1 using Equation (6) to account for the effects of AGE and BW using the reference curve for the typical individual with BW = 1.5 and AGE = 11 days with DOSE = 38.8 mg and INF = 50 mg/h or

With BSV correction prior to the steps performed in (a), applying Equation (1) to account for the effects of BSV.

Transformation of exposure as in Step 3 of the scaling process was not required in either of these since covariates affected only the dependent variable.

Figure 3a shows the simulated data, with high variability and a trend of decreasing maximum concentration with time after dose. Figure 3b shows the simulated data stratified by age and birth weight; the sparsity of data and variability are evident when trellised by subgroup, and, although there is no clear relationship between concentration, group, and time, the trend of decreasing maximum can be seen in some groups. In Figure 3c, after V2ACHER transformation with scaling for AGE and BW CEs, the random distribution of the transformed data around the curves reveal that they are consistent with the covariate‐dependent PK time profile. Figure 3d shows how accounting for effects of BSV in clearance and central volume of distribution in the scaling on top of the steps taken in Figure 3c shifts data points further toward the red (re‐estimated) and blue (simulated) model curves, which helps to identify the extent of remaining residual variability in the model. Also, although it is not clear from Figure 3a,b that the simulated data are sufficient to support model re‐estimation, the proximity (not dependent on V2ACHER) of re‐estimated and simulated curves in Figure 3c show that the data are sufficient to re‐estimate the model consistently.

FIGURE 3

Visualization of raw and transformed pharmacokinetic data for Example 2. Raw data from a simulation of the phenobarbital concentration over time after a single intravenous dose in neonates are shown (a) combined into a single graph without adjusting for covariate effects, (b) trellised by AGE (1, 11, and 21 postnatal days) and body weight at birth (BW; 0.5, 1.5, and 3.5 kg), (c) after scaling for covariate effects where gray symbols represent the original data and colored data points the data points scaled (using V2ACHER) to the reference typical individual (BW 1.5 kg and AGE 11 days), and (d) after both adjustment for between‐subject variability (BSV) and scaling for covariate effects with gray and colored symbols as described in (c). No independent variable scaling was performed in any of the plots, that is, data points were only moved in vertical direction. The blue curve shows the predicted phenobarbital pharmacokinetics for a typical reference individual (see the text) based on the original model and the red curve (near or over the blue curve) shows the predicted phenobarbital pharmacokinetics for a typical reference individual after re‐estimation of the model parameter using the smaller, simulated data set

Example 3: model with three discrete covariates

An analysis data set was assembled from 17 clinical trials containing three key covariates (meta‐data summarized in Table 1). An MBMA model was created to describe the disease incidence rates as a function of PD biomarker values, as shown in Equation (8).

TABLE 1

Properties of the 54 summary‐level published data points from 17 clinical trials by population, disease severity, and compound formulation showing the sparsity and diversity of the data

	Population A		Population B		Population C		All Populations		Total
	CF1	CF2	CF1	CF2	CF1	CF2	CF1	CF2
Severity 1	5a	1	2	1	13	1	20	3	23
Severity 2	0	0	8	2	6	1	14	3	17
Severity 3	0	0	4	0	9	1	13	1	14
Sum	5	1	14	3	28	3	47	7	54
Total	6		17		31		54

Abbreviation: CF, compound formulation.

Each data point represents an aggregated incidence rate of one arm of a clinical trial.

IR (the response) represents an incidence rate, IRmax represents the maximum IR, IRmin represents the minimum possible IR, C represents the PD biomarker value of titer, suppr represents the inhibition depending on C, EC50 represents the C at which 50% inhibition is achieved; and γ determines the maximum slope of the sigmoidal curve.

CEs were modeled as follows: where EC50x is the effect of covariates on clearance for covariate x = FORMULATION, POPULATION, or (for COV) both. IRmaxy is the CE for y=POPULATION, DISEASE (disease severity 1, 2, or 3), or a combination thereof (for COV). Random effects for each trial were modeled to impact IRmax proportionally with IRmaxtrial = IRmax. , where IRmax is the typical trial maximum response and reflects the random effect for trial i (between‐trial variability). We can arbitrarily choose IRmaxx = 1 for one population and for one disease severity level.

Because the data depend on three covariates and on random effects as well as C, understanding and communicating the model and data required the creation and use of V2ACHER.

In this example, there is a mechanistic reason to expect a proportional variability between trials, and the key results of interest involve the relative (within‐trial) effect of active versus placebo arms. Thus, the raw data were adjusted for between‐trial variability using Equation (1) and for differences in formulation using Equation (5a,b) of Step 3 (this adjustment has no impact on the key results of interest). IR was then scaled to account for population attributes and disease severity using Equation (6a,b). The results, shown in Figure 4, were plotted on a semilogarithmic scale to improve visualization of proportional differences, especially at low IR.

FIGURE 4

Visualization of raw and transformed biomarker–response data for Example 3. Raw data for the level of a pharmacodynamic biomarker value and the disease incidence rate in humans are shown (a) combined into a single graph without adjusting for covariate effects. Each data point represents a clinical trial arm with error bars for their 95% confidence interval (CI). The symbol indicates the disease severity level, and the symbol's size indicates the size of the trial arm. (b) Trellised by population and disease severity level. (c) The final output is shown superimposed on the predicted, typical biomarker‐incidence rate relationship for reference Population B at reference disease severity Level 2, with an associated 95% CI

Figure 4a shows the sparsity of the data and challenge in interpretation without consideration of CEs; the data do not show an obvious trend, and variability is high. Figure 4b shows that trellising the IR versus C data by population and disease severity yields sparse plots that obfuscate any relationship. Figure 4c shows V2ACHER‐transformed data overlaid on model prediction, demonstrating that the transformation helps to more clearly reveal the relationship IR(C; covariates).

Example 4: V3PC—creation of a VPC using V2ACHER‐transformed data

The method presented here is related to the prediction‐corrected VPC (pcVPC) method described by Bergstrand et al. The VPC method was developed to allow modelers to visually assess model quality using a simulation approach. In pcVPC, the value of the resulting VPC may be limited by binning across important covariates. This fourth set of examples illustrates how transforming the data first using V2ACHER before creating a VPC plot (“V3PC” for V2ACHER VPC) can lessen this limitation.

Since VPCs are used to help assess attributes of BSV, Step 2 of the scaling, which shifts the data to account for BSV, should be omitted so that any biases introduced by or into estimates of BSV will not be hidden by the visualization.

The use and relative merits of V3PC and pcVPC can be illustrated using a simple Emax model:where Rmin, Rmax, and R correspond to minimum, maximum, and simulated response, respectively; C is the concentration of the drug; and EC50 is the concentration leading to half‐maximum increase of the response. The CE was modeled using Equation (3) with EC50COV = 0 for Population 1 and EC50COV = Θpopulation‐2 for Population 2.

A rich (Data Set 1) and a sparse data set (Data Set 2) were created; both simulations include two populations to allow for simulations with two different EC50 values. A multiplicative (relative) CE was used on the EC50 for the second population (Equation (2)). For simplicity, this was the only population difference in the simulation so that no response scaling (Step 4) is needed. (The example could be easily extended to include response scaling.)

Data Set 1 provided 120 observations from 24 subjects, Data Set 2 included 30 observations from 10 subjects. Responses were simulated in NONMEM based on Equation (9) using Rmin and Rmax values of 1 and 10. EC50 of Population 1 (reference population) was 0.5, and the EC50 of the second population was 3. A BSV term shared between Rmin and Rmax was included, assuming normal distribution with mean = 1 and SD = 0.2. In addition, random additive (SD = 0.6) and proportional (SD = 0.2) noise were added. Parameters were re‐estimated for each simulated data set using NONMEM. The resulting parameter estimates are listed in Supplementary Material S2.

A total of 500 replicates of each model were simulated for VPC using NONMEM. VPC analysis was then carried out using tidyvpc in both normal and prediction‐correction (pc) modes. V3PC plots were obtained by applying Steps 1 and 3 of the V2ACHER transformation to both the original and simulated data sets, with Population 1 as the reference population, followed by a normal VPC analysis on the V2ACHER‐transformed data. The analysis included automatic binning into five bins with each containing an equal number of observations and calculation of the median and respective quantiles as well as 95% confidence intervals, just as was done for the VPC and pcVPC.

A comparison between VPC, stratified pcVPC, pcVPC, and V3PC plots for the rich and sparse data sets is provided in Figure 5. As expected, the median of simulated data in normal VPC and pcVPC is an average of the two original simulation curves, that is, it does not reflect the properties of either of the simulated populations. In contrast, V3PC represents the properties of both populations (using representation based on Population 1 as a reference). All methods perform well on the rich data set if it is not stratified (Figures 5a,e,f), with the confidence intervals somewhat more separated in the V3PC than in pcVPC and the relatively symmetric scatter of the points relative to model prediction preserved in V3PC and not in the VPC or pcVPC. The larger (better) separation of confidence intervals in V3PC could be because the method does not bin across covariates and might thus better represent the true uncertainty in the models. Stratifying by population in the pcVPC (Figure 5c) results in very little prediction correction and a plot similar to stratified VPC (not shown). The stratification leads to larger uncertainty in the predictions, which is reflected in wider, overlapping confidence intervals even with the rich data set. For the sparse data set, the appropriateness of the model is harder to assess with VPC and pcVPC: due to the binning across covariates, profiles seem erratic, with a pattern not representing the expected shape from a simple Emax model (Figures 5b,f). The stratified pcVPC (Figure 5d), however, helps to identify which population contributes most to the estimation of the lower plateau of the curve, as red data end before the plateau is reached. In contrast, with sparse data V3PC reveals its strength by giving an integrated view of the results: it retains (Figure 5f) the expected shape of response and the relationship of the data to the underlying model.

FIGURE 5

Comparison of prediction‐corrected visual predictive check (pcVPC) and visualization resulting from applying V2ACHER to the data used in a visual predictive check (resulting in the “V3PC”) for rich and sparse data. Plots use data from Example 4, and black dotted and dashed lines correspond to the median of observations and empirical 90% observation interval, respectively; gray solid lines represent the median of simulated data and 90% prediction interval. Gray‐shaded areas show 95% confidence intervals around respective simulated quantiles: values for each bin are plotted at the center of the bin and then plotted as interpolated between bin centers with a piecewise linear function (assumed to have a constant value over the lowest and highest one‐half bin width). The prediction curves for Populations 1 and 2 from the re‐estimated models are included as red and blue solid lines, respectively, to illustrate where the responses would typically be expected. These would normally not be included in traditional visual predictive check (VPC) or pcVPC plots and are not intended to be used to assess model appropriateness in this example. Data points are color coded by population. (a, b) For VPC plots of rich and sparse data, red and blue points represent the original observed data from Populations 1 and 2, respectively. (c–f) For stratified (c, d) and unstratified (e, f) pcVPC plots of rich and sparse data, light red points represent the observed data from Population 1, and light blue points represent the observed data from Population 2 after prediction correction. (The red and blue solid lines are not prediction corrected.) (g, h) For V3PC plots of rich and sparse data, red points represent the observed data from reference Population 1. Purple points indicate data from Population 2 after V2ACHER transformation (adjustment for the covariate effect with Step 3); they are positioned relative to the (reference) red curve as though they had been obtained from Population 1

DISCUSSION

V2ACHER is a PMx scaling method that adjusts for CEs on the dependent and independent variables and, optionally, accounts for random effects, thereby facilitating interpretation and assessment of a population PMx model using a single graph. The four examples demonstrate that the technique is flexible and can be applied to multiple types of PMx data and models, including exposure response, PK‐PD, biomarker response, MBMA, and other models. The technique can be applied to data sets with a simple covariate structure as well as complex population data sets derived from multiple trials. The resulting plot can help communicate modeling results to team members who might be most comfortable interpreting the traditional pooled statistical analysis of simple data sets. The main purpose of this communication would be to highlight population commonalities, that is, how the different populations with different properties can be described with a single underlying model. Other visualization tools might be more appropriate for communicating covariate impact resulting in differences between population responses. In addition, the version of V2ACHER shown here cannot be used on nonstatic (e.g., compartmental) models, although this is the subject of future work.

Visualization of modeling results is critical to understanding, evaluating, interpreting, and communicating them. It is often also critical for engaging a preclinical or clinical development team in model evaluation, for obtaining their feedback to use in model development, and, ultimately, for engaging the model to support decision making. By applying the V2ACHER stepwise transformation method, complex, dense, or sparse data sets with BSV and multiple covariate dependencies can be presented in a single plot and overlaid intuitively on a reference curve. This visualization enhances insight into whether and/or how the data might be integrated into a coherent model. Furthermore, the use of V2ACHER to transform data to a reference in the context of performing a VPC results in an intuitively interpretable visualization (V3PC) for evaluating (detecting) potential model misspecification.

V2ACHER may improve on the visualizations proposed by Overgaard et al., who showed the discrepancy between the simplistic summary of the combined data and the model prediction and proposed a visualization to mitigate this discrepancy. Example 1 showed that V2ACHER‐transformed data match the model prediction and can provide a more easily interpreted visualization that demonstrates that data from females can be described with the same model structure as the one used for males.

PMx models can be particularly useful in the analysis of sparse data. Drug development guidelines require that medications be studied in children, and these studies typically produce sparse data. As a result, V2ACHER could be especially useful for pediatric drug development and dose prediction. Example 2 demonstrates an example of this, showing that—after scaling for covariates included in the drug clearance model—V2ACHER can be applied to a sparse neonatal PK data set. The resulting plots (Figure 3) show intuitively that the underlying structural model is identical for all of the subgroups included in the parameter estimation and illustrate how the CEs in the model allow prediction of the response with that structural model. The visualization thus demonstrated that all neonatal subgroups share common properties that can be described with that model: other visualization methods might highlight differences between the populations. Figure 3 also illustrates that V2ACHER can visualize how a complex model can be reliably estimated from a relatively sparse data set. These aspects of V2ACHER can be used to increase the trust of development teams in using such models for further decision making.

MBMA models can also rely on relatively small data sets, and, again, covariate‐based grouping of the data can result in sparsity that can be particularly challenging. The third example demonstrates how V2ACHER can mitigate (although not eliminate) the effects of that sparsity in an MBMA by producing an integrated visualization across the groups while maintaining the group‐based effects relative to the reference.

The fourth example demonstrates how applying V2ACHER to data before running a VPC (generating a “V3PC plot”) can be used in the context of model qualification to improve on some limitations in pcVPC arising when binning across covariates. The example demonstrates that V3PC is sometimes able to improve visualization of confidence intervals around predicted response quantiles (Figure 5g vs. 5e). For sparse data, it also removes (Figure 5h) erratic, potentially artifactual trends that are visible in both VPC and pcVPC profiles (Figure 5b,f). pcVPC might trigger questions on the interpretability of the shape of the response curve from development teams, whereas V3PC retains the response curve of a selected reference population. Retaining the response curve can be achieved for both VPC and pcVPC by stratifying by population (shown for pcVPC; Figure 5c). One advantage of this stratified view is that it helps to highlight which population contributes to the end of an estimated curve. However, due to the stratification, visualization of the response curve can come at the cost of overlapping confidence intervals, even in the rich data scenario, which, again, can limit the assessment of model appropriateness for the pcVPC method. Example 4 included a covariate on the independent variable. Differences between pcVPC and V3PC might be less apparent if the underlying model only included covariates affecting the dependent variable. In such a scenario, pcVPC would most likely enable modelers equally well to assess model appropriateness. However, additional value might be provided by V3PC, as the response curve would correspond to that of a chosen reference population.

Another important question is if the scaling applied in V3PC would hide model misspecifications. To address this question, a model known to be misspecified (simple linear regression) was fit to the data used in Example 4. As shown in Figure S3, both pcVPC and V3PC would reveal the model misspecification in the rich data scenario, as the central tendency of observed data is not fully contained in the confidence interval of median model prediction. This is less apparent for the sparse data example, as confidence intervals are very wide: V3PC would be expected to perform at least as well as pcVPC, although neither visualization can overcome a data set too sparse to determine misspecification. Although the main purpose of VPC is usually not to reproduce the shapes of the prediction curves, as an experienced modeler is well equipped to correctly interpret pcVPC plots, V3PC is easy to interpret and might therefore be better suited for direct communication with nonmodelers. The example illustrates how V3PC can improve the utility of VPC, especially for complex data sets. V3PC thus enhances the use of VPC for model assessment with similar principles to those of pcVPC, extending the interpretability to data with significant covariate considerations.

The utility of V2ACHER is not limited to applications demonstrated by the four examples. For example, disease progression may be characterized by various scales to quantify the severity, such as the American College of Rheumatology ACR20, ACR50, and ACR70 scores for rheumatoid arthritis. This is sometimes also combined (in a model) with treatment as a covariate altering that progression over time., The V2ACHER methodology is potentially suitable, by analogy with Example 3, to align (via the model) and visualize coherence of the measured and predicted scores. Discrete response data (specifically, ordered categorical data) described by a single model also require clear visualization to demonstrate consistency between model and data, and this could also be supported by V2ACHER.

Finally, although both V3PC and V2ACHER are demonstrated without stratification, if useful to do so, modelers can also decide to apply V2ACHER techniques for some covariates while stratifying by one or more others—combining the benefits of stratification with V2ACHER’s visual alignment across selected covariates.

The primary limitation of the V2ACHER method is that, as described here, it will not work on nonstatic (e.g., differential equation‐based) models for which time‐dependence shifts in a nonlinear manner between covariate values. The simplest example might be a compound for which the time of maximum concentration shifts with dose level and or other covariates. The scaling proposed for V2ACHER (as with other existing methods) will not align times of data points properly to compare their position relative to a reference time course. This is under investigation: a simplified version of dynamic time warping, is likely to enable generalization to nonstatic models.

It is also possible, under very restricted circumstances, for the visualization to mask potential problems with an underlying model, including artifacts created by missing or inaccurate data, or failure to include all the relevant covariates. This effect could emerge, for example, when data are combined across trials with substantially different designs so that the raw data include sources of bias difficult to evaluate in a single combined plot. To avoid overinterpretation of the data with a scaling procedure, it is always important to evaluate the underlying data and model, including mathematical and scientific assumptions in the chosen model and the statistical indicators of model goodness of fit. V2ACHER contributes to only one part of a decision‐making strategy; the underlying biological assumptions in the data, such as the proposed link between a surrogate biomarker and a clinical outcome, must be supported with information independent from the resulting visualization.

In addition, models including BSV can be affected by shrinkage. As Step 2 of the scaling method makes use of individual predictions for correction of data points, the extent of model parameter shrinkage should be assessed before applying this step. In models with high shrinkage, V2ACHER visualizations including Step 2 are expected to be influenced by this phenomenon in a similar way as, for example, plots of individual predictions versus observations, and the evaluation of model adequacy might be impaired. Using V2ACHER without Step 2 does not depend on individual predictions and is not influenced by parameter shrinkage. The specific application and shrinkage properties will determine if V2ACHER plots should be created with BSV correction, without BSV correction, or both ways.

As planning and analysis of clinical trials increasingly require consideration of MBMAs, more covariates, smaller data sets, and adaptive designs, it is increasingly important to enable analyses and reporting that support such work. V2ACHER’s advantages over previous techniques have already enabled it to support these kinds of tasks and decisions., , Future application is, thus, likely to enable improvement in the efficiency of future drug discovery and development.

CONFLICT OF INTEREST

J.R.S. is an employee of Merck Sharp & Dohme Corp., a subsidiary of Merck & Co., Inc., Kenilworth, NJ, USA and owns shares of Merck & Co., Inc., Kenilworth, NJ, USA. J.L., N.P., and S.Y.A.C. are employees of Certara Strategic Consulting. S.Y.A.C. and J.L. own shares of Certara Strategic Consulting, Princeton, NJ, USA.

AUTHOR CONTRIBUTIONS

J.L., N.P., S.Y.A.C., and J.R.S. wrote the manuscript. J.L., N.P., S.Y.A.C., and J.R.S. designed the research. J.L., N.P., S.Y.A.C., and J.R.S. performed the research. J.L., N.P., and J.R.S. analyzed the data.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.