Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › A Worked Joint Analysis Of Two Pharmacodynamic Endpoints

Joint Models of Multiple Endpoints · Tutorial 394 of 598

A Worked Joint Analysis Of Two Pharmacodynamic Endpoints

Work through a joint model for a continuous biomarker and a binary symptom response, including its likelihood, predictions, diagnostics, and interpretation.

Advanced 12 min read

What You'll Learn

  • Define endpoint-specific and joint targets for paired biomarker and symptom outcomes.
  • Specify a shared-random-effect model for continuous and binary endpoints.
  • Construct the joint likelihood and interpret its dependence parameter.
  • Calculate marginal symptom-response predictions from a probit model.
  • Compare observed patterns with model predictions and identify failure modes.
  • Translate endpoint-specific estimates and shared dependence into a cautious decision.

Why Analyze the Endpoints Together?

A treatment may change a laboratory biomarker and improve a symptom, but those outcomes need not move together in every participant. Separate analyses can estimate an average treatment effect on each endpoint. They do not, by themselves, describe how the paired outcomes co-vary or estimate whether the people with a larger biomarker response also tend to have symptom improvement. A joint analysis makes the paired structure explicit while preserving a distinct treatment effect for each endpoint.

Consider a hypothetical randomized study with 120 participants, 60 assigned to placebo and 60 to an active dose. At week 8, each participant has a continuous biomarker measurement and a binary symptom outcome: responder or non-responder. The analysis below focuses on this paired visit. It does not assume that a biomarker change is a substitute for symptom benefit; the two outcomes retain separate clinical interpretations.

Definition: An endpoint-specific estimand is a precisely defined target about one outcome, such as the treatment-versus-placebo difference in mean log biomarker change or in symptom-response probability. A joint estimand is a target involving the paired outcomes or their dependence, such as the probability of symptom response at a specified biomarker value or the association between the two outcomes. A joint model may estimate these targets together without collapsing them into a single endpoint.

Here define \(A_i=1\) for active treatment and \(A_i=0\) for placebo. Let \(B_i\) be participant \(i\)'s log week-8-to-baseline biomarker ratio; lower values indicate a greater biomarker reduction. Let \(S_i=1\) indicate symptom response at week 8 and \(S_i=0\) otherwise. Endpoint-specific targets include the treatment contrast in mean \(B_i\) and the treatment contrast in \(P(S_i=1)\). A dependence target describes whether unobserved participant-level variation associated with \(B_i\) is also related to the chance of symptom response.

A Joint Model for the Paired Outcomes

Use a shared random effect \(b_i\) to represent participant-level variation that may link the endpoints. The continuous biomarker model is normal, while the binary symptom model uses a probit latent-variable representation:

$$ \begin{aligned} B_i \mid b_i,A_i &\sim N(\beta_0+\beta_A A_i+b_i,\sigma^2),\\ S_i &= \mathbf{1}(L_i>0),\\ L_i \mid b_i,A_i &= \gamma_0+\gamma_A A_i+\kappa b_i+u_i,\\ b_i &\sim N(0,\tau^2),\qquad u_i\sim N(0,1). \end{aligned} $$

The indicator \(\mathbf{1}(L_i>0)\) equals 1 when its condition is true and 0 otherwise. The random effect \(b_i\) shifts the biomarker mean; \(\kappa\) determines how that same shift relates to the latent symptom-response propensity. Thus \(\kappa=0\) means the shared effect does not link the endpoints in this specification. The residuals for the biomarker and symptom latent variable are independent conditional on \(b_i\) and treatment. Participants are independent, and the model assumes the stated normal distributions, a common biomarker residual variance, and a correct probit link. In this example, outcomes are complete and treatment is randomized; missing outcomes or departures from randomization would require additional assumptions and handling.

Conditional on \(b_i\), the endpoint contributions multiply. The likelihood contribution for a participant with observed pair \((B_i,S_i)\) is therefore:

$$ L_i=\int \frac{1}{\sigma}\phi\!\left(\frac{B_i-\beta_0-\beta_A A_i-b}{\sigma}\right) \Phi(\gamma_0+\gamma_A A_i+\kappa b)^{S_i} \left[1-\Phi(\gamma_0+\gamma_A A_i+\kappa b)\right]^{1-S_i} \frac{1}{\tau}\phi\!\left(\frac{b}{\tau}\right)\,db. $$

Here \(\phi\) and \(\Phi\) are the standard normal density and distribution functions. The integral averages over the unobserved participant effect, rather than pretending it was observed. The joint likelihood is the product of \(L_i\) across participants. This is what makes the analysis joint: each participant's paired measurements enter one likelihood contribution through the same \(b_i\). Fitting two endpoint models separately would not estimate this shared dependence.

Worked Joint Analysis

Worked Example: Define the Targets and Summarize the Paired Data

Suppose the hypothetical study has the following week-8 summaries. The biomarker entries are group means of log ratios; symptom counts give responders out of 60.

GroupMean log biomarker ratioSymptom respondersObserved response proportion
Placebo-0.0427 of 6027 / 60 = 0.450
Active dose-0.2341 of 6041 / 60 = 0.683

The endpoint-specific biomarker contrast is \(-0.23-(-0.04)=-0.19\) log units. Exponentiating gives \(\exp(-0.19)=0.827\), so the active-to-placebo ratio of geometric biomarker ratios is about 0.827, or 17.3% lower under active treatment. The observed symptom-response risk difference is \(41/60-27/60=14/60=0.233\), a 23.3 percentage-point difference.

These are marginal group summaries; they do not show whether the participants with larger biomarker reductions are more likely to respond symptomatically. The joint target therefore includes the within-person association, represented here by \(\kappa\), in addition to the two treatment contrasts.

Worked Example: Interpret a Hypothetical Joint-Model Fit

Suppose maximum-likelihood fitting to all 120 paired records gives the following hypothetical estimates: \(\hat{\beta}_0=-0.04\), \(\hat{\beta}_A=-0.19\), \(\hat{\sigma}=0.30\), \(\hat{\tau}=0.25\), \(\hat{\gamma}_0=-0.10\), \(\hat{\gamma}_A=0.60\), and \(\hat{\kappa}=-1.20\). These are illustrative fitted results, not estimates reconstructed from the group summaries alone.

The estimate \(\hat{\beta}_A=-0.19\) is the treatment contrast in the model's mean log biomarker ratio. On the ratio scale, \(\exp(-0.19)=0.827\); this corresponds to an estimated 17.3% lower geometric ratio for active treatment than placebo. This is an endpoint-specific treatment effect, not evidence by itself of symptom benefit.

The negative \(\hat{\kappa}\) indicates that a participant-level shift toward a higher biomarker value is associated with a lower latent propensity for symptom response, conditional on treatment. Equivalently, the shared variation links lower biomarker values with greater symptom-response propensity. It is an association encoded by the model, not proof that changing the biomarker causes symptom improvement.

For an endpoint-specific coefficient check, the Wald statistic for the conditional symptom treatment coefficient is \(z=\hat{\gamma}_A/SE(\hat{\gamma}_A)\). If its hypothetical standard error is 0.28, then \(z=0.60/0.28=2.14\), with a two-sided normal-approximation \(p\)-value of about 0.032. This tests \(H_0:\gamma_A=0\) against a nonzero conditional probit treatment coefficient. It does not directly test every possible marginal symptom estimand.

From Conditional Effects to Marginal Predictions

A treatment decision often concerns the population response probability, not a coefficient conditional on an unobserved \(b_i\). In this probit model, the marginal probability can be derived directly. Since \(L_i=\gamma_0+\gamma_A A_i+\kappa b_i+u_i\), and independent normal variables satisfy \(b_i\sim N(0,\tau^2)\) and \(u_i\sim N(0,1)\), the random part \(\kappa b_i+u_i\) is normal with mean zero and variance \(\kappa^2\tau^2+1\). Therefore:

$$ P(S_i=1\mid A_i=a) =P(L_i>0\mid A_i=a) =\Phi\!\left( \frac{\gamma_0+\gamma_A a}{\sqrt{1+\kappa^2\tau^2}} \right). $$

This derivation matters because substituting \(b_i=0\) gives a conditional prediction for a participant at the average random effect, not the population-average response probability. Integrating over the random effect changes the scale of the probit argument.

Worked Example: Calculate Model-Based Symptom Predictions

Use the hypothetical estimates \(\hat{\gamma}_0=-0.10\), \(\hat{\gamma}_A=0.60\), \(\hat{\kappa}=-1.20\), and \(\hat{\tau}=0.25\). The marginal probit denominator is:

$$ \sqrt{1+\hat{\kappa}^{\,2}\hat{\tau}^{\,2}} =\sqrt{1+(1.20)^2(0.25)^2} =\sqrt{1.09}=1.044. $$

For placebo, the standardized argument is \(-0.10/1.044=-0.0958\), so the predicted response probability is \(\Phi(-0.0958)\approx0.462\). For active treatment, the argument is \((-0.10+0.60)/1.044=0.479\), giving \(\Phi(0.479)\approx0.684\). The model-based risk difference is \(0.684-0.462=0.222\), or 22.2 percentage points. These predictions are close to the observed proportions, 0.450 and 0.683, but the comparison is descriptive: fitted predictions and observations are not independent evidence.

Checking Fit and Choosing a Dependence Structure

Before using the estimates, check whether the model reproduces features of the paired data, not merely each endpoint's marginal distribution. Compare observed and predicted biomarker distributions and symptom proportions by treatment group. Also examine symptom-response rates across biomarker ranges within each group, and compare observed paired patterns with simulations from the fitted joint model. The earlier material on residual and predictive checks for joint models provides the general diagnostic framework; here the important question is whether the particular normal/probit shared-effect structure captures the observed pairing.

For example, if the model predicts that participants with lower biomarker values are more likely to respond, but observed response rates show no such pattern or show the reverse within treatment groups, the shared-effect structure may be inadequate. Other warning signs include systematically poor biomarker tail predictions, group-specific symptom calibration errors, estimates of \(\kappa\) that change sharply under reasonable model variations, or convergence and boundary problems. A good marginal fit alone does not establish a good joint fit.

A likelihood-ratio comparison can assess whether this model's shared link improves fit over the nested independence model \(\kappa=0\), provided regularity conditions hold and the models use the same marginal specification. Suppose the hypothetical fitted values are \(-2\ell=214.6\) for the shared model and \(-2\ell=220.8\) for the independence model. Then \(2(\ell_{\text{shared}}-\ell_{\text{independent}})=220.8-214.6=6.2\). With one added parameter, the reference \(\chi^2\) distribution has \(df=1\), giving \(p\approx0.0128\). This is evidence against independence within this model family; it does not establish that the chosen dependence structure is correct.

The shared-random-effect model is one practical choice, not a universal solution. A residual covariance model is more natural when both outcomes have continuous repeated measurements; a copula can represent other dependence patterns; and separate endpoint fits may be adequate when dependence is not a target and does not materially affect inference. These alternatives encode different assumptions. With one continuous and one binary endpoint, the shared-effect formulation has a clear interpretation, but may be restrictive if the association varies by treatment, is nonlinear, or has tail behavior not captured by normal latent variation.

Turning Results Into a Decision

Suppose a prespecified development rule requires at least a 10% reduction in the biomarker ratio and at least a 15 percentage-point increase in symptom response. The model's point estimates, a 17.3% lower ratio and a 22.2 percentage-point response increase, meet both numerical thresholds. The proper conclusion is that the fitted estimates meet these illustrative point-estimate criteria. A final decision should also account for uncertainty in the two treatment effects, the joint dependence estimate, safety evidence, and the clinical consequences of false-positive or false-negative decisions.

The shared effect informs whether the endpoints tend to move together across participants; it is not an extra treatment benefit to add to the two endpoint effects. Nor should the biomarker result be treated as proof of clinical response. If the biomarker improved but symptom response did not meet its criterion, the endpoint-specific findings would disagree and the decision rule should be applied as written. If both improve, their uncertainty and dependence matter when assessing a joint success rule.

Common pitfall: “The biomarker effect is significant, so the symptom endpoint must also benefit; the two endpoints are jointly modeled.” This is incorrect. Joint modeling represents dependence and can inform each estimate, but it does not make the endpoint-specific effects interchangeable. Report each target on its own scale, state the shared-association result separately, and apply the prespecified decision rule to the relevant estimates and uncertainty.

The main limitations are structural and practical. The normal random-effect distribution and probit link may not describe the data; a single shared effect may not capture all dependence; and estimates of association can be weakly identified when paired outcomes contain little information about their relationship. Missingness, intercurrent events, or deviations from randomized treatment assignment also affect the interpretation unless addressed in the estimand and analysis assumptions. Diagnostics and sensitivity analyses should therefore support, rather than replace, a clinically defined decision.

Check Your Understanding

Answer each question, then compare with the concise solutions.

  1. What does \(\beta_A\) represent in the biomarker model, and what is its unit?
  2. Why is \(\Phi(\gamma_0+\gamma_A a)\) not generally the marginal symptom-response probability?
  3. For \(\kappa=-1.20\), what direction of association does the shared effect imply between biomarker values and latent response propensity?
  4. If the observed symptom proportions are close to fitted proportions, is that alone enough to establish an adequate joint fit?

Solutions. (1) \(\beta_A\) is the active-versus-placebo difference in mean log biomarker ratio, in log-ratio units. (2) It conditions on \(b_i=0\); the population prediction integrates over \(b_i\), producing the denominator \(\sqrt{1+\kappa^2\tau^2}\). (3) Higher shared-effect shifts increase the biomarker but, because \(\kappa\) is negative, decrease latent response propensity; lower biomarker values are associated with greater response propensity. (4) No. Marginal agreement does not show that the model reproduces paired dependence, so joint predictions and within-group paired patterns must also be checked.

Clinical Biostats