Estimating Mann–Whitney-Type Causal Effects for Right-Censored Survival Outcomes

2.1 Notation and assumptions

Let A denote the actual treatment received by an individual subject in a study, which may be a randomized clinical trial or an observational study. Assuming consistency or stable unit treatment value, the actual survival time is then T=T(A)=AT(1)+(1−A)T(0). Let W be a vector of relevant covariates measured at baseline (before treatment). We assume that treatment assignment is ignorable [22] given W in the sense that

The ignorability assumption is trivially true in a randomized clinical trial. In an observational study, the assumption requires that W contain enough information to explain any association between A and the potential outcomes. We also assume positivity:

which is trivially true in a clinical trial but not trivial in an observational study.

The actual survival time T may be right-censored by a censoring time C. We assume independent censoring:

where ⊥ denotes independence. The observed data for an individual subject can be represented as O=(W,A,X,Δ), where X=T∧C and Δ=I(T≤C). The observed data for the whole study will be conceptualized as independent copies of O and denoted by Oi=(Wi,Ai,Xi,Δi), i=1,…,n.

Our goal is to use the observed data to estimate θ for a general (specified) function h. For identifiability, we assume that h(t1,t0) is determined by t1∧τ and t0∧τ for some τ>0 such that P(C>τ|A,W)>0 almost surely. Accordingly, we will restrict attention to the interval (0,τ] in discussions of distribution, survival and hazard functions. For theoretical reasons, we also assume that E{h(Ti(1),Tj(0))2}<∞, which is satisfied by all examples mentioned earlier. The efficient influence function in this estimation problem is given in Appendix A. In the rest of this section, we propose several estimators of θ, whose asymptotic properties are studied in Appendix B. Asymptotic variance formulas can be derived but may be cumbersome to use in variance estimation, depending on the specific forms of the working models involved. For ease of implementation, we use bootstrap standard errors and confidence intervals for inference.

2.2 Outcome regression (OR) estimator

The OR method is based on the fact that

where i≠j and Fa(t|W)=P{T(a)≤t|W} for a=0,1. Here and in the sequel, we write h(ν1,ν0)=∬h(t1,t0)dν1(t1)dν0(t0), h(ν1,t0)=∫h(t1,t0)dν1(t1) and h(t1,ν0)=∫h(t1,t0)dν0(t0) for any measures (ν1,ν0). Assumption (3) implies that Fa(t|W)=P(T≤t|A=a,W), which, by assumptions (4) and (5), can be identified from the observed data.

To deal with the curse of dimensionality, we will assume a parametric or semiparametric model for Fa(t|W), say Fa(t|W;α), where α may be finite- or infinite-dimensional. This will be referred to as the OR model. A prominent example of a semiparametric model for Fa(t|W) is the proportional hazards model [23], [24]:

where λ(·|A,W;α) is the conditional hazard function of T given (A,W), α=(η,λ0), λ0 is the baseline hazard function, and V is a vector-valued function of (A,W). Under this model,

where Λ and Λ0 are cumulative versions of λ and λ0, respectively.

Whether Fa(t|W;α) is parametric or semiparametric, it is usually convenient to work with the corresponding hazard function λ(t|A,W;α). The likelihood for α based on the observed data may be written as

At least for parametric models and the proportional hazards model, α can be estimated by maximizing the above likelihood. Let αˆ be an estimate of α, which may be obtained by maximum likelihood or other means. Then θ can be estimated by

If the model Fa(t|W;α) is correct and αˆ is consistent for α (in a suitable sense), then θˆor is consistent for θ under mild regularity conditions. In Appendix B, we show that n(θˆor−θ) converges to a zero-mean normal distribution under general conditions. A key condition we assume is that αˆ is n-consistent and asymptotically linear in a suitable sense.

2.3 Inverse probability weighted (IPW) estimator

It follows from assumptions (3)–(5) and a conditioning argument that

where Sc(t|A,W)=P(C>t|A,W). If Sc(t|A,W) and p(W) are known, then the above identity suggests that θ can be estimated by

In reality, although p(W) is known by design in a clinical trial, Sc(t|A,W) is generally unknown and must be estimated. We assume a model for Sc(t|A,W), say Sc(t|A,W;β), which may be parametric or semiparametric so β may be finite- or infinite-dimensional. The model Sc(t|A,W;β) may be specified using the same considerations for specifying Fa(t|W;α) except that we are now modeling the censoring time instead of the survival time. Let βˆ be an estimate of β, which may be obtained by maximizing the likelihood:

where λc(·|A,W;β) is the conditional hazard function of C given (A,W) under the model Sc(t|A,W;β), and Λc(·|A,W;β) is the cumulative version of λc(·|A,W;β).

In an observational study, we also need to estimate p(W). In fact, even if p(W) is known, using an estimate of p(W) instead of the known value usually leads to better efficiency in the IPW approach [25]. Let p(W) be modeled as p(W;γ), where γ is a finite- or infinite-dimensional parameter. Because A is binary, a typical choice for p(W;γ) would be a logistic regression model. Let γˆ be an estimate of γ, which may be obtained by maximizing the likelihood:

Once βˆ and γˆ are obtained, the IPW estimator of θ is readily available as a weighted average:

The denominator here serves to normalize the inverse probability weights. If the models Sc(t|A,W;β) and p(W;γ) are correct and (βˆ,γˆ) estimate (β,γ) consistently, then θˆipw is consistent for θ and asymptotically linear under appropriate conditions (see Appendix B).

2.4 First doubly robust (DR1) estimator

The DR1 method is adapted from the DR method of Zhang and Schaubel [19] for estimating the difference in mean restricted survival time. As an important by-product, Zhang and Schaubel [19] propose a DR estimator of Fa(t), which can be described as follows. Let us define

Then the Zhang-Schaubel estimator of Λa(t), the cumulative hazard function of T(a), is given by

and the corresponding estimator of Fa(t) is Fˆadr(t)=1−exp{−Λˆadr(t)}. The superscript “dr” in these estimators indicates that the estimators are doubly robust in the sense of being consistent and asymptotically normal if (i) Fa(t|W;α) is correctly specified, (ii) Sc(t|A,W;β) and p(W;γ) are both correctly specified, or (iii) both (i) and (ii) hold.

Now θ can be estimated by θˆdr1=h(Fˆ1dr,Fˆ0dr), and it is easy to see that θˆdr1 is doubly robust in the same sense described above. Because the Fˆadr (a=0,1) are not guaranteed to be probability measures, θˆdr1 is not guaranteed to lie in the range of h. When θˆdr1 falls outside the range of h, it can be truncated into the range of h without changing its consistency or asymptotic normality. It is not immediately clear whether θˆdr1 is locally efficient, that is, whether it attains the nonparametric information bound when all three working models are correct. This is the main motivation for our developing a second DR estimator of θ based directly on the efficient influence function for estimating θ.

2.5 Second doubly robust (DR2) estimator

In Appendix A, we show that the efficient influence function for estimating θ is

where ℏ1(t,F0|W)=E{h(T,F0)|A=1,W,T≥t}, ℏ0(F1,t|W)=E{h(F1,T)|A=0,W,T≥t}, and Mc(t)=(1−Δ)I(X≤t)−Λc(X∧t|A,W). Motivated by this result, we propose to estimate θ by setting the sample average of an empirical version of ϕeff(O) to 0. The resulting estimator is

where

The use of Fˆadr in θˆdr2 to replace Fa in (6) is an attempt to provide maximal protection against model misspecification. In Appendix B, we show that θˆdr2 is indeed doubly robust and, furthermore, locally efficient. Like θˆdr1, θˆdr2 may fall outside the range of h but can be truncated back into the range of h.

3.1 Simulation

Here we report a simulation study of the finite-sample performance of the methods described in Section 2. Data are generated according to the following mechanism:

where 0=(0,0)′, I is the 2-by-2 identity matrix, expit(x)=1/{1+exp(−x)}, and ηA is either 0 (no treatment effect) or −0.5 (protective treatment effect). The function h in the definition of θ is take to be (2) with τ=2. The true value of θ is 0 when ηA=0 and approximately 0.15 when ηA=−0.5. We consider two sample sizes: n=200,500. In each case, 1,000 replicate samples are generated.

Each simulated sample is analyzed using the four methods (OR, IPW, DR1 and DR2) described in Section 2 as well as a naive method that ignores censoring and possible confounding and estimates θ with the U-statistic

where n1=∑i=1nAi and n0=n−n1. The OR, IPW, DR1 and DR2 methods are implemented with correct and incorrect working models. The correct models are:

where αa=(ηa,λa0), βa=(ηac,λa0c), and λa and λac are unspecified baseline hazard functions. The incorrect models result from replacing (W1,W2) with (I(W1>0),I(W2>0)) in the correct models. All models are estimated using the maximum likelihood approach.

Table 1 shows a summary of the simulation results: empirical bias and standard deviation for estimating θ. As expected, the naive method is severely biased. The OR estimator becomes biased when the model λ(t|A,W;α) is misspecified, as does the IPW estimator when the models λc(t|A,W;β) and p(W;γ) are misspecified. In contrast, the two DR estimators are nearly unbiased unless all models are misspecified, demonstrating double robustness. Regardless of model (in)correctness, the estimators that adjust for confounding and censoring generally follow the pattern OR<DR<IPW in terms of variability. The efficiency comparison of the two DR estimators does not seem to follow a clear pattern. At n=200, DR2 appears more efficient than DR1 when all models are correct, but this difference appears to diminish with increasing sample size.

3.2 Application

We now apply the methods compared in Section 3.1 to a study of hospitalized pneumonia in young children, described in Section 1.13 of Klein and Moeschberger [26]. The study is based on 3,470 newborn children from the National Longitudinal Survey of Youth [27]. The research question is whether the mother’s feeding choice (breast feeding or not) protects the infant against hospitalized pneumonia in the first year of life. The outcome variable of interest is the time to hospitalized pneumonia, which is restricted to one year and possibly censored before one year. The available baseline covariates are infant race (black, white or other), indicators of normal birthweight (at least 5.5 pounds) and having at least one sibling, and some maternal and family characteristics: age, years of schooling, alcohol use, cigarette use, region of the country (Northeast, Northcentral, South and West), poverty, and urban environment.

We are interested in the causal effect of breast feeding on hospitalized pneumonia as measured by θ with h defined by (2). This effect measure can be interpreted as a win-lose probability difference in a hypothetical comparison of the restricted times to hospitalized pneumonia of two randomly chosen infants, who are randomly assigned to breast feeding and no breast feeding as in a clinical trial. The propensity score model in our analysis is a logistic regression model based on all infant characteristics as well as maternal age, cigarette use, and years of schooling. The OR model is specified as a proportional hazards model based on infant characteristics, maternal cigarette use, and urban environment, with separate parameters for each treatment group. The model for censoring is also a separate proportional hazards model for each treatment group, with infant race and maternal age as covariates.

Table 2 shows the point estimates of θ obtained using the five methods compared in Section 3.1 together with nonparametric bootstrap standard errors based on 1,000 bootstrap samples. The IPW method clearly stands out with a large positive point estimate and a large standard error, both of which are likely due to the large variability of the IPW estimator. Considering the large standard error, the IPW estimate of θ may well be due to random variation and does not establish a positive effect of breast feeding. The naive method, on the other hand, has the smallest (i. e., most negative) point estimate and also a relatively large standard error. The results are somewhat similar for the other three methods (OR, DR1 and DR2). DR2 is the only method that reaches statistical significance (with or without adjusting for multiplicity); however, the estimated effect is rather small in magnitude under the DR2 method. Taken together, the point estimates and standard errors in Table 2 provide no clear evidence for a substantial effect of breast feeding on hospitalized pneumonia.