Matched design for marginal causal effect on restricted mean survival time in observational studies

1.1 Causal inference for observational survival data

In biomedical studies, time-to-event is a commonly used outcome measure, and the statistical analyses for such data are usually referred to as survival analysis. Investigating the causal relationship between exposure and the time-to-event outcome is an important topic, with either randomized trials or observational studies. Causal inference for observational survival data has several challenges. First, since not all subjects can be observed for the full duration of time to event, the survival data suffer from censoring, which is a type of missing data problem. Therefore, standard statistical methods are usually not sufficient to handle both censoring and the missingness of potential outcomes. Second, the hazard ratio is a popular choice for measuring the association of survival outcomes between two groups, for convenience and easy interpretation. However, the hazard ratio is generally not an appropriate marginal causal effect measure due to its noncollapsibility property [1–3]. Other effect measures need to be considered to warrant valid marginal causal interpretation for survival data. Third, confounding is a major challenge in observational studies, which includes both measured and unmeasured confounders. Propensity score adjustments are popular tools for controlling the observed confounding [4]. But even with successful adjustment of observed confounding, observational data are still vulnerable to unmeasured confounding. Thus, appropriate sensitivity analysis needs to be developed to assess the impact of hidden bias [5].

The issues of using the hazard ratio as a marginal causal effect measure have been discussed extensively in the literature. Greenland et al. [1] pointed out that the hazard ratio has the noncollapsibility property when the treatment effect is nonzero. Hernán [6] argued that using the hazard ratio as a treatment effect measure may not have valid causal interpretation even in randomized studies, since the hazard ratio has a built-in selection bias and may change over time. Martinussen and Vansteelandt [2] studied the estimation of the treatment effect in the presence of confounders and found that the amount of confounding due to noncollapsibility in the Cox proportional hazards (PH) model would be very difficult to quantify. Aalen et al. [7] offered a more theoretical perspective on the conditions under which the hazard has a valid causal interpretation. They suggested that the hazard function h ( t , x , z ) must satisfy an additive assumption h ( t , x , z ) = a ( t , z ) + b ( t , x ) to yield a causal interpretation, where a ( t , z ) is a function of survival time t and treatment assignment z and b ( t , x ) is a function of survival time t and covariates x . Ni et al. [8] further illustrated that even under a PH model, the marginal hazard ratio is not a constant, after integrating out covariates. Thus, a valid and simple-to-use causal effect measure for survival outcomes is highly desirable.

1.2 RMST difference as a marginal causal effect measure

The restricted mean survival time (RMST) has been used in randomized clinical studies to evaluate treatment effects [9,10]. The RMST difference is more advantageous than the hazard ratio as a marginal effect measure. First, the RMST has an intuitive interpretation as the area under the survival curve over a certain time horizon. Second, the RMST difference is the difference of truncated mean survival time between two groups, which is essentially a mean difference. So it is collapsible, meaning that the marginal and conditional effects are compatible. Third, the treatment effect measured by the RMST difference can be asymptotically unbiasedly estimated without PH assumption, while the conventional Cox model heavily relies on such assumption.

To take advantage of the collapsibility of RMST difference, we can construct RMST regression by including covariates to control for confounding or increase estimation efficiency. Several methods of regressing RMST on multiple covariates have been developed. Karrison [11] examined the RMST as an index for comparing survival in two groups and proposed to model the hazard with piece-wise exponential models assuming covariates have a multiplicative effect on the hazard. Zucker [12] further simplified the implementation procedure for Karrison’s method and provided an extended version to achieve robustness against model misspecification. Andersen et al. [13] compared several regression analysis methods of mean survival time and RMST, and they proposed a regression method based on pseudo-observations. Tian et al. [14] developed an RMST regression model with adjustment for baseline covariates. They constructed an estimating equation with the inverse probability of censoring weighting (IPCW) to obtain consistent estimates. Wang and Schaubel [15] modeled the RMST using generalized estimating equation methods, which allows censoring to depend on both baseline covariates and time-dependent factors.

Although RMST differences have been reported in many randomized clinical studies, there is only limited discussion of using RMST in observational studies, probably due to the challenge of confounding adjustment. Propensity score weighting and stratification methods have been explored in the literature. Zhang and Schaubel [16] derived a double-robust estimator for RMST difference based on the inverse probability of treatment weighted (IPTW) estimating equation with augmentation term. To adjust for confounding factors, they built three working models for survival time, treatment assignment, and censoring, then incorporated them into the augmentation term. They assumed the PH assumption in outcome modeling, which might be violated in practice. Conner et al. [17] proposed a weighted method to compare the adjusted RMST difference directly. Unlike Zhang and Schaubel’s work, Conner et al. estimated the RMST based on the Kaplan–Meier (KM) estimator rather than the Nelson–Aalen estimator. They adjusted the KM estimator with IPTW and derived the adjusted RMST by integrating the IPTW-adjusted KM estimator. Ni et al. [8] proposed a propensity score stratified RMST difference estimation strategy to examine the marginal causal effect with observational survival data, which can combine stratification with further regression adjustment. Some existing methods still rely on modeling assumptions, and no matching based method has been proposed. More importantly, none of them touches on unmeasured confounding assessment, which is a big issue in observational studies. The matched design provides more flexibility in confounding adjustment. Observed confounding can be controlled via matching and additional regression modeling. Unobserved confounding can be explored via sensitivity analysis in matched datasets. To address the void of literature, we propose a matching-based RMST difference estimation strategy, which also facilitates the implementation of sensitivity analysis for hidden bias.

1.3 A motivating example: atherosclerosis risk in communities (ARIC)

In the United States, stroke is a severe disease that causes serious disability for adults and is a leading cause of death [18,19]. Several previous studies have shown that smoking is an important risk factor for stroke [20,21], and even passive smoking could increase the risk of stroke [22]. Although the causal pathway between smoking and stroke is unclear, Shah and Cole [23] found that the more people smoke, the more likely they were to have a stroke, and people who quit smoking showed a significantly lower risk of stroke, which provides some evidence for the causal relationship between smoking and stroke.

The ARIC study [24] is a prospective cohort study conducted in four U.S. communities. Four thousand adults aged 45–64 years were randomly sampled from each of four U.S. communities, and the final dataset contains information of 15,792 individuals. After a baseline examination during 1987–1989, subjects were followed up for the development of incident ischemic stroke, and the first definite or probable hospitalized stroke. Due to the length of follow-up, not all event times were observed, so the data were subjected to censoring. One primary outcome is the time to first stroke or death (whichever comes first), and a subject is censored if the incidence of stroke or death is not observed by the end of the study. We try to answer a causal question, using this ARIC dataset: how smokers’ stroke-free survival would change had they not smoked at baseline. Matched design is a natural choice to address this as it is about the causal effect of those being exposed to smoking, rather than for the entire population.

Existing literature mostly used the Cox PH model to analyze the ARIC data. Kwon et al. [25] and Ding et al. [26] studied the association between smoking status and risk of stroke, using the Cox PH model to estimate the hazard ratio (HR) of smoking status on the risk of stroke. Thus, the estimated effects were interpreted as conditional rather than marginal. Moreover, it is possible that some prognostic factors were not included in the confounder adjusted regression model, which would lead to biased estimates of conditional effects. An analysis with RMST as the effect measure may provide new insight into this research.

In this article, we propose a propensity score matching-based RMST difference estimator and develop a corresponding sensitivity analysis strategy for assessing the impact due to unmeasured confounding. We apply this method to the ARIC study to examine the causal effect of smoking on stroke-free survival. The rest of this article is organized as follows: In Section 2, we set up the notation and assumptions and describe the proposed RMST estimator with its theoretical properties. In Section 3, we conduct a simulation study to examine the empirical performance of our proposed method under different scenarios and also compare it with several commonly used methods in practice. In Section 4, we develop a sensitivity analysis strategy by adapting the E -value approach to matched data. In Section 5, we present the analysis results of the ARIC data. Section 6 concludes this article with some discussions.

2.1 Notation and assumptions

We follow the potential outcomes framework proposed by Rubin [27] to define the causal effects. In a two-arm survival analysis study, let A be the treatment assignment indicator (or more generally, the exposure status), such that A = 1 indicates being exposed to the treatment and A = 0 indicates being exposed to the control. Let T a denote the potential event time and S a ( t ) denote the corresponding survival function for a subject if under treatment value a . The following two assumptions are extensions of commonly used assumptions for causal inference in observational studies [28].

The potential restricted event time is defined as Z a = min ( T a , τ ) , where τ is the truncation time point, which is usually prespecified at the design stage based on clinical relevance and study feasibility. Both T a and Z a are subject to censoring by a random variable C . We introduce two additional assumptions for survival data.

This assumption ensures that we can asymptotically unbiasedly estimate the survival function via the KM approach within the matched sample. This also implies the conditional independence between censoring and the truncated survival times, Z a .

Assumption 4 is a technical one to ensure that the prespecified τ is clinically meaningful, and RMST can be asymptotically unbiasedly estimated.

Let δ a = I ( Z a < C ) denote the censoring indicator, and then the observed restricted time is defined as Y a = min ( Z a , C ) = ( Z a ) δ a C ( 1 − δ a ) . For a subject under treatment value a , the potential outcome of RMST is defined as μ a ( τ ) = E ( Z a ) = ∫ 0 τ S a ( t ) d t , then the average treatment effect (ATE) on RMST, denoted by Δ A T E , can be defined as follows:

Then the average treatment effect for the treated (ATT) on RMST, denoted by Δ ATT , can be defined as follows:

which is more meaningful for many observational study applications, where only a portion of the population, not everyone, could have been exposed to the treatment. Since Z a = min ( T a , τ ) and τ is a fixed constant, ( T 0 , T 1 ) ⊥ ⊥ A ∣ X implies ( Z 0 , Z 1 ) ⊥ ⊥ A ∣ X , and similar conclusions could be made about assumption 3. Following Theorem 3 in Rosenbaum and Rubin [4], we can establish the strong ignorability based on propensity score e ( X ) = P ( A = 1 ∣ X ) for survival outcomes in proposition 1 (proof provided in Appendix A).

2.2 Matched RMST difference estimator

In randomized trials, the marginal causal effect of the treatment on RMST can be asymptotically unbiasedly estimated [29] by direct contrast of group-specific RMST estimates since confounding effects are eliminated by design. In observational studies, however, additional adjustments are needed for confounding control. Propensity score-based approaches are popular for this purpose, which may take the form of matching, stratification, or weighting [4,30]. Among different propensity score adjustment strategies, matching is a design tool that selects comparable control units to match with treated units, and it often results in more robust causal effect estimates as it does not rely on outcome model specification. Usually, matching uses all treated and a subset of control units, so it estimates the ATT [28].

Our proposed propensity score matched RMST estimation includes the following steps:

1. Propensity score estimation. The propensity score is defined as the conditional probability of treatment given a vector of observed covariates [4]. We estimate the propensity score by fitting a logistic regression on A with X , though other estimation options, either parametric or nonparametric, are also available [31, 32].

2. Propensity score matching. We use the optimal matching algorithm by Hansen and Klopfer [33] to create pair matches without replacement based on the estimated propensity score, and the unmatched controls will be removed from the matched sample. Matching quality is assessed by checking the postmatching covariate balance. Any substantial covariate imbalance would lead to a recalibration of the propensity score model. We will proceed to the next step only after a satisfactory balance is achieved. Note that matching is used as a design procedure, and the specific propensity score values are not used in the subsequent estimation process (e.g., not a part of any estimating equations). So the uncertainty of propensity score estimation is not considered in subsequent analysis and variance calculation.

3. Treatment effect estimation. Suppose we obtain n pairs of data through matching, where each pair contains exactly one treated and one control subject. We estimate the RMST, μ ( τ ) , by μ ˆ ( τ ) = ∫ 0 τ S ˆ ( t ) d ( t ) , where S ˆ ( t ) is estimated by the nonparametric KM method. Let S ˆ 0 ( t ) and S ˆ 1 ( t ) denote the KM estimates of survival function for control and treated groups in the matched sample, respectively. On the basis of the matched sample, our estimator for the averaged treatment effect on the treated (ATT) is

   Δ ˆ ATT = μ ˆ 1 ( τ ) − μ ˆ 0 ( τ ) = ∫ 0 τ [ S ˆ 1 ( t ) − S ˆ 0 ( t ) ] d t .

2.3 Variance estimation

The matching process may introduce correlations between the two subjects in the same pair, as they are matched on similar propensity scores. Therefore, the variance calculation of Δ ˆ ATT needs to account for such correlation:

The overall variance has two components, the marginal variance of RMST estimates and their covariance. For two dependent event times with independent censoring and no competing risk, Murray and Cole [34] provided closed-form asymptotic covariance formulas for KM survival estimates and corresponding RMST estimates. To address the dependence structure introduced in the matching process, we adapt their formulas to compute the covariance between the control and treated group RMST estimates in the matched sample.

Specifically, let T 0 be the event time for a subject from the control group with marginal hazard function h 0 ( ⋅ ) , and T 1 be the event time for a subject from the treatment group with marginal hazard function h 1 ( ⋅ ) , then the event times for a matched pair of control and treated subject can be denoted as ( T 0 , T 1 ) . Let C 0 and C 1 be the censoring variables for the control and treated subject, respectively. Then the observed time can be denoted as T ˜ 0 = min ( T 0 , C 0 ) for control group with censoring indicator δ 0 = I ( T 0 < C 0 ) and T ˜ 1 = min ( T 1 , C 1 ) for treated group with censoring indicator δ 1 = I ( T 1 < C 1 ) . Then, the joint hazard function is h i j ( u , v ) = lim Δ u , Δ v → 0 1 Δ u Δ v P ( u ≤ T ˜ i < u + Δ u , v ≤ T ˜ j < v + Δ v , δ i = 1 , δ j = 1 ∣ T ˜ i ≥ u , T ˜ j ≥ v ) where i , j ∈ { 0 , 1 } , and the conditional hazard function is h i ∣ j ( u ∣ v ) = lim Δ u → 0 1 Δ u P ( u ≤ T ˜ i < u + Δ u , δ i = 1 ∣ T ˜ i ≥ u , T ˜ j ≥ v ) , where i , j ∈ { 0 , 1 } . Then, the covariance between two RMSTs can be computed as follows:

where G 01 ( u , v ) = P ( T ˜ 0 ≥ u , T ˜ 1 ≥ v ) P ( T ˜ 0 ≥ u ) P ( T ˜ 1 ≥ v ) [ h 01 ( u , v ) − h 0 ∣ 1 ( u ∣ v ) h 1 ( v ) − h 1 ∣ 0 ( v ∣ u ) h 0 ( u ) + h 0 ( u ) h 1 ( v ) ] . Details about the computation of function G 01 ( u , v ) are included in Appendix C.

For the marginal variances, two methods may be considered:

1. Murray’s method: The aforementioned covariance formulas can be used to compute the marginal variance, since the marginal variance of RMST could be written as the covariance with itself, that is, var ∫ 0 τ S ˆ ( t ) d t = cov ∫ 0 τ S ˆ ( t ) d t , ∫ 0 τ S ˆ ( t ) d t .

2. Hosmer’s method: we may also consider the computation method introduced in the study by Hosmer et al. [35]. Let t 1 < t 2 < ⋯ < t D represent distinct event times. For each k = 1 , … , D , let Y k be the number of surviving units just prior to event time t k , and let d k be the number of events at t k . Let S ˆ ( t k ) = ∏ l = 1 k 1 − d l Y l denotes the KM estimate of the survival function at event time t k , and let N τ be the number of t k values that are less than truncation time point τ , and then the RMST is estimated by

   ∫ 0 τ S ˆ ( t ) d t = ∑ k = 1 N τ S ˆ ( t k − 1 ) ( t k − t k − 1 ) + S ˆ ( t N τ ) ( τ − t N τ ) ,

   and the marginal variance of RMST can be estimated as follows:

   var ∫ 0 τ S ˆ ( t ) d t = m m − 1 ∑ k = 1 N τ d k A k 2 Y k ( Y k − d k ) ,

   where A k = ∫ t k τ S ˆ ( t ) d t = ∑ l = k N τ S ˆ ( t l ) ( t l + 1 − t l ) + S ˆ ( t N τ ) ( τ − t N τ ) and m = ∑ l = 1 N τ d l .

3.1 Data generation

To assess the empirical performance of the proposed method, we simulate an observational dataset with known confounders. Several existing methods for causal inference with survival outcomes are compared.

We generate ten independent baseline covariates denoted by X 1 to X 10 . Among them, X 1 , X 3 , … , X 9 are five binary covariates following Bernoulli distribution with parameters 0.2, 0.4, 0.6, 0.8, and 0.5, respectively, and X 2 , X 4 , … , X 10 are five continuous covariates following the standard normal distribution. We then generate potential survival time T 1 as the outcome under treatment and potential survival time T 0 as the outcome under control from Weibull distribution [36]. Specifically, we simulate a uniform random variable Q on [0,1] and then generate the potential survival time as follows ( j = 1 , 0 ):

where A is the treatment indicator and β A is the conditional multiplicative treatment effect on the hazard function given covariates, and ν j and λ 0 j are the shape and baseline scale parameters of Weibull distribution for treatment group j , respectively. When ν 0 = ν 1 , we have the PH model, otherwise the model is nonproportional hazards (NP). The treatment indicator A is generated from Bernoulli distribution with P ( A = 1 ∣ X ) defined by the logistic model logit ( P ( A = 1 ∣ X ) ) = − 1.95 + log ( 1.2 ) X 1 + log ( 1.1 ) X 2 + log ( 1.4 ) X 3 + log ( 1.2 ) X 4 + log ( 1.6 ) X 5 + log ( 1.3 ) X 6 + log ( 1.8 ) X 7 . Thus, X 1 , X 4 , X 6 , and X 7 are true confounders. This setup allows about 20% of the population to be exposed to treatment.

In the simulation, we assume censoring variable C is marginally independent of T a and X over treatment indicator A for simplicity, and C is generated from an exponential distribution with rate parameter γ , which is chosen to create four different levels of censoring. For simplicity, we use the same censoring variable for both arms in the simulation.

Let τ be the prespecified truncation time point, and the observed event time is T = T 0 ( 1 − A ) + T 1 A . We generate the restricted event time Z = min ( T , τ ) and the observed restricted time Y = min ( Z , C ) = min ( T , C , τ ) . The restricted event time Z is censored if the observed time C < Z with censoring status δ Z = I ( Z < C ) , otherwise it is noncensored.

We simulate 500 datasets of sample size 2,500 for each scenario and set the truncation time point τ to 100. The true RMST difference is determined by calculating the empirical difference between the potential RMSTs under treated and control conditions, and we compute both ATT and ATE versions of true RMST difference to serve as benchmarks for different methods as appropriate. In the j th simulated dataset, we calculate Δ j = ∑ i = 1 n Z i 1 − Z i 0 n , where Z i A = min ( T i A , τ ) is the potential restricted event time for the i th individual and n is the sample size of the treated group (for ATT) or the entire sample (for ATE). Then, the true marginal effect on RMST is calculated as Δ 0 = ∑ j = 1 500 Δ j 500 .

Both PH and NP settings are examined. Under both settings, we set β A to five different values: 0 , − 0.4 , − 0.8 , − 1.2 , and − 2 . For each treatment effect value, we also consider four different levels of censoring rates (CRs), which are 0, 20, 40, and 60%. Detailed parameter setup for PH and NP scenarios in observational studies are summarized in Table 1.

3.2 Estimation strategies

The proposed method is compared with three existing estimation strategies:

1. Propensity score matched RMST estimation. This is our proposed method as described in the previous section, and the propensity score is estimated using the correct model specification. The estimated treatment effect is compared to the ATT version of the true RMST difference in our simulation.

2. Conner’s IPTW RMST estimation. This method is proposed by Conner et al. [17], and they estimated the RMST based on the inverse probability treatment weighting (IPTW) adjusted KM estimator. In our simulation, we use the ATT version of weight to adjust for observed confounding, so it is compared to the true ATT RMST difference. The propensity score is estimated using the correct model specification.

3. Tian’s RMST regression. This method is proposed by Tian et al. [14], which uses the IPCW estimating equation with identity link function to estimate the treatment effect on RMST with adjustment for covariates. The estimated treatment effect is compared to ATE version of the true RMST difference. We consider four different outcome models in Tian’s RMST regression: (1) outcome model using the treatment indicator only; (2) outcome model using the true covariate set; (3) outcome model using all covariates; and (4) outcome model using a wrong covariate set. Due to space limitation, only the results of RMST regression with true covariates are summarized in the following section, which has the best performance among the four models. An important caveat is that the RMST regression model with the true covariate set does not represent the true outcome model since the data are generated based on a hazard model.

4. Inverse probability treatment weighting (IPTW) Cox regression. This method estimates β A . The propensity score is estimated using the correct model specification. We use the ATT weight to fit a weighted Cox regression model and regard β A as the truth to calculate the bias and coverage probabilities since there is no single value true marginal hazard ratio. We consider four different outcome models in the IPTW Cox regression : (1) outcome model using the treatment indicator only; (2) outcome model using the true covariate set; (3) outcome model using all covariates; (4) outcome model using a wrong covariate set. Due to space limitation, only the results of IPTW Cox regression with the true covariate set are summarized in the following section, which has the best performance among the four models. We understand that the results here are not directly comparable to the first three methods, as they are based on different effect measures. Due to the high popularity of the IPTW Cox model in practice, however, we think there is some value in presenting the results as a reference.

3.3 Performance assessment

We summarize treatment effect estimates from 500 Monte Carlo iterations into four measures: (1) percentage bias (Bias %), which is the bias divided by the true value for nonzero treatment effect scenarios. For the zero treatment effect scenario, we just report the bias. The bias is computed as the average of 500 treatment effect estimates minus the truth; (2) coverage probability (CP), which is the proportion of 500 95% confidence (CIs) intervals that cover the truth; (3) model-based standard error (SEM), which is the average of the 500 estimated standard errors from the model-based formula; and (4) empirical standard error (SEE), which is the standard error of the 500 point estimates of the treatment effect.

3.4 Results

Simulation results under the PH setting are summarized in Table 2. The proposed matched RMST method generates unbiased estimates of the target parameters under most scenarios, and the coverage probabilities are around 95%. For a small effect size ( β A = − 0.4 ), the bias is a bit large for a high CR. Conner’s method has a similar performance, with moderately larger biases. Averaging across all scenarios, bias from Conner’s method is 65% higher than our method. The results of the IPTW Cox model are mostly good since we use the correct outcome model. The CP may be a bit lower than the nominal level, sometimes, which may be due to the underestimated standard error. Tian’s RMST regression method shows a relatively large percentage bias and lower CP, especially under scenarios with large treatment effects. This is likely due to the incorrect covariate functional form specification in the model even though we include the right covariate set.

Simulation results under the NP setting are summarized in Table 3. Both our matched RMST method and Conner’s method have similar performance (with the latter having more bias) as under the PH setting since these methods do not rely on the PH assumption. Tian’s RMST regression method performs somewhat worse, with a bigger bias and much lower than ideal coverage probabilities. Because the PH assumption does not hold here, the IPTW Cox model completely misses the target with large bias and very small coverage probabilities.

4.1 An overview of E -value

Propensity score adjustment can only control for observed confounding. Unmeasured confounding is likely to be present in observational studies since researchers have no control over the treatment assignment. Thus, sensitivity analysis is important to assess the impact of hidden bias.

Ding and VanderWeele [37,38] developed a new sensitivity analysis strategy, known as the E -value method. It assumes a hypothetical unmeasured confounder, U , and provides a lower bound of the strength of association on the risk ratio scale that U would have both the exposure and the outcome, to explain away the observed association. Below is a brief review of the conventional E -value method to set the stage for our sensitivity analysis of RMST difference.

Let A denote a binary exposure and D denote a binary outcome, X is a vector of measured confounders, and U is a binary unmeasured confounder with levels k = 0 , 1 . The observed relative risk of exposure A on the outcome D within stratum of X = x is expressed as follows:

Then the relative risk of exposure on level k of the unmeasured confounder U within the stratum of X = x is expressed as follows:

Since U is not observed, to facilitate the analysis, we take the maximal relative risk of A on U within stratum X = x , denoted as RR A U ∣ x = max k RR A U , k ∣ x . Similarly, we can define an upper bound for the relative risk between U and D as RR U D ∣ x = max ( RR U D ∣ A = 0 , x , RR U D ∣ A = 1 , x ) , where RR U D ∣ A , x is an upper bound of the relative risk between U and D in exposed or unexposed group, within stratum X = x . If X and U are sufficient to control for all confounding effects, the true causal relative risk is

The relative risk pair ( RR A U ∣ x , RR U D ∣ x ) are used to measure the strength of confounding between the exposure A and the outcome D induced by the confounder U within the stratum of X = x . Even though we cannot estimate the true relative risk, its ratio with the observed relative risk is bounded by the following quantity, which is a function of the sensitivity parameters RR A U ∣ x and RR U D ∣ x :

For given values of RR A U ∣ x and RR U D ∣ x , we can identify a range of possible values for the true relative risk. If the range covers one, the observed significant association would be explained away by the presence of unmeasured confounding at the given magnitude.

4.2 Sensitivity analysis on RMST difference with matched data

This E -value method can be adapted to conduct sensitivity analysis for our RMST difference estimator in matched design. There are a series of propositions to justify the theoretical validity of using the E -value for the RMST difference estimator. In the interest of space, we just illustrate the main idea in this subsection and present the propositions and their detailed proofs in Appendix B.

Let A be the treatment indicator and Z = min ( T , τ ) be the RMST outcome, where T is the event time and τ is the truncation time point. Let e ( X ) be the propensity score and U be a binary unmeasured confounder with levels k = 0 , 1 . The relative risk of treatment A on level k of the unmeasured confounder U with a given propensity score value e ( X ) = e ( x ) is defined as follows:

The maximal relative risk of A on U with e ( X ) = e ( x ) is RR A U ∣ e ( x ) = max k RR A U , k ∣ e ( x ) . We define the expectations of the RMST outcome Z given U = u and e ( X ) = e ( x ) with and without treatment as follows:

Then, the mean ratios of U on Z with and without treatment with e ( X ) = e ( x ) are defined as follows:

As shown in proposition 4 in Appendix B, since both unmeasured confounder parameters RR A U ∣ e ( X ) = e ( x ) and MR U Z ∣ e ( X ) = e ( x ) are no less than 1, we can identify the bounding factor as follows:

where ( RR A U ∣ e ( X ) = e ( x ) , MR U Z ∣ e ( X ) = e ( x ) ) are prespecified sensitivity analysis parameters. In theory, one can identify a separate bounding factor for each matched pair and choose the maximal value to facilitate the calculation. But this could be quite cumbersome in practice. Instead, we follow the idea of using Γ in the conventional Rosenbaum’s sensitivity analysis, which is a prespecified upper bound of the association. Denote RR A U = max e ( x ) ( RR A U ∣ e ( X ) = e ( x ) ) and MR U Z = max e ( x ) ( MR U Z ∣ e ( X ) = e ( x ) ) , then the maximal bounding factor can be calculated as follows:

Because B F U ∗ is an increasing function of both RR A U and MR U Z , taking prespecified upper bounds of both associations leads to an upper bound of bounding factors. In practice, one can identify a range of possible values for ( RR A U , MR U Z ) and calculate the upper bound of treatment effects for each combination using the following formulas:

Let ACE A Z true denote the true average causal effect. When the treatment effect is positive, we have

When the treatment effect is negative, we have

The sensitivity analysis is generally done by checking the behavior of the 95% CI bounds. Since our real data analysis shows a negative treatment effect, we focus on equation (4) and use a one-sided 95% CI for illustration purposes. Equation (4) implies that the treatment effect estimate should be bounded by a function of the bounding factor and the mean survival times from each treatment group. Denote the right-hand side quantity as right-hand side (RHS). We can calculate the 95% CI for RHS using normal approximation and denote the upper confidence bound as RHS u b . Because ACE A Z true ≤ RHS , we have P ( ACE A Z true ≤ RHS u b ) ≥ 0.95 . Therefore, we can regard RHS u b as an upper bound for the upper 95% confidence bound of the treatment effect estimate. If this value is less than zero, we would reject the null hypothesis. Otherwise, we would retain the null.

4.3 Interpreting the sensitivity analysis

For an unmeasured confounding with a prespecified magnitude of associations ( RR A U , MR U Z ) , the bounding factor B F U ∗ can be computed by equation (2). For a positive treatment effect based on the observed data, we can compute a lower bound for the lower 95% confidence bound of the treatment effect estimate by equation (3). A positive lower bound indicates that there is still a positive treatment effect with an unmeasured confounding effect of magnitude ( RR A U , MR U Z ) . A nonpositive lower bound indicates that the positive treatment effect could be explained away by the unmeasured confounding of magnitude ( RR A U , MR U Z ) . For a negative treatment effect based on the observed data, we can compute an upper bound for the upper 95% confidence bound of the treatment effect estimate by equation (4) as described earlier, and similar interpretations can be made. A negative upper bound indicates that there is still a negative treatment effect with an unmeasured confounding effect of magnitude ( RR A U , MR U Z ) . A nonnegative upper bound indicates that the negative treatment effect could be explained away by the unmeasured confounding of magnitude ( RR A U , MR U Z ) . A detailed numerical example is presented in Section 5.