Causal Inference for First Non‐Fatal Events With the Competing Risk of Death: A Principal Stratification Approach

TL;DR

This paper introduces a new statistical method to estimate the direct effect of treatments on nonfatal events in clinical trials, accounting for the competing risk of death.

Contribution

The novel contribution is the development of a principal stratification framework to estimate direct treatment effects on nonfatal events in the presence of death.

Findings

The proposed model estimates the principal stratum hazard ratio, which reflects the direct treatment effect on nonfatal events.

Simulation studies confirm the reliability of the new estimators.

The method is illustrated using a heart failure clinical trial.

Abstract

In clinical trials involving both mortality and morbidity, an active treatment can influence the observed risk of the first nonfatal event either directly, through its effect on the underlying nonfatal event process, or indirectly, through its effect on the death process, or both. Discerning the direct effect of treatment on the underlying first nonfatal event process holds clinical interest. However, with the competing risk of death, the Cox proportional hazards model that treats death as non‐informative censoring and evaluates treatment effects on time to the first nonfatal event provides an estimate of the cause‐specific hazard ratio, which may not correspond to the direct effect. To obtain the direct effect on the underlying first nonfatal event process, within the principal stratification framework, we define the principal stratum hazard and introduce the proportional principal stratum hazards model. This model estimates the principal stratum hazard ratio, which reflects the direct effect on the underlying first nonfatal event process in the presence of death and simplifies to the hazard ratio in the absence of death. The principal stratum membership is identified probabilistically using the shared frailty model, which assumes independence between the first nonfatal event process and the potential death processes, conditional on per‐subject random frailty. Simulation studies are conducted to verify the reliability of our estimators. We illustrate the method using the Carvedilol Prospective Randomized Cumulative Survival trial, which involves heart‐failure events.

Figures

Tables

TABLE 1. Comparison of estimation results between the PS and the cause‐specific approach.

	Summary statistics												Estimation of β
	Placebo						Active treatment						Hypothetical		Cause‐specific		Principal stratum
γ	Dead	Censored	LOF	D‾	D˜	Event	Dead	Censored	LOF	D‾	D˜	Event	Bias (SE)	HR	Bias (SE)	HR	γ˜	Bias (SE)	HR
ν0=0.25
0.5	29%	67%	5%	1.6	2.0	60%	25%	70%	5%	1.7	2.0	39%	0.002 (0.005)	0.5	0.054 (0.005)	0.53	0.5	−0.003 (0.006)	0.50
																	2.0	0.039 (0.006)	0.52
																	5.0	0.046 (0.005)	0.52
2.0	35%	60%	5%	1.6	2.0	80%	30%	65%	5%	1.6	2.0	58%	−0.003 (0.004)	0.5	0.008 (0.005)	0.50	0.5	−0.044 (0.005)	0.48
																	2.0	−0.010 (0.005)	0.50
																	5.0	0.001 (0.005)	0.50
5.0	37%	58%	5%	1.5	2.0	84%	31%	64%	5%	1.6	2.0	67%	−0.001 (0.004)	0.5	0 (0.004)	0.50	0.5	−0.045 (0.004)	0.48
																	2.0	−0.018 (0.004)	0.49
																	5.0	−0.007 (0.004)	0.50
ν0=0.4
0.5	37%	58%	4%	1.5	2.0	57%	25%	70%	5%	1.7	2.0	40%	0.002 (0.005)	0.5	0.107 (0.005)	0.56	0.5	0.002 (0.006)	0.50
																	2.0	0.071 (0.006)	0.54
																	5.0	0.088 (0.005)	0.55
2.0	48%	48%	4%	1.4	1.8	75%	30%	65%	5%	1.6	2.0	59%	−0.003 (0.004)	0.5	0.034 (0.005)	0.52	0.5	−0.067 (0.005)	0.47
																	2.0	−0.008 (0.005)	0.50
																	5.0	0.015 (0.005)	0.51
5.0	51%	45%	4%	1.4	1.7	80%	31%	64%	5%	1.6	2.0	67%	−0.001 (0.004)	0.5	0.012 (0.004)	0.51	0.5	−0.078 (0.005)	0.46
																	2.0	−0.030 (0.004)	0.49
																	5.0	−0.006 (0.004)	0.50

TABLE 2. Sensitivity analysis comparing the principal stratum and cause‐specific approaches.

	Summary statistics												Estimation of β
	Placebo						Active treatment						Hypothetical		Cause‐specific		Principal stratum
γ	Dead	Censored	LOF	D‾	D˜	Event	Dead	Censored	LOF	D‾	D˜	Event	Bias (SE)	HR	Bias (SE)	HR	γ˜	Bias (SE)	HR
ν0=0.25
0.5	30%	65%	5%	1.6	2.0	71%	26%	69%	5%	1.7	2.0	43%	−0.887 (0.005)	0.41	0.078 (0.005)	0.45	0.5	0.016 (0.006)	0.42
																	2.0	0.070 (0.005)	0.44
																	5.0	0.081 (0.005)	0.45
2.0	35%	60%	5%	1.6	2.0	82%	30%	65%	5%	1.6	2.0	60%	−0.748 (0.004)	0.47	0.020 (0.005)	0.48	0.5	−0.032 (0.005)	0.46
																	2.0	0.004 (0.005)	0.48
																	5.0	0.017 (0.005)	0.48
5.0	37%	58%	5%	1.5	2.0	85%	31%	64%	5%	1.6	2.0	67%	−0.713 (0.004)	0.49	0.009 (0.004)	0.49	0.5	−0.036 (0.005)	0.47
																	2.0	−0.008 (0.004)	0.49
																	5.0	0.004 (0.004)	0.49
ν0=0.4
0.5	40%	56%	4%	1.5	2.0	68%	26%	69%	5%	1.7	2.0	44%	−0.887 (0.005)	0.41	0.124 (0.005)	0.47	0.5	0.009 (0.006)	0.42
																	2.0	0.095 (0.005)	0.45
																	5.0	0.116 (0.005)	0.46
2.0	48%	48%	4%	1.4	1.8	77%	30%	65%	5%	1.6	2.0	60%	−0.748 (0.004)	0.47	0.040 (0.005)	0.49	0.5	−0.061 (0.005)	0.45
																	2.0	0.002 (0.005)	0.47
																	5.0	0.026 (0.005)	0.49
5.0	51%	45%	4%	1.4	1.7	80%	31%	64%	5%	1.6	2.0	68%	−0.713 (0.004)	0.49	0.019 (0.004)	0.50	0.5	−0.071 (0.005)	0.46
																	2.0	−0.021 (0.005)	0.48
																	5.0	0.003 (0.004)	0.49

TABLE 3. Estimation of treatment effects for the first nonfatal event from the principal stratum (PS) and the cause‐specific (CS) approach.

Approach	γ˜	HR	95% CI	Proportionality
PS	0.25	0.781	(0.665, 0.919)	p=0.64
	0.5	0.791	(0.676, 0.927)	p=0.61
	1	0.804	(0.692, 0.934)	p=0.68
	2	0.813	(0.702, 0.941)	p=0.76
	5	0.818	(0.708, 0.948)	p=0.81
	10	0.820	(0.711, 0.949)	p=0.84
CS	∞	0.821	(0.716, 0.943)	—

Full text

Introduction

1

In clinical trials involving both mortality and morbidity, the treatment effect on the observed risk of the first nonfatal event is commonly evaluated using the cause‐specific (CS) hazard ratio. This ratio is estimated using the Cox proportional hazards model that treats death as non‐informative censoring. Since the death process and the first nonfatal event process within the same subject are correlated—and because death prevents the occurrence of a nonfatal event—the active treatment may influence the CS hazard ratio directly by affecting the underlying nonfatal event process, indirectly through the death process, or both. This correlation complicates and introduces uncertainty in the interpretation of the CS hazard ratio. Therefore, except in cases where it is confidently known that the active treatment does not have differential effects on mortality, the “direct effect” on the underlying first nonfatal event process that is not mediated by competing risks of death—a modification of the underlying mechanism that produces the first nonfatal event—is desired [1]. A formal mathematical definition of direct effect will be provided in the following section.

Frangakis and Rubin [2] proposed the principal stratification framework to account for post‐randomization outcomes such as death. We focus on the principal stratum (PS) estimand for addressing the competing risk of death in a randomized clinical trial. Principal stratification is a partition of the population into subpopulations (strata) based on joint values of the competing outcomes under all treatment conditions. In our context, the PS estimand targets subjects who would survive to a specific time point t, irrespective of their treatment assignment, within the intended study population. This subpopulation constitutes the PS of interest, termed “Always Survivors” defined at t, and is denoted as 𝒜(t). The treatment effect within 𝒜(t) is often referred to as the survivor average causal effect (SACE) [3].

Typically, 𝒜(t) is defined using a limited number of time points, and the SACE is presented across 𝒜(t) defined at varying t [4, 5]. This mirrors milestone survival analysis, which compares survival probabilities at prespecified time points rather than considering the entire curve [6]. While milestone analysis is often advocated in immunotherapy trials, where non‐proportionality is common and significantly reduces the power of the Cox model, simulation studies suggest the Cox model generally has higher power unless there is a significant treatment effect delay [7]. Moreover, presenting snapshot effects at specific time points does not provide decision‐makers with insights into the sensitivity of conclusions to the potentially arbitrary selection of t. In trials where violation of the proportionality assumption is not obvious, we expect the Cox model, capturing treatment effects over time, to have higher power and be more informative than snapshot effects. To capture treatment effects across 𝒜(t) defined at different t, we extend the hazard function and define the PS hazard as the instantaneous probability of experiencing the nonfatal event given that the subject has not experienced the nonfatal event and death and would survive to t under the counterfactual arm. Additionally, we propose the proportional principal stratum hazards (PPSH) model, which maintains the structure of the Cox model but replaces the hazard function with the PS hazard function. The PPSH model assumes and estimates a constant hazard ratio shared across 𝒜(t) defined at different t. If the PPSH model is correctly specified, the PS hazard ratio reflects the direct effect on the underlying first nonfatal event process, conditional on being in 𝒜(t). We refer to this as the “conditional direct effect,” as will be discussed in more detail later.

Unfortunately, we cannot directly observe 𝒜(t) due to its counterfactual nature: if a subject survived to time t under one arm, we cannot determine if they would have survived to time t under the other arm. However, correctly estimating the probability of each subject belonging to 𝒜(t) (referred to as PS probability) enables us to estimate the PS hazard ratio. Estimating PS probability requires further assumptions. One commonly used assumption is the monotonicity assumption [8], which, in our context, implies that for every subject, their potential time to death in the active treatment arm is no shorter than in the placebo arm. While the monotonicity assumption may be scientifically plausible in certain situations, it is important to note that an active treatment delaying mortality for the population does not necessarily mean it delays mortality for every subject. It remains uncertain whether the monotonicity assumption is plausible in the presence of death. Moreover, the monotonicity assumption alone is not sufficient to estimate PS probabilities. Additional assumptions are typically required [9, 10, 11].

The challenge in estimating PS probabilities is that the potential death time from the counterfactual arm and the first nonfatal event time cannot be jointly observed, preventing direct identification of the joint survival function. To address this issue, the shared frailty model offers a potential solution. In this model, each subject is assigned a random effect, referred to as frailty. The potential death time from the counterfactual arm and the first nonfatal event time are assumed to be independent, conditional on the frailty. Similar assumptions have been adopted in previous studies to factorize the joint density of potential outcomes [4, 5].

The paper is structured as follows: Section 2 introduces the PS estimand within our proposed causal estimand framework. Section 3 discusses the PPSH model and the identification of PS membership. Section 4 presents a simulation study to verify our proposed method. Section 5 demonstrates the application of our method in a real cardiovascular trial. Finally, we conclude with a discussion in Section 6.

Estimand

2

We focus on clinical trials where patients are randomly assigned to an active treatment or a placebo. In this paper, with no loss of generality and to avoid extra notation, we will subsequently assume that we are already within cells defined by a specific value of covariates.

Assume there are n subjects (i=1,…,n) in the study. Each subject i is assigned to either the active treatment arm (Zi=1) or the placebo arm (Zi=0). Following causal inference conventions, we impose the Stable Unit Treatment Value Assumption (SUTVA) [12]. SUTVA posits that there is only one version of the active treatment; otherwise, we would need to define potential outcomes corresponding to each version. Additionally, SUTVA assumes no interference between subjects, implying that each potential outcome for the ith subject is determined solely by the treatment received by that subject and does not vary with treatments assigned to others. Let Yiz and Tiz be independent and identically distributed (IID) realizations of Yz and Tz, representing the potential time to death and potential time to the first nonfatal event under treatment z, respectively. Since death precludes the occurrence of a nonfatal event, Tiz is subject to truncation of Yiz [13]. Let Yi and Ti be the IID realizations of Y and T, representing the observed outcomes.

Assumption 1Under SUTVA, the observable variables Yi and Ti satisfy Yi=Yiz and Ti=Tiz if Zi=z.

Causal Estimand

2.1

Hernán and Robins [14] define a causal effect as a contrast of any function of the distributions of counterfactual outcomes under different actions or treatment values. This definition can be expressed mathematically as follows. A parameter β is a causal estimand if there exists a monotone function h(β) and functional G(F) such that:

$$[eqn]$$

where F1 and F0 represent the cumulative distribution functions (CDFs) of counterfactual outcomes of interest, possibly as functions of some baseline covariates. This definition requires that the counterfactual outcome of interests must be totally ordered and defined for all subjects. Different choices of G(·) result in different causal estimands.

Let Fz(t)=P(Tz≤t) represent the CDF of the potential first nonfatal event time in the absence of death. Fz(t) can be functions of baseline covariates. According to the definition (1), G{F1(t)}−G{F0(t)} is a causal estimand for treatment effects on the first nonfatal event in the absence of death and for direct effects on the underlying first nonfatal event process in the presence of death.

In the absence of death, the hazard ratio in the Cox model that treats the first nonfatal event T as the outcome is a causal estimand in a randomized trial, provided that the Cox model is correctly specified. Due to randomization, the pair (T1,T0) is independent of the treatment assignment Z, then in the Cox model, we have Fz(t)=F(t|Z)=1−exp{−exp(βZ)Λ0(t)} where Λ0(t) is the baseline cumulative hazard function. Let G{F(t|Z)}=log[−log{1−F(t|Z)}]=logΛ0(t)+βZ. Then G{F1(t)}−G{F0(t)}=G{F(t|Z=1)}−G{F(t|Z=0)}=β. This may not hold in an observational study, where a patient's prognosis may influence the chosen treatment option. If Fz(t) depends on baseline covariates and proportional hazards hold, then β represents the conditional hazard ratio, which also has a causal interpretation.

Cause‐Specific Hazard Ratio

2.2

In the presence of death, Tiz is subject to truncation of Yiz. A subject is considered at risk for the first nonfatal event at time t if they have not experienced the nonfatal event (T≥t) and have not died (Y>t). Accordingly, the hazard can be expressed as limΔt→0P(t≤T<t+Δt|T≥t,Y>t)/Δt, termed as the CS hazard. The CS hazard ratio can then be represented as:

$$[eqn]$$

In a randomized trial, this ratio can be estimated by:

$$[eqn]$$

Let Fz∗(t)=P(Tz≤t|Yz>t) represent the CDF of the potential first nonfatal event time in the presence of death. By defining G{Fz∗(t)}=log[−log{1−Fz∗(t)}], the contrast G{F1∗(t)}−G{F0∗(t)} corresponds to the CS hazard ratio in the Cox model. Due to correlation between Yz and Tz, Fz∗(t)≠Fz(t). As a result, with the same transformation G(·), G{F1∗(t)}−G{F0∗(t)}≠G{F1(t)}−G{F0(t)}. Thus, the CS hazard ratio in the presence of death does not correspond to direct effects on the underlying first nonfatal event process, unless the population is homogeneous, implying Yz and Tz are independent and Fz∗(t)=Fz(t).

If the active treatment exhibits no differential effects on survival (i.e., Y1=Y0), then

$$[eqn]$$

Since the potential death times (Y1,Y0) are unaffected by the actual treatment assignment, conditioning on Y1>t and Y0>t is akin to conditioning on a baseline covariate. Thus, the CS hazard ratio reflects the direct effect conditional on being in 𝒜(t) (the conditional direct effect). However, if the active treatment has differential effects on survival, this hazard ratio does not correspond to the conditional direct effect and lacks a causal interpretation. Nonetheless, this observation motivates us to define the PS hazard ratio, as outlined below.

Principal Stratum Hazard Ratio

2.3

Under SUTVA, at each time t, subjects can be classified into four strata according to a pair of counterfactual death times (Yi1,Yi0):

• Always Survivors: {i|Yi1>t,Yi0>t}, subjects who would survive to t regardless of treatment assignment.
• Never Survivors: {i|Yi1≤t,Yi0≤t}, subjects who would not survive to t regardless of treatment assignment.
• Active Survivors: {i|Yi1>t,Yi0≤t}, subjects who would survive to t under the active treatment arm, but would not survive to t under the placebo arm.
• Placebo Survivors: {i|Yi1≤t,Yi0>t}, subjects who would survive to t under the placebo arm, but would not survive to t under the active treatment arm.

Among these four strata at t, our primary focus is on the “Always Survivors” stratum 𝒜(t). Because subjects in 𝒜(t) can survive to t regardless of the treatment, there is no competing risk of death. We define the PS hazard within 𝒜(t), denoted as λP(t), as the instantaneous probability of experiencing the nonfatal event, given that the subject has not experienced either the nonfatal event or death, and would survive to time t under the counterfactual arm. Mathematically, this is represented as:

$$[eqn]$$

Note that if there is no mortality in the trial, i.e., P(Y1>t)=P(Y0>t)=1, then λP(t) becomes the well‐known hazard function limΔt→0P(t≤T<t+Δt|T≥t)/Δt.

The PS hazard ratio can be expressed as

$$[eqn]$$

In the absence of death, when P(Y1>t)=P(Y0>t)=1, the PS hazard ratio simplifies to the hazard ratio. In the presence of death, (Y1>t,Y0>t) defines the hazard ratio within 𝒜(t), where there is no competing risk of death. The PS membership at t is determined based on counterfactual outcomes Y1 and Y0, which are not affected by treatment assignment and therefore can be regarded as a baseline covariate. Thus, the PS hazard ratio in 𝒜(t) at time t, where there is no competing risk of death, can be considered as the conditional hazard ratio with PS membership at time t as the baseline covariate, and therefore, reflects the conditional direct effect on the underlying first nonfatal event process.

The PS hazard ratio shares a similar idea with the concept of the causal hazard ratio introduced by Martinussen et al. [15], which is defined within the subpopulation of those who will potentially survive until a specific time, regardless of the treatment received. While Martinussen et al. [15] focused on the treatment effect on mortality, our focus is on the direct effect on the underlying first nonfatal event process. Both definitions are related to the SACE [3]. Further discussion on SACE can be found in Section 6.

Following the structure of the Cox model, the PPSH model specifies that

$$[eqn]$$

where λ0P(t) is the (marginal) baseline PS hazard function, and

$$[eqn]$$

which provides an estimate of the PS hazard ratio in a randomized trial. If the PPSH model is correctly specified, the PS hazard ratio, which reflects the conditional direct effect on the underlying first nonfatal event process, has a causal interpretation.

Estimation of the Principal Stratum Hazard Ratio

3

Principal Stratum Membership

3.1

The actual identification of 𝒜(t) is generally not possible. However, with assumptions and modeling, we can estimate the probability of each subject i being in 𝒜(t) (PS probability). Here, we introduce one potential approach to estimate this probability.

Let Ci be the IID realizations of the censoring time C. We assume censoring time is independent of potential outcomes:

Assumption 2 (T1,T0,Y1,Y0)╨C

Let 𝒮 represent the set of the first nonfatal event times from all subjects. Suppose there are m event times within the set 𝒮. We denote these unique event times as 0<t1<t2<⋯<tm. At each tj∈𝒮, we define ℛj as the at‐risk set at tj, which comprises subjects who remain in the study (Yi>tj,Ci>tj) and have not yet experienced a nonfatal event by time tj (Ti≥tj). Each subject i belonging to ℛj, whose last known follow‐up time is denoted by Di=min(Yi,Ci) with realization di, where di>tj, falls into one of the following two cases:

1. Ti=tj: The first nonfatal event occurred at tj.
2. Ti>tj: The first nonfatal event has not yet occurred by tj.

Subject i belonging to ℛj is associated with a PS probability of being in 𝒜(tj), denoted by pij. Suppose, without loss of generality, that a subject i is assigned to arm z.

If the subject i falls into case 1, the PS probability at tj is:

$$[eqn]$$

As previously mentioned, the joint distribution of Y1−z and Tz cannot be directly identified because these potential outcomes are correlated but never jointly observed. The correlation arises because both Y1−z and Tz are outcomes from the same subjects. The correlation can be modeled by a subject‐level random effect. Specifically, we assume the correlation between Y1−z and Tz is captured by the variability of this random effect (e.g., its variance), allowing us to model them as conditionally independent given the random effect. In time‐to‐event analysis, such random effects are known as frailties, and the resulting models are referred to as shared frailty models [16]. To factorize the joint distribution appearing in the numerator of Equation (3), we adopt the following conditional independence assumption:

Assumption 3Conditional on the frailty θ, the potential death time from the counterfactual arm is independent of the potential first nonfatal event time: Y1−z╨Tz|θ, for z=0,1.

While this conditional independence assumption differs slightly from those in Lyu et al. [4] and Comment et al. [5], it shares a similar motivation. In their frameworks, the assumption is that the potential outcomes under different treatment arms are independent conditional on θ, i.e., (Y1−z,T1−z)╨(Yz,Tz)|θ, which appears to be a somewhat stronger assumption than ours.

Using the imposed conditional independence assumption, we can factorize Equation (3), as follows:

$$[eqn]$$

where θ is the frailty and p(θ) is the density function of θ.

Let ηYz(t) and ηTz(t) denote the conditional cumulative hazard of Yz and Tz given θ=1, then

$$[eqn]$$

For an explicit expression for the PS probability, with practical considerations, we assume θ∼Γ(γ,γ) where the mean is 1, and the variance is 1/γ. Then:

$$[eqn]$$

Similarly, if the subject i falls into case 2, the PS probability at tj is:

$$[eqn]$$

Assuming a distribution for θ may allow us to derive conditional cumulative hazards ηYz(t) and ηTz(t) from marginal functions P(Yz>t) and P(Tz>t|Yz>t). As above, assuming θ∼Γ(γ,γ) and let SYz(t)=P(Yz>t), we have:

$$[eqn]$$

Then, the conditional cumulative hazard of Yz given θ = 1 is

$$[eqn]$$

The conditional cumulative hazard of Tz given θ = 1 is derived similarly, but we must account for the competing risk of death. Let STz(t|Yz>t)=P(Tz>t|Yz>t), we have

$$[eqn]$$

If we assume P(Tz>t|θ,Yz>t)=P(Tz>t|θ), then the expression simplifies to

$$[eqn]$$

where p(θ|Yz>t) is available in closed form only if θ follows a gamma distribution [16].

Using this result, the conditional cumulative hazard of Tz given θ = 1 is

$$[eqn]$$

Therefore, obtaining a closed‐form expression for ηTz(t) requires the assumption P(Tz>t|θ,Yz>t)=P(Tz>t|θ), which implies that, conditional on θ, the distribution of Yz is independent of Tz. This leads to the following assumption:

Assumption 4Conditional on the frailty θ, the potential death time is independent of the potential first nonfatal event time: Yz╨Tz|θ, for z=0,1.

This assumption has also been used in Lyu et al. [4] and Comment et al. [5] to factorize their joint likelihood functions.

Assumptions 3 and 4 together imply that the potential first nonfatal event time is conditionally independent of the potential death times under both treatment arms, given that the frailty appropriately accounts for the correlation. Notably, we do not impose any conditional independence assumption between the two potential death times, Yz and Y1−z, across arms.

The marginal functions SYz(t) and STz(t|Yz>t) are not directly observable but, in a randomized trial, they can be estimated by SY(t|Z=z) and ST(t|Y>t,Z=z), respectively. This is because randomization guarantees independence between treatment assignment Z and potential outcomes. Therefore, we have:

$$[eqn]$$

Similarly, we have ST(t|Y>t,Z=z)=STz(t|Yz>t).

To estimate SY(t|Z=z), in practice, we fit a Cox model where we view death Y as the outcome and the treatment assignment Z as the covariate. ST(t|Y>t,Z=z) can be estimated non‐parametrically as ∑Zi=zI(Ti>t)/∑Zi=zI(Yi>t), where I(·) is the indicator function.

We assume θ∼Γ(γ,γ) and use θ to induce correlation between Tz and Y1−z. However, Tz and Y1−z are never observed jointly, so we need to prespecify some values of γ when calculating the PS probabilities. The estimation of PS hazard ratio, which will be addressed in Section 3.2, needs estimated PS probabilities as input. Therefore, we can expect the estimated PS hazard ratio to vary with specified γ. While one might consider estimating γ—since θ also induces correlation between the observed nonfatal event and death times—the goal of the PPSH model is not to produce a single point estimate, but rather to complement the CS hazard ratio with a range of plausible PS hazard ratios based on prespecified values of γ. Moreover, estimating γ faces numerous technical issues, as detailed in Appendix A.

In summary, we present one potential approach to estimate PS probabilities under the conditional independence assumption and the independent censoring assumption in a randomized trial. The estimation consists of three steps:

1. Fit the Cox model for death time Y to obtain the estimated ŜY(t|Z=z). Estimate ST(t|Y>t,Z=z) non‐parametrically as ŜT(t|Y>t,Z=z)=∑Zi=zI(Ti>t)/∑Zi=zI(Yi>t). ŜY(t|Z=z) and ŜT(t|Y>t,Z=z) serve as the estimates of the marginal functions SYz(t) and STz(t|Yz>t), respectively.
2. Estimate the conditional cumulative hazard functions ηYz(t) and ηTz(t) from the estimated marginal functions ŜYz(t) and ŜTz(t|Yz>t) using Equations (4) and (5).
3. Collect the first nonfatal event times from all subjects into set 𝒮. At each tj∈𝒮, we calculate the PS probability for each subject in ℛj={i|Yi>tj,Ci>tj,Ti≥tj}, using the proposed formula with η^Yz(t) and η^Tz(t) from step 2 and the prespecified γ. The formula to use depends on which case the subject belongs to.

Proportional Principal Stratum Hazards Model

3.2

The likelihood function of the PPSH model resembles that of the Cox model but incorporates certain modifications. Following the notation in Section 3.1, suppose there are m unique event times 0<t1<t2<⋯<tm within the set 𝒮. We first consider the scenario where there are no tied events. That is, only one subject had a nonfatal event at tj, and we index this subject by (j). Additionally, we assume that this subject belongs to the treatment arm z. Let 𝒱={tj∈S:Y(j)1−z>tj}. In other words, for each tj, where j=1,…,m, if the subject (j) who experienced a nonfatal event at tj belongs to 𝒜(tj), tj is a member of 𝒱; otherwise, it is not. Since the PS hazard is defined within 𝒜(tj), the partial likelihood (PL) is constructed within 𝒱. Suppose the value of pij is known, we have:

$$[eqn]$$

The second equation holds due to the structure of the PPSH model (2), where the baseline λ0P(t) is dispensed. Then the log PL is

$$[eqn]$$

In the case of tied events (more than one subject had a nonfatal event at tj), Breslow suggests using the same expression for pl(β) [17].

Let

$$[eqn]$$

The PL score function is

$$[eqn]$$

where I{Y(j)1−z>tj} is the indicator that the subject (j) belongs to 𝒜(tj).

We can estimate β by solving U(β)=0. However, 𝒱 is non‐identifiable because we cannot observe I{Y(j)1−z>tj}. To address this non‐identifiability issue, we introduce the modified score function:

$$[eqn]$$

We can demonstrate that, for j:tj∈𝒮,

$$[eqn]$$

where the second equation holds because of the conditional independence assumption.

Therefore, instead of solving U(β)=0, we can estimate β by solving U∗(β)=0. The standard Newton‐Raphson algorithm can be employed to find a solution to U∗(β)=0. Specifically, starting with an initial guess β^(0), the algorithm iteratively computes

$$[eqn]$$

until convergence, where ℒ(β) is the negative second derivative of log PL:

$$[eqn]$$

Given that the objective function exhibits strict concavity, the algorithm is likely to maintain numerical stability and converge quickly toward β^. Convergence problems are very rare using the default initial value of β^(0)=0 [18].

In summary, the estimation of the PS hazard ratio consists of two stages. In stage 1, at each nonfatal event time, we estimate PS probabilities for each at‐risk subject. In stage 2, a PPSH model, essentially a weighted Cox model with PS probabilities as weights, is fitted. Since PS probabilities are estimated rather than observed, we replace them with the estimates from stage 1.

To address the variability introduced by estimating PS probabilities, we construct bootstrap confidence intervals (CI) using the following steps:

1. Resample subjects with replacement to generate a random bootstrap sample of the same size as the original dataset. Resampling with replacement will include duplicate records from the same subject, and we treat each sampled record as coming from a distinct subject for bootstrap purposes. Accordingly, we assign a unique new subject ID to each record in the bootstrap sample, regardless of whether it originates from the same subject in the original data.
2. Estimate the PS hazard ratio from the bootstrap sample. This includes estimating PS probabilities and fitting the PPSH model using those estimates.
3. Repeat steps 1 and 2 for a total of B bootstrap replicates. Let {β^(1),β^(2),…,β^(B)} denote the bootstrap estimates of the PS hazard ratio. Denote β(α/2)∗ and β(1−α/2)∗ as the 100(α/2)% and 100(1−α/2)% of the bootstrap estimates. The resulting 100(1−α)% bootstrap CI is (β(α/2)∗,β(1−α/2)∗), using the percentile method [19].

Empirical evidence supporting the validity of this bootstrap procedure is provided in Appendix B.

The estimated PS probabilities in Section 3.1 are conditional on specified γ. As a result, the estimated PS hazard ratio is also conditional on γ. With different specified γ, we can expect different estimates of β. Sensitivity analysis is recommended to provide researchers with a range of causal estimates β^ under different values of γ.

A key assumption of the PPSH model is proportionality. A formal proportionality test is presented in Appendix C.

Simulations

4

Settings

4.1

Assuming that θi follows a gamma distribution with a mean of one and a variance of 1/γ, we model the potential death time Yiz, given θi, for simulation purposes, as an exponential distribution with a mean of 1/θiνz. This yields a survival function denoted as S(t;θi)=exp(−θiνzt). The rates for the placebo and active treatment arms are denoted as ν0 and ν1, respectively. Additionally, we consider the independent loss to follow‐up C, which follows an exponential distribution with a mean of 1/νc. The maximum follow‐up time for patients is denoted as τ. The last known follow‐up time is then calculated as D=min(Y,C,τ).

The model (2) focuses on the marginal PS hazard function λP(t|Z). By defining ηT0(t)=ϕt, and

$$[eqn]$$

where νY=ν0+ν1, we achieve λP(t|Z=1)λP(t|Z=0)=r, representing a marginal PPSH ratio of r for the first nonfatal event. For detailed derivation of ηT1(t), please refer to the Appendix D.

To generate nonfatal event time Tiz from ηTz(t), we solve θηTz(t)+logU=0, where U∼Unif(0,1). This approach works because −logU∼Exp(1), and the term 1−exp{−θηTz(t)} represents the CDF of an Exp(1) distribution. If the generated time Ti>Di, there is no observed nonfatal event during the follow‐up for subject i.

This simulation algorithm introduces correlations between (Yz,Y1−z,Tz) through θ. The three processes are then correlated but are mutually independent conditional on θ.

We also generate a separate hypothetical dataset where there is no mortality by setting ν0=ν1=0, following the same simulation algorithm as described earlier, with the only exception that each subject's last follow‐up time is now calculated as D=min(C,τ). As discussed, when there is no mortality, the PS hazard function simplifies to the well‐known hazard function. Therefore, for this hypothetical dataset, we could expect the Cox model to provide an estimated hazard ratio of r.

Results

4.2

In Table 1, we present a comparison of estimation results between our proposed PS approach and the CS approach, with the estimation results from the hypothetical dataset as a reference. Estimation on the hypothetical dataset is conducted by treating the nonfatal event as the outcome and applying the Cox model. To obtain the CS estimate, we also use the Cox model and treat the nonfatal event as the outcome, but relabel death as censoring.

We use common parameter values for this comparison: ν1=0.2, νc=0.03, τ=2, ϕ=2, r=0.5. The sample size is 300 with equal allocation. We explore different values of ν0 and γ. With ν1=0.2, ν0=0.4 indicates remarkable effects on mortality (hazard ratio of 0.5 conditional on θ), while ν0=0.25 represents a more realistic yet impressive effect (hazard ratio of 0.8 conditional on θ). γ=5 denotes a highly homogeneous population, while γ=0.5 signifies high heterogeneity. For each parameter configuration, we conduct 1000 simulations and provide summary statistics for the datasets used to estimate β.

In Table 1, when increasing ν0 from 0.25 to 0.4—thereby raising the death rate in the placebo arm—we observe a smaller proportion of subjects experiencing a nonfatal event in that arm, for the same value of γ and with r=0.5 held constant. This occurs because the nonfatal event process is subject to the competing risk of death. As ν0 increases, the time to death in general decreases, leading to more subjects dying before they have a chance to experience a nonfatal event. This demonstrates that even if the treatment effect on the underlying nonfatal event process remains the same, altering the treatment effect on the death process can change the observed risk ratio for the first nonfatal event.

The hypothetical dataset is generated separately using the same parameters except that ν0=ν1=0. We define the Monte Carlo bias as the Monte Carlo sample mean of estimates minus the true theoretical value log(0.5). The HR refers to the estimated hazard ratio, computed as the exponential of the Monte Carlo sample mean of estimates. As expected, the Monte Carlo bias from the hypothetical dataset is minimal, and the HR is close to its nominal value.

In the ideal scenario where the true value of γ is known, the Monte Carlo bias in our proposed PS approach is also minimal across all parameter configurations, likely due to random variability. In reality, where the true value of γ is unknown, misspecifying γ˜ in the PS estimation leads to bias, and the bias increases as γ˜ deviates further from the true value.

The bias increases for both PS and CS methods with greater treatment effects on mortality. The relative magnitude of bias between the two approaches is similar across different ν0 but varies with different γ. When the population is highly heterogeneous (γ=0.5), bias in CS estimates is pronounced and greater than that from the PS estimates, even when specifying γ˜=5, far from the true value of 0.5. However, bias in CS estimates is less pronounced as γ increases. In a highly homogeneous population (γ=5), bias in CS estimates becomes minimal and comparable to that from the PS approach.

This phenomenon is understandable considering the origin of bias in the CS approach. Each subject in our simulations has a frailty, with those having higher frailty experiencing mortality and the nonfatal event earlier. The CS approach compares subjects conditional on survival, but since treatment prolongs survival, those in the treatment arm are expected to have higher average frailty and thus a greater risk of the nonfatal event. Consequently, CS estimates tend to be greater than the true value. As treatment effects on mortality become more substantial, the disparity in frailty between the two arms after conditioning on survival becomes more pronounced, leading to increased bias. However, in highly homogeneous populations, the frailty of each subject is less distinguishable, resulting in comparable frailty between the two arms even after conditioning on survival, and thus, comparable risk of the nonfatal event. Therefore, CS estimates are closer to the true value in such scenarios.

In Table 1, with a heterogeneous population (γ=0.5), the PS hazard ratio approaches the CS hazard ratio as γ˜ increases. Setting a large value for γ˜ in the PS approach assumes a highly homogeneous population with similar frailty. In a homogeneous population, even after conditioning on the death time, survivors remain comparable in frailty between both arms. Therefore, for this assumed highly homogeneous population, the PS approach would yield an estimated HR similar to that of the CS approach.

Sensitivity Analysis

4.3

When calculating the PS probabilities, we assume the frailty follows a gamma distribution. This choice is common for several reasons. First, the conditional distribution of frailty, given survival, is only available in closed form when frailty follows a gamma distribution [16]. Second, theoretical results indicate that the gamma distribution is the limiting distribution of the frailty of long‐time survivors, irrespective of the frailty distribution at baseline [20]. However, estimators of the PPSH model may be biased if the functional form of the frailty distribution is misspecified. Nonetheless, empirical evidence suggests that this bias is generally minimal. Simulation studies indicate that the gamma distribution is robust to such misspecification in terms of bias and efficiency. Specifically, when the true frailty distribution is either inverse Gaussian or positive stable, the estimates from a Cox model with an assumed gamma frailty distribution are minimally affected [21, 22, 23].

In this sensitivity analysis, we generate datasets assuming frailty follows an inverse Gaussian distribution at baseline, with a mean of 1 and a variance of 1/γ. All other parameter settings, simulation algorithms, and estimation procedures, including assuming a gamma distribution for frailty in estimation, remain unchanged.

It is worth noting that, if the frailty used in data generating follows a gamma distribution, the nonfatal event generating algorithm described in Section 4.1 by theory will provide a PPSH, and when there is no mortality, it should provide a proportional hazard. However, if the frailty used in data generating follows an inverse Gaussian distribution, the data generating algorithm described in Section 4.1 does not, in theory, provide a marginally PPSH, and there is no explicit formula by which we can generate a marginally PPSH either. Nonetheless, the inverse Gaussian distribution shares a similar shape with the gamma distribution, so we could expect that even if we generate the datasets using the algorithm in Section 4.1 with inverse Gaussian distributed frailty, the proportionality assumption may not be severely violated. In fact, we use this algorithm to generate datasets from inverse Gaussian distributed frailty in the absence of mortality with other parameter values specified in Section 4.1 and test the proportional hazards assumption with the R function cox.zph. At the α=0.05 level, there is no violation of the proportional hazards assumption 84%, 93%, and 95% out of 10 000 simulations when γ=0.5, γ=2, and γ=5, respectively. As a comparison, if the gamma‐distributed frailty is used in data generating, 95% out of 10 000 simulations, there is no violation of the proportionality assumption across all γ, at the α=0.05 level. Therefore, with inverse Gaussian distributed frailty, although the proportionality assumption is not perfectly protected, we could still expect the estimated PS hazard ratio to provide an informative summary of the treatment effects.

Table 2 presents results from the sensitivity analysis. Compared to Table 1, when the population is highly heterogeneous (γ=0.5), the average proportion of deaths and subjects with a nonfatal event is lower when the frailty follows a gamma distribution. This discrepancy arises because, with the same mean and variance, the first quartile of the gamma distribution is generally smaller than that of the inverse Gaussian distribution. This suggests that the relatively healthier subjects in the gamma‐generated population have a lower probability of death and nonfatal events compared to those in the inverse Gaussian‐generated population.

Estimation results from the hypothetical dataset are presented as a reference where ν0=ν1=0. We set r=0.5 in the simulation, but since the frailty now follows an inverse Gaussian distribution, logr is no longer the theoretical true value of β. “Est” is the Monte Carlo sample mean of estimates from the hypothetical dataset, serving as the empirical true value of β. Then the bias is defined as the Monte Carlo sample mean of estimates minus the corresponding empirical true value. In a homogeneous population (large γ), the difference in frailty distribution has minimal impact. The empirical HR from the hypothetical dataset is close to 0.5, and the bias from both CS and PS approaches is minimal. As population heterogeneity increases, the empirical HR deviates further from 0.5, and the bias of the CS approach becomes pronounced. The bias of the PS is minimal when γ˜=γ, and is in general smaller than that from the CS approach when γ˜ is misspecified. Similar to the results in Table 1, when the population is heterogeneous, the PS hazard ratio approaches the CS hazard ratio as γ˜ increases, for the reason discussed in Section 4.2.

A Real Example

5

We illustrate our method using data from the Carvedilol Prospective Randomized Cumulative Survival (COPERNICUS) trial [24], a double‐blind, placebo‐controlled study examining the impact of carvedilol on morbidity and mortality among patients with severe heart failure. The trial included 2289 patients in total over a mean period of 10.4 months. Primary outcome analysis showed a 35% decrease in all‐cause mortality with carvedilol. 700 hospitalizations occurred in the carvedilol arm and 848 in the placebo arm. Notably, 38% had no hospitalizations, with 25%, 13%, 10%, and 14% experiencing one, two, three, and more than three hospitalizations, respectively. After excluding hospitalizations after the first, 382 hospitalizations remained in the carvedilol arm and 429 in the placebo arm. Compared to trials with mixed‐risk populations [25], the COPERNICUS, focusing on high‐risk patients, has a higher proportion experiencing multiple hospitalizations.

Table 3 presents a comparison of estimation results between the PS and the CS approach. γ˜ represents the value of γ employed in the PS approach, indicating the assumed heterogeneity level within the study population. HR of the PS approach is calculated as exp(β^). 95% CI of the PS approach is the percentile interval from 1000 bootstraps: (exp(β0.025∗),exp(β0.975∗)), where β(0.025)∗ and β(0.975)∗ represent the 2.5th and 97.5th percentiles of the bootstrap estimates.

The PS hazard ratio is smaller than the CS hazard ratio regardless of γ˜. Carvedilol generally extended survival, implying that subjects surviving to time t in the placebo arm were generally healthier at baseline compared to those in the carvedilol arm, and hence expected to be at lower risk of the nonfatal event. The CS hazard ratio is thus closer to 1.

In this example, the PS hazard ratio is not very sensitive to the choice of γ˜, ranging from 0.781 in an assumed highly heterogeneous population (γ˜=0.25) to 0.820 in an assumed highly homogeneous population (γ˜=10). The PS hazard ratio approaches the CS hazard ratio as γ˜ increases, because the CS approach essentially treats the population as homogeneous, which corresponds to γ˜→∞.

The CI from the PS approach narrows with larger γ˜ but remains wider than that from the CS approach due to the variability introduced in estimating PS probabilities. Using only the first nonfatal event from each subject might seem to diminish the power compared to approaches utilizing all nonfatal events. However, this is not necessarily true, as the PS approach is based on PL, which yields optimal power. Recurrent event approaches that focus on the marginal feature rely on estimating equations [26], which exhibit lower efficiency unless the variance structure is correctly specified. Sun and Cook [27] analyzed COPERNICUS data using all nonfatal events but reported wider CIs. The lack of power may be because risk ratios were derived from estimating equations, which are less efficient than likelihood‐based methods.

The proportionality test is conducted using the method proposed in the Appendix C with g(t)=t. The p values are reported in Table 3 and did not reveal violations of the proportionality assumption.

Discussion

6

In clinical trials involving both mortality and morbidity, the CS hazard ratio for the first nonfatal event remains commonly used [28]. However, it is crucial to acknowledge that this ratio does not reflect the direct effect on the underlying first nonfatal event process, unless the treatment has no differential effects on survival or the population is entirely homogeneous (within each given value of covariates). When treatment has differential effects on survival, the CS hazard ratio should be interpreted with caution. The difference between the CS hazard ratio and the direct effect grows with population heterogeneity, and cardiovascular studies reported significant variability among subjects [29]. Simonetto et al. [30] estimated inter‐subject heterogeneity in coronary heart disease risk using mortality data, suggesting that frailty with a variance between 1 and 4 is appropriate to capture total heterogeneity. In our setting, this corresponds to γ˜ between 0.25 and 1. While not directly applicable, their findings offer insights into a plausible range for γ˜ in cardiovascular trials. By specifying a range of γ, particularly small values, researchers can gauge the potential difference between the CS hazard ratio and the conditional direct effect. In the COPERNICUS trial, the CS hazard ratio is estimated at 0.82, slightly conservative but close to the PS hazard ratio. When utilizing the CS hazard ratio to characterize treatment effects on the first nonfatal event, providing the PS hazard ratio alongside can offer researchers a clearer understanding of the reliability of the CS hazard ratio and aid in understanding the direct effect on the underlying first nonfatal event process.

Over the past decade, it has been argued that for all‐cause mortality, the hazard ratio in the Cox model does not provide a causal interpretation when factors influencing the at‐risk process are not controlled for. Specifically, the Cox hazard ratio tends to underestimate the causal effect of a beneficial treatment because the active treatment arm often includes relatively more frail subjects compared to the placebo group at each time t>0 [15, 31, 32]. Although the Cox model is never correctly specified in practice (since we cannot control for all factors affecting the at‐risk process), it still serves as a useful approximation of the causal effect. Nonetheless, this paper does not focus on providing a single point estimate of the conditional direct effect. Instead, we aim to establish a plausible range for the conditional direct effect, illustrating how the CS hazard ratio may differ from the conditional direct effect. This difference arises from the dependence between death and the first nonfatal event, even when conditioned on covariates that might partially explain this dependency. To quantify this difference, we use the frailty approach, tuning the distribution of the frailty according to our beliefs about the extent of this dependency.

The practical application of SACE has faced criticism. For each subject i, let Mi=min(Yi1,Yi0). For j=0,…,m, if Mi>tj, then i∈𝒜(tj). Therefore, 𝒜(0)⊇𝒜(t1)⊇⋯⊇𝒜(tm), where 𝒜(0) represents the entire population, as all subjects can survive at t=0. Since Mi is not observable due to its counterfactual nature, it has been argued that the subpopulation 𝒜(tj), for which the SACE at tj is relevant, is an unidentifiable subset of the population. This argument holds if we consider the SACE at a specific single tj. However, under the assumption of proportionality, all subjects contribute to the estimation of the PS hazard ratio, unless the subject who experiences a nonfatal event at t1 would not survive to t1 in the counterfactual arm. Therefore, every subject may contribute to the estimation to varying degrees. Additionally, some argue that identifying the SACE relies on untestable assumptions. Despite these criticisms, the value of SACE lies in its role as the only well‐defined estimator that measures the direct effect on the underlying first nonfatal event process in the presence of the competing risk of death.

The last known follow‐up time for each subject may contribute to estimating the PS probabilities. One approach is to use the frailty θ to capture the pairwise correlation between the first nonfatal event time and potential death times from each treatment arm. This approach implicitly assumes that the correlation between death times is the same as the correlation between the death time and the first nonfatal event time. However, the correlation between death times across treatment arms is likely stronger than the correlation between death and the nonfatal event within the same arm. Intuitively, given a death time Yi from the placebo arm, we expect the potential death time in the active treatment arm to be close to F1−1{F0(Yi)}, where F1 and F0 represent the CDF of Y. In contrast, uncertainty remains regarding when the first nonfatal event occurs. To accommodate varying degrees of correlations when estimating PS probabilities, one possible approach is to employ a trivariate copula model with two dependence parameters α0 and α1 [33], where α1 measures dependence between Yz and Y1−z and α0 measures dependence between Tz and Yz as well as between Tz and Y1−z. Appendix E details the estimation of PS probabilities using trivariate copulas and presents the corresponding estimation results for the COPERNICUS trial. Regardless of the method chosen to estimate the PS probabilities, the proposed structure of the PPSH model remains applicable.

In practice, baseline covariates X are typically collected. X can be incorporated into the estimation of PS probabilities. In step 1 outlined in Section 3.1, baseline covariates can be used to estimate ŜY(t|Z=z,X) and ŜT(t|Y>t,Z=z,X), allowing us to compute ηYz(t|X) and ηTz(t|X) in step 2. However, the practical benefit of adjusting for covariates in this context is limited. Including baseline covariates may appear to relax our conditional independence assumptions to Y1−z╨Tz|X,θ and Yz╨Tz|X,θ, implying that X accounts for some of the correlation between potential outcomes. Nevertheless, we prespecify a range of γ values, and the correlation strength implied by this range generally covers what can be explained by X. Thus, covariate adjustment in the estimation of PS probabilities may be redundant. When fitting the PPSH model, covariates X can be included in Equation (2) under a multiplicative effect assumption, similar to the Cox model. In this case, the interpretation of the PS HR shifts from marginal to conditional due to non‐collapsibility. Appendix E presents the estimated PS hazard ratios, conditional on baseline covariates, for the COPERNICUS trial.

The PPSH model is designed to address the time to the first nonfatal event, wherein subjects who experience the nonfatal event are no longer at risk for that event. However, there is a possibility of expanding our framework to recurrent event scenarios, similar to the Anderson‐Gill model, which serves as an extension of the Cox model.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data files: sim70311‐sup‐0001‐DataFiles.zip.