Home | About | Journals | Submit | Contact Us | Français |

**|**Biostatistics**|**PMC2800161

Formats

Article sections

- Abstract
- INTRODUCTION
- AVERAGE CAUSAL EFFECT ESTIMAND AND IDENTIFIABILITY ASSUMPTIONS
- ESTIMATION OF THE ACE
- SIMULATION STUDY
- EXAMPLE
- CONCLUDING REMARKS
- SUPPLEMENTARY MATERIAL
- FUNDING
- Supplementary Material
- References

Authors

Related links

Biostatistics. 2010 January; 11(1): 34–47.

Published online 2009 October 8. doi: 10.1093/biostatistics/kxp034

PMCID: PMC2800161

Peter B. Gilbert^{*}

Department of Biostatistics, University of Washington, Seattle, WA 98105, USA and Fred Hutchinson Cancer Research Center, 1100 Fairview Avenue North, Seattle, WA 98109, USA ; Email: pgilbert/at/scharp.org

Department of Biostatistics, University of Washington, Seattle, WA 98105, USA

Received 2008 December 30; Revised 2009 July 29; Accepted 2009 July 30.

Copyright © The Author 2009. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oxfordjournals.org.

This article has been cited by other articles in PMC.

In the past decade, several principal stratification–based statistical methods have been developed for testing and estimation of a treatment effect on an outcome measured after a postrandomization event. Two examples are the evaluation of the effect of a cancer treatment on quality of life in subjects who remain alive and the evaluation of the effect of an HIV vaccine on viral load in subjects who acquire HIV infection. However, in general the developed methods have not addressed the issue of missing outcome data, and hence their validity relies on a missing completely at random (MCAR) assumption. Because in many applications the MCAR assumption is untenable, while a missing at random (MAR) assumption is defensible, we extend the semiparametric likelihood sensitivity analysis approach of Gilbert *and others* (2003) and Jemiai and Rotnitzky (2005) to allow the outcome to be MAR. We combine these methods with the robust likelihood–based method of Little and An (2004) for handling MAR data to provide semiparametric estimation of the average causal effect of treatment on the outcome. The new method, which does not require a monotonicity assumption, is evaluated in a simulation study and is applied to data from the first HIV vaccine efficacy trial.

For randomized placebo-controlled efficacy trials of an HIV vaccine evaluated in HIV-uninfected volunteers, a primary objective is to evaluate the effect of vaccination on the incidence of HIV infection. Another objective, which was secondary in 2 trials of antibody-based vaccines (Flynn *and others*, 2005; Pitisittithum *and others*, 2006) and co-primary in 2 trials of T cell–based vaccines (Buchbinder *and others*, 2008), is to evaluate the effect of vaccination on HIV viral load measured after HIV infection. Two causal approaches of interest for assessing the latter objective are intent-to-treat (ITT) methods that assess the burden of illness (a composite end point that combines the infection end point and the postinfection end point) (Chang *and others*, 1994) and conditional methods that assess the postinfection end point in the principal stratum of subjects who would be HIV-infected under either treatment assignment (the so-called “always infected” stratum). As discussed by Gilbert *and others* (2003) [henceforth GBH] and others, the 2 approaches address different substantive questions. For example, in registrational/licensure trials, the ITT approach may be most appropriate for the primary analysis and the conditional approach would be used in secondary analyses of the mechanistic vaccine effect on the postinfection outcome, whereas in preregistrational test-of-concept trials, the conditional approach might be used for the primary analysis (Mehrotra *and others*, 2006).

Here, we develop methodology for the conditional approach to evaluate the causal vaccine effect on the postinfection end point. For concreteness, we refer to the first postrandomization event as infection and the outcome of interest measured after this event as viral load. However, the method has general application to evaluating causal treatment effects on outcomes measured after a postrandomization event, including studies of quality of life (Rubin, 2000), prostate cancer severity (Shepherd *and others*, 2008), kidney disease, and cancer screening (Joffe *and others*, 2007).

Because the set of trial participants who are in the always infected stratum is unknown, causal treatment effects for this group are not identified from randomized trial data. Two approaches for addressing the nonidentifiability have been to derive sharp bounds for causal effects (Jemiai and Rotnitzky, 2003; Hudgens *and others*, 2003; Zhang and Rubin, 2003) and to estimate causal effects under an additional set of identifiability assumptions that include models describing the nature and degree of possible selection bias, with a sensitivity analysis to explore how the inferences vary over a range of the selection models (GBH; Hayden *and others*, 2005; Hudgens and Halloran, 2006; Jemiai and Rotnitzky, 2005 (henceforth JR); Jemiai *and others*, 2007; Shepherd, Gilbert, Jemiai *and others*, 2006; Shepherd *and others*, 2007, 2008). All these methods assume that viral load is observed for all infected subjects. Therefore, the validity of their inferences depends on a missing completely at random assumption (MCAR). While the MCAR assumption is often untenable, trials may collect sufficient data on participant characteristics to make a missing at random (MAR) assumption plausible. This article extends the approach of GBH and JR to accommodate MAR missingness of the viral load end point.

To illustrate the new method, we focus on the first efficacy trial (Flynn *and others*, 2005) and define the end point *Y*_{p} to be the pre-antiretroviral therapy (pre-ART) log_{10} viral load at the month 12 visit post HIV infection diagnosis. Only subjects who have not started ART by the month 12 visit contribute a value *Y*_{p}. We exclude viral load values measured after ART initiation because ART strongly effects viral load levels (Gilbert *and others*, 2003). Inferences about *Y*_{p} apply to a population where ART is not prescribed during the first 12 months after infection diagnosis.

Of the 368 subjects who became HIV-infected during the trial, 121 had *Y*_{p} observed, 138 had missing data because they initiated ART prior to the month 12 study visit, and 109 had missing data because they dropped out prior to initiating ART and prior to the visit. Figure 1 shows the VaxGen data on *Y*_{p} = *Y*_{3} in relationship to 2 key covariates: *Y*_{1} is early pre-ART square root CD4 cell count, and *Y*_{2} is early pre-ART log_{10} viral load; these variables average all pre-ART values measured within 2 months of infection diagnosis. Viral load early and at 12 months were positively correlated (Figures 1(e) and (f)), with Spearman rank correlation 0.68 (0.55) for the vaccine (placebo) group. There is evidence of a negative correlation between early CD4 cell count and viral load at 12 months (Figures 1(c) and (d)) for the vaccine group (Spearman correlation = − 0.33) but not for the placebo group (Spearman correlation = 0.16).

VaxGen trial data: jittered pre-antiretroviral (pre-ART) viral loads at the month 12 postinfection diagnosis visit (*Y*_{3}) ((a) and (b)); pre-ART viral loads at the month 12 postinfection diagnosis visit versus square root CD4 cell counts early after infection **...**

In addition, lower levels of *Y*_{1} and higher levels of *Y*_{2} were highly predictive of ART initiation (*P* values <0.001 in a multivariate Cox model; Gilbert *and others*, 2005), showing that MCAR was badly violated. MAR may be plausible, however, because physicians base decisions to prescribe ART on the monitoring of CD4 cell counts, viral loads, HIV-related clinical events, and comorbidities (Hammer *and others*, 2008). In fact, MAR missingness due to ART initiation may approximately hold “systematically” in settings where ART is offered to all infected participants when their biomarkers are observed to cross prespecified thresholds or they present with prespecified symptomatic illnesses.

Two popular approaches to making valid inferences under MAR are inverse probability weighted (IPW)–based methods (e.g. Robins *and others*, 1995) and likelihood-based methods (e.g. Little and Rubin, 2002). For some HIV vaccine trials, the former methods are expected to be relatively inefficient and possibly unstable because the estimated weights for some subjects are expected to be near zero. Specifically, for a subject with *Y*_{p} observed, the weight in the denominator of the estimating equation equals the estimated probability that the subject did not drop out by the month 12 visit multiplied by the estimated probability that the subject did not start ART by the visit conditional on not dropping out. These estimates can be computed as fitted values in regression models based on the subject's covariates. The latter conditional probability may be near zero for subjects whose CD4 cell counts drop below the prespecified threshold at which ART initiation is recommended (currently the recommended threshold is between 200 and 500 cells/mm^{3} depending on the country). In fact, for ethical reasons any efficacy trial will offer ART to all infected participants who meet treatment criteria, such that the more successful the treatment coverage the closer some estimated weights in the denominator may be to zero. Therefore, for some HIV vaccine trials, IPW methods are expected to provide relatively imprecise inferences.

While likelihood methods are less subject to the instability problem, they are susceptible to misspecification of the model relating *Y*_{p} to covariates. To partially ameliorate this problem, we use the robust likelihood–based method of Little and An (2004) (henceforth LA) that is based on penalized splines of the propensity score. Our approach only handles monotone missing data (i.e. dropout); an alternative approach would handle nonmonotone missingness using parametric multiple imputation and Monte Carlo Markov chain, for example, by extending the method of Mogg and Mehrotra (2007) for analyzing viral load data to address postrandomization selection bias. The article is organized as follows. Section 2 describes the causal estimand of interest and identifiability assumptions. Section 3 shows how to combine the methodologies of GBH/JR and LA into a procedure for consistently estimating the causal estimand under an MAR assumption. Section 4 evaluates the new method in a simulation study. Section 5 applies the method to the first HIV vaccine efficacy trial, and Section 6 offers concluding remarks.

Let *Z* be treatment assignment, and let *X* be a *q*-vector of baseline covariates fully observed for everyone. Let *S* be the indicator of the postrandomization event. Subjects experiencing *S* = 1 are subsequently evaluated at *V* visits, where variables *Y*_{1},…, *Y*_{n1} are collected at visit 1, variables *Y*_{n1 + 1},…, *Y*_{n1 + n2} are collected at visit 2, and so on, with variables *Y*_{∑i = 1V − 1ni + 1},…, *Y*_{p} collected at visit *V*, where *p* = ∑_{i = 1}^{V}*n*_{i}. The entire collection of *p* variables measured after *S* = 1 is *Y* (*Y*_{1},…, *Y*_{p})^{′}, where *Y*_{p} is the outcome variable of interest. For *j* = 1,…, *p*, let *M*_{j} be the indicator of whether *Y*_{j} is missing and set *M* = (*M*_{1},…, *M*_{p})^{′}. The variables *Y* are only meaningful if *S* = 1; thus, *Y* and *M* are undefined if *S* = 0 and we denote this by *Y* = *M* = *. For HIV vaccine trials, *Z* is vaccination assignment (*Z* = 1, vaccine; *Z* = 0, placebo) and *S* is HIV infection diagnosis during the trial.

Each participant has potential infection outcome *S*(1) if assigned vaccine and *S*(0) if assigned placebo. For *Z* = 0,1, the potential outcomes *Y*(*Z*) (*Y*_{1}(*Z*),*Y*_{2}(*Z*),…, *Y*_{p}(*Z*))^{′} and *M*(*Z*) (*M*_{1}(*Z*),*M*_{2}(*Z*),…, *M*_{p}(*Z*))^{′} are defined if *S*(*Z*) = 1; otherwise *Y*(*Z*) * and *M*(*Z*) *. With *μ*_{z} *E*(*Y*_{p}(*z*)|*S*(0) = *S*(1) = 1) for *z* = 0,1, the “average causal effect (ACE)” estimand of interest is ACE *μ*_{1} − *μ*_{0}. Our goal is to estimate the ACE based on assumptions and the observed i.i.d. data (*Z*_{i},*X*_{i},*S*_{i},*M*_{i},*Y*_{i}), *i* = 1,…, *N*.

For subjects with *S* = 1, let *Y*_{obs} denote the components of *Y* that are observed and *Y*_{mis} denote the components of *Y* that are missing. Let *f* be the conditional cumulative distribution function (cdf) of *M* given *Y* and *S* = 1, *f*(*M*|*Y*,*S* = 1,*ν*), where *ν* denotes unknown parameters. MAR states that missingness depends only on the observed values *Y*_{obs}, that is,

For simplicity, we develop the methods for a setting where *Y*_{1},…, *Y*_{p − 1} are fully observed in infected subjects (with *S* = 1) and only *Y*_{p} has missing values. In Section 6, we describe how to extend the method to the case of monotone missing data.

Throughout, we make the following 3 assumptions.

- A1: Stable unit treatment values assumption (Rubin, 1978).
- A2: The treatment assignment
*Z*is independent of (*X*,*S*(0),*S*(1),*M*(0),*M*(1),*Y*(0),*Y*(1)). - A3: For infected subjects (with
*S*= 1), the missing data mechanism for*Y*is MAR.

All the papers for studying the causal vaccine effect in the always infected stratum cited in Section 1 assume A1 and A2; we refer the reader to these articles for discussion about their justification. As discussed in Section 1, the MAR assumption A3 may be quite plausible when missingness is mainly due to ART initiation and ART guidelines are used. If the missing data mechanism is believed to be non-MAR, then the methodology here may mislead and should be used with caution. An advantageous feature of the MAR assumption is that investigators can design clinical trials to collect the covariate data that make it plausible.

We next postulate additional assumptions that identify the ACE and are indexed by fixed sensitivity parameters. Following JR, we suppose 3 models, which we refer to collectively as A4:

(2.1)

(2.2)

(2.3)

where *g*_{0} and *g*_{1} are known invertible link functions whose inverses are continuous in *α*_{0} and *α*_{1}, *α*_{0} and *α*_{1} are unknown parameters to be estimated, and *β*_{0}, *β*_{1}, and are known sensitivity parameters that are varied over plausible ranges, where subject matter experts can help define the plausible ranges. With logit links *g*_{0} and *g*_{1}, *β*_{0} is interpreted as the difference in the log odds of infection in the vaccine group given infection in the placebo group with *y* versus *y* − 1 viral load, and *β*_{1} is interpreted similarly reversing the role of vaccine and placebo. The parameter is interpreted as the probability that a subject infected in the vaccine group would also be infected in the placebo group. Except for the method of JR, all of the previously developed methods cited above assume = 1 (i.e. monotonicity, that the vaccine does not increase the risk of infection for any subject), in which case the selection model (2.2) is superfluous and only model (2.1) is used, and the only sensitivity parameter is *β*_{0}. In this case, the methods of GBH and JR are equivalent. For greater applicability of the method, here, we allow for nonmonotonic settings by considering the trio of models (2.1–2.3).

For the case that *Y*_{p} is measured from all subjects with *S* = 1, JR proved that A1, A2, and A4 identify the ACE and developed unbiased estimating equations that can be solved to obtain a consistent asymptotically normal estimator of the ACE. Shepherd *and others* (2008) summarized this result together with a computational procedure for solving the equation and for obtaining the consistent sandwich variance estimator. With *θ* (*p*_{0},*p*_{1},*α*_{0},*α*_{1},*μ*_{0},*μ*_{1})^{′}, JR's estimating equation ∑_{i = 1}^{N}*U*_{i}(*θ*) = 0 is detailed using our notation in Appendix A ([A.1]–[A.6]) in the supplementary material available at *Biostatistics* online.

JR's estimating equation was derived based on the following relationships, which form the basis for integrating the LA and JR methods:

(3.1)

(3.2)

(3.3)

(3.4)

We show how to augment JR's estimating equations to handle missing data via LA's method. Following Section 5 of LA, we suppose *q* + *p* ≥ 3 so that there are at least 2 covariates. For infected subjects, separately for groups *Z* = 0 and *Z* = 1, define the logit of the propensity score for *Y*_{p} to be observed, given the covariates *X*,*Y*_{1},…, *Y*_{p − 1}:

(3.5)

LA observed that, conditional on the propensity score and assuming MAR, the missingness of *Y*_{p} does not depend on *Z*, *X*, *Y*_{1},…, *Y*_{p − 1} (Rosenbaum and Rubin, 1983). Consequently, each of the expectations in (3.1–3.4) can be written as

(3.6)

where *h*(*Y*_{p}) denotes a function of *Y*_{p}. Therefore, JR's estimating function *U*_{i}(*θ*) can be modified to accommodate MAR missingness, where the modified function *U*_{i}^{M}(*θ*) uses (3.6) instead of (3.1)–(3.4). The estimating function *U*_{i}^{M}(*θ*) is detailed in Appendix A ([A.7]–[A.12]) in the supplementary material available at *Biostatistics* online. The function *U*_{i}^{M}(*θ*) includes predicted values [*h*(*Y*_{p})] for *h*(*Y*_{p})= *g*_{z}^{−1}(*Y*_{p}; α_{0}, β_{0}) or *Y*_{p}*g*_{z}^{ − 1}(*Y*_{p};*α*_{0},*β*_{0}), for *z* = 0,1, and *U*_{i}^{M}(*θ*) becomes fully specified once LA's technique is applied to obtain the predicted values [*h*(*Y*_{p})].

We use LA's (Section 5) penalized spline propensity prediction method to predict *Y*_{p} for infected subjects with missing data. First, a logistic regression model is fit relating *M*_{p} to (*X*_{1},…, *X*_{q},*Y*_{1},…, *Y*_{p − 1})^{′}, which yields estimated propensities for all subjects with *S* = 1. Second, a spline regression model of *Y*_{p} on is fit using subjects for whom *Y*_{p} is observed. Beyond , additional covariates that predict *Y*_{p} may be entered into the regression model parametrically. From this fitted model, the value *E*[*Y*_{p}|*Y*_{1}^{*},*X*_{2},…, *X*_{q},*Y*_{1},…, *Y*_{p − 1},*S* = 1,*Z* = *z*] is predicted for each subject with *S**M*_{p} = 1 using his or her estimated propensity and other covariates. Third, as described below, fitted values and are computed using these regression fits and either numerical integration (for a general link function) or a closed-form equation (for an inverse probit link function). These steps may be done in the vaccine and placebo groups separately.

Applied to our setting, LA's specific propensity spline prediction model replaces one of the predictor variables, say *X*_{1}, by *Y*_{1}^{*} and supposes

(3.7)

for *z* = 0,1, where for *j* = 2,…, *q*, *s*_{zj}^{X}(*Y*_{1}^{*}) = *E*(*X*_{j}|*Y*_{1}^{*},*S* = 1,*Z* = *z*) is a spline for the regression of *X*_{j} on *Y*_{1}^{*}, where *X*_{j}^{*} = *X*_{j} − *s*_{zj}^{X}(*Y*_{1}^{*}); and for *j* = 1,…, *p* − 1, *s*_{zj}^{Y}(*Y*_{1}^{*}) = *E*(*Y*_{j}|*Y*_{1}^{*},*S* = 1,*Z* = *z*) is a spline for the regression of *Y*_{j} on *Y*_{1}^{*}, where *Y*_{j}^{*} = *Y*_{j} − *s*_{zj}^{Y}(*Y*_{1}^{*}). Furthermore, for *z* = 0,1, *r*_{z} is a parametric function with unknown parameter vector *γ*_{z} that satisfies *r*_{z}(*Y*_{1}^{*},0,…, 0,*γ*_{z}) = 0 for all *γ*_{z}. For subject *i* in group *z* with *S*_{i}*M*_{pi} = 1, the predicted value of *Y*_{p} is obtained as , where denotes the sample estimate of the spline *s*_{zj}^{X}, denotes the sample estimate of the spline *s*_{zj}^{Y}, *x*_{ij} is the realization of *X*_{ij}, and *y*_{ij} is the realization of *Y*_{ij}. The model (3.7) enters *X*_{j}^{*} and *Y*_{j}^{*} as covariates; in follow-up work, Zhang and Little (2005) showed that the LA method—and hence our sensitivity analysis method—remains valid if the *X*_{j}^{*},*Y*_{j}^{*} are replaced with the *X*_{j},*Y*_{j}.

Next, using (3.7), for a general link function, we compute the fitted values as , where Φ is the cdf of the standard normal distribution. The predicted value is also computed using numerical integration. For the case that *g*_{z} is the inverse probit function, numerical integration is not needed and the fitted values and can instead be computed as simple functions of [*Y*_{p}]. This procedure is described in Appendix A in the supplementary material available at *Biostatistics* online.

The estimate is computed by solving ∑_{i = 1}^{N}*U*_{i}^{M}(*θ*) = 0 (*U*_{i}^{M}(*θ*) is defined in Appendix A [A.7]–[A.12] in the supplementary material available at *Biostatistics* online) with the following steps.

- Step 1: Estimate
*p*_{0}Pr(*S*(0) = 1) by solving ∑_{i = 1}^{N}*U*_{1i}^{M}(*p*_{0}) = 0, and estimate*p*_{1}Pr(*S*(1) = 1) by solving ∑_{i = 1}^{N}*U*_{2i}^{M}(*p*_{1}) = 0. - Step 2: Plug and the fitted values into
*U*_{3i}^{M}(*θ*) and solve for*α*_{0}in ∑_{i = 1}^{N}*U*_{3i}^{M}(*α*_{0}) = 0 with a 1D line search. Similarly, plug and the fitted values into*U*_{4i}^{M}(*θ*) and solve for*α*_{1}. - Step 3: Plug the estimates of
*p*_{0}and*α*_{0}and the fitted values into*U*_{5i}^{M}(*θ*) and solve ∑_{i = 1}^{N}*U*_{5i}^{M}(*μ*_{0}) = 0 for*μ*_{0}. Similarly, solve ∑_{i = 1}^{N}*U*_{6i}^{M}(*μ*_{1}) = 0 for*μ*_{1}.

Under the monotonicity assumption (i.e. = 1), the computational algorithm is the same except that the second part of Step 2 is omitted (*α*_{1} is no longer relevant), and *U*_{6i}^{M}(*θ*) simplifies (given in expression (A.13) of Appendix A in the supplementary material available at *Biostatistics* online). If *g*_{z} is inverse probit, then Steps 2 and 3 use closed-form formulas (A.14) and (A.15) in Appendix A in the supplementary material available at *Biostatistics* online; otherwise they use numerical integration as described in Section 3.3.

Deriving asymptotic-based standard error estimators for and confidence intervals (CIs) for ACE is difficult, given the smoothing and tuning parameter selection involved with the penalized splines. We develop a bootstrap approach. As LA did not develop CIs nor standard error estimators, the performance of this approach is of interest even for the setting without a postrandomization event.

Within each treatment group *Z* = *z* separately, *B* bootstrap data sets are constructed by sampling with replacement *N*_{z} realizations of (*X*_{i},*S*_{i},*M*_{i},*Y*_{i})|*Z*_{i} = *z*, *z* = 0,1. The estimation procedure described in Section 3.4 is carried out for each bootstrap data set. Then, standard errors for , , and are estimated by the sample standard deviations of the bootstrap estimates, and (1 − *α*) × 100% CIs are obtained as the *α*/2 and 1 − *α*/2 percentiles of the bootstrap estimates.

LA's Theorem 1 states a double robustness–type property of the propensity spline prediction procedure, see Kang and Schafer (2007) and comments within for further discussion. In our setting, this robustness property can be stated as follows.

Theorem 3.1

Assume A1–A4, and let and be the estimators of *μ*_{0} and *μ*_{1}, respectively, obtained by solving ∑_{i = 1}^{N}*U*_{i}^{M}(*θ*) = 0. Then, for *z* = 0,1, is a consistent estimator of *μ*_{z} if either (i) the mean of *Y*_{p} conditional on (*Y*_{1}^{*},*X*_{2},…, *X*_{q},*Y*_{1},…, *Y*_{p − 1},*S* = 1,*Z*) in (3.7) is correctly specified or (ii) the propensity *Y*_{1}^{*} is correctly specified, and (iii) *E*[*X*_{j}^{*}|*Y*_{1}^{*},*S* = 1,*Z* = *z*] = *s*_{zj}^{X}(*Y*_{1}^{*}) for *j* = 2,…, *q* and *E*[*Y*_{j}^{*}|*Y*_{1}^{*},*S* = 1,*Z* = *z*] = *s*_{zj}^{Y}(*Y*_{1}^{*}) for *j* = 1,…, *p*.

LA comment that the robustness property of (iii) is that *r*_{z} does not need to be correctly specified and suggest that (b2) is a weak assumption because of the flexibility of the spline regression models.

We briefly describe the simulation study of an HIV vaccine trial, with expanded details in Appendix B in the supplementary material available at *Biostatistics* online.

Four simulation experiments were carried out such that A1–A4 hold with = 1 (i.e. monotonicity holds) and *g*_{0}^{ − 1}(*y*;*α*_{0},*β*_{0}) = Φ(*α*_{0} + *β*_{0}*y*). We created moderate or extreme selection bias specified by *β*_{0} = 1 or *β*_{0} = 3. Based on the VaxGen trial, we consider the 3 variables *Y*_{1}, *Y*_{2}, and *Y*_{3} defined in Section 1 (the outcome *Y*_{3} is month 12 pre-ART log_{10} viral load; *Y*_{1} is early pre-ART square root CD4 cell count; *Y*_{2} is early pre-ART log_{10} viral load).

Each simulation generated data with an even chance of assignment to vaccine or placebo, an overall placebo infection rate of 25%, and *Y*(0) = (*Y*_{1}(0),*Y*_{2}(0),*Y*_{3}(0))^{′} multivariate normal, with means and variances equal to the sample estimates obtained in the VaxGen trial (Gilbert *and others,* 2005). The correlations were set to cor(*Y*_{1}(0),*Y*_{2}(0)) = *c**o**r*(*Y*_{1}(0),*Y*_{3}(0)) = − 0.5 and cor(*Y*_{2}(0),*Y*_{3}(0)) = 0.8. For infected vaccine group subjects, *Y*(1) was set to *Y*(0) + (0,0, − 0.5)^{′}, which specifies the true ACE = − 0.5.

We assumed complete data for *Y*_{1} and *Y*_{2} in infected subjects and created missing data for *Y*_{3} that was either MCAR or MAR. In either case, *Y*_{3} is missing in about 50% of infected subjects. The MAR scenario generated missingness indicators from a logistic regression model of *M*_{3} on *Y*_{1} and *Y*_{2} with coefficients chosen from a fit to the VaxGen data.

Four methods were used to estimate the ACE with a 95% bootstrap percentile CI. The first was ordinary least squares complete case (OLSCC), which simply compares sample averages of *Y*_{3} between groups in all infected subjects that have *Y*_{3} measured. The next 3 methods are causal: before deletion (BD), which applies GBH/JR to *Y*_{3} using the complete data (an unknowable gold standard for comparison); CC, a complete case analysis that applies GBH/JR; SPPL, the new method with penalized spline propensity prediction based on a regression of *Y*_{3} on the spline of *Y*_{1}^{*}, where *Y*_{1}^{*} is the linear predictor of the estimated propensity to observe *Y*_{3} from a linear logistic regression of *M*_{3} on *Y*_{1} and *Y*_{2}. This SPPL method uses 15 equally spaced fixed knots and a truncated linear basis and uses model (3.7) with a linear additive model for *Y*_{2}^{*}:*r*_{z}(*Y*_{1}^{*},*Y*_{2}^{*},*γ*_{z}) = *γ*_{z1}*Y*_{1}^{*} + *γ*_{z2}*Y*_{2}^{*}, for *z* = 0,1. Note that the OLSCC method is equivalent to the GBH/JR method for the special case that = 1 and *β*_{0} = *β*_{1} = 0.

Table 1 shows bias and root mean squared error (RMSE) of the estimated ACE and the coverage probability of the nonparametric bootstrap CIs. In all cases, the OLSCC method is badly biased because it ignores the postrandomization selection bias. As expected, the CC method is unbiased under MCAR missingness but biased under MAR missingness, whereas the new SPPL method is unbiased under both MCAR and MAR data. Corresponding to this, the CIs about the ACE for the CC method have too-low coverage probability under MAR data, whereas the CIs for the SPPL method are always near the correct level. The RMSEs for the SPPL method are 21.1% lower on average than those of the CC method, showing that the new method provides substantially better estimators than the existing methods.

Figure 2 shows distributions of the CI lengths, demonstrating that the SPPL method provides much more precise estimates than the CC method. Furthermore, under MAR missingness the bootstrap variance estimates of had a similar distribution as the Monte Carlo sample variances of the ′s for the SPPL method but not for the CC method, indicating that the bootstrap variance estimator for the SPPL method is unbiased (not shown). Additional simulations were conducted under the above design except that *Y*(1) was set to *Y*(0) + (0, − 0.5, − 0.5)^{′}. The results were similar, demonstrating that the new method performs well when there is an average causal effect on a postinfection covariate used to predict the outcome of interest.

We now apply the newly developed SPPL method to the VaxGen data. The same method as evaluated in the simulations was applied, except with 3000 bootstrap replications. Our outcome of interest *Y*_{3} and covariates *Y*_{1} and *Y*_{2} are the same as considered above. Of the 368 subjects who became HIV-infected, there is nearly complete data for *Y*_{1} (*n* = 334) and *Y*_{2} (*n* = 341), and, as discussed above, one-third of the infected subjects have *Y*_{3} observed (*n* = 121) (Figure 1). The *R*^{2} value of a linear regression of *Y*_{3} on *Y*_{1} and *Y*_{2} is 0.40, suggesting that accounting for the predictive information in *Y*_{1} and *Y*_{2} with the new method is expected to provide greater efficiency than the CC method that does not leverage this positive correlation.

Because the estimated incidence of HIV infection was slightly lower in the vaccine group than the placebo group (Flynn *and others*, 2005), the data support (but of course do not prove) the monotonicity assumption that = 1. Accordingly, we first conduct a sensitivity analysis assuming = 1. For the sensitivity parameter *β*_{0} ranging from − 2 to 2, which in our opinion covers the plausible range for potential selection bias, Figure 3 shows point and 95% bootstrap CI estimates of the ACE by the CC method and by the new method. Both analyses support no causal vaccine effect on viral load, and the new method provides slightly narrower CIs.

Estimated ACE and 95% CIs with = 1. Solid lines are for the SPPL method and dashed lines are for the CC method.

Acknowledging that monotonocity may not hold, we next performed a sensitivity analysis setting = 0.8. This value was chosen based on the fact that the upper 95% confidence limit for the hazard ratio (vaccine/placebo) was 1.17. A value = 0.8 assumes that a vaccinated subject who becomes infected would have a 20% chance of avoiding infection had he/she been assigned placebo, which approximately specifies a 17% plausible elevation of infection risk in the vaccine group. The shaded regions in Figure 4 show values of (*β*_{0}, *β*_{1}) under which a 95% CI for the ACE excludes 0. Figure 4 shows that only under conditions of large selection bias in certain directions are the data consistent with a beneficial or harmful effect. We conclude that the data are consistent with the null hypothesis of ACE = 0.

(a) Estimated ACE, (b) lower 95% confidence limits, and (c) upper 95% confidence limits with = 0.8. Solid lines are for the SPPL method and dashed lines are for the CC method.

We conducted the example using an inverse probit link function *g*_{z}. With this link, interpreting the sensitivity parameters *β*_{z} requires thinking in the “*Z*-metric” (normal quantitle metric), where a one unit increase in *Y*_{p} leads to an increase in the probit score Φ^{ − 1}(Pr(*S*(1) = 1|*S*(0) = 1,*Y*_{p}(0) = *y*) by *β*_{z} standard deviations. For those uncomfortable thinking in the *Z*-metric, it may be preferable to use an alternative link function that they find more interpretable, for example, logit. Furthermore, it is good practice in sensitivity analysis to repeat the analysis for multiple link functions.

Recently developed methods for estimating causal treatment effects on an outcome measured after a postrandomization event, including the methods of GBH and JR, are valid under an MCAR assumption. We have extended the GBH and JR methods to allow for MAR missingness, by augmenting it with LA's robust likelihood approach that uses penalized spline propensity prediction. The simulation study showed that the new method corrects for the bias of the existing methods and improves precision by incorporating covariates that predict the outcome.

We have focused on accommodating missing data on the outcome of interest *Y*_{p}, and as such have supposed complete data for the postrandomization event *S*. The method could be extended to accommodate MAR missingness of *S* via IPW approach of Shepherd *and others*, 2008 to address this issue.

If the missing outcome data follow a monotone pattern, such that *Y*_{j + 1},…, *Y*_{p} are missing for all subjects for whom *Y*_{j} is missing, then the method can be extended using the approach described in Section 7 of LA. Specifically, the propensity spline model (3.7) can be applied sequentially to each block of missing values. With this approach, the fitted values and appearing in the estimating function *U*_{i}^{M}(*θ*) are computed by replacing missing covariate values by predicted values in sequential regression models.

Extensions of interest include accommodating nonmonotone missing data and expanding the sensitivity analysis to evaluate the impact on inferences of deviations from the MAR assumption.

Supplementary material is available at http://biostatistics.oxfordjournals.org.

National Institutes of Health (RO1 AI054165-06) to P.B.

*Conflict of Interest*: None declared.

- Buchbinder SP, Mehrotra DV, Duerr A, Fitzgerald DW, Mogg R, Li D, Gilbert PB, Lama JR, Marmor M, Del Rio C.
*and others*. The Step study: the first test-of-concept efficacy trial of a cell-mediated immunity HIV vaccine. Lancet. 2008;372:1881–1893. [PMC free article] [PubMed] - Chang MN, Guess HA, Heyse JF. Reduction in the burden of illness: a new efficacy measure for prevention trials. Statistics in Medicine. 1994;13:1807–1814. [PubMed]
- Flynn NM, Forthal DN, Harro CD, Judson FN, Mayer KH, Para MF. The rgp120 HIV Vaccine Study Group. Placebo-controlled phase 3 trial of recombinant glycoprotein 120 vaccine to prevent HIV-1 infection. The Journal of Infectious Diseases. 2005;191:654–665. [PubMed]
- Gilbert PB, Ackers ML, Berman PW, Francis DP, Popovic V, Hu DJ, Heyward WL, Sinangil F, Shepherd B, Gurwith M. HIV-1 virologic and immunologic progression and antiretroviral therapy initiation among HIV-1-infected participants in an efficacy trial of recombinant glycoprotein 120 vaccine. The Journal of Infectious Diseases. 2005;192:974–983. [PubMed]
- Gilbert PB, Bosch R, Hudgens MG. Sensitivity analysis for the assessment of causal vaccine effects on viral load in HIV vaccine trials. Biometrics. 2003;59:531–541. [PubMed]
- Hammer SM, Eron JJ, Jr, Reiss P, Schooley RT, Thompson MA, Walmsley S, Cahn P, Fischl MA, Gatell JM, Hirsch MS.
*and others*. Antiretroviral treatment of adult HIV infection: 2008 recommendations of the International AIDS Society—USA panel. Journal of the American Medicial Association. 2008;300:555–570. [PubMed] - Hayden D, Pauler DK, Schoenfeld D. An estimator for treatment comparisons among survivors in randomized trials. Biometrics. 2005;61:305–310. [PubMed]
- Hudgens MG, Halloran ME. Causal vaccine effects on binary postinfection outcomes. Journal of the American Statistical Association. 2006;101:51–64. [PMC free article] [PubMed]
- Hudgens MG, Hoering A, Self SG. On the analysis of viral load endpoints in HIV vaccine trials. Statistics in Medicine. 2003;22:2281–2298. [PubMed]
- Jemiai Y, Rotnitzky A. Semiparametric methods for inferring treatment effects on outcomes defined only if a post-randomization event occurs, [Unpublished PhD. Dissertation]. Harvard School of Public Health, Department of Biostatistics. 2005
- Jemiai Y, Rotnitzky A, Shepherd B, Gilbert PB. Semiparametric estimation of treatment effects given base-line covariates on an outcome measured after a post-randomization event occurs. Journal of the Royal Statistical Society, Series B. 2007;69:879–902. [PMC free article] [PubMed]
- Joffe MM, Small D, Hsu C-Y. Defining and estimating intervention effects for groups that will develop an auxiliary outcome. Statistical Science. 2007;22:74–97.
- Kang JDY, Schafer JL. Demystifying double robustness: a comparison of alternative strategies for estimating a population mean from incomplete data. Statistical Science. 2007;22:523–539. [PMC free article] [PubMed]
- Little R, An H. Robust likelihood-based analysis of multivariate data with missing values. Statistica Sinica. 2004;14:949–968.
- Little RJA, Rubin DB. Statistical Analysis with Missing Data. New York: Wiley; 2002.
- Little R, Zhang G. Extensions of the penalized spline of propensity prediction method of imputation. Biometrics. 2008 doi: 10.1111/j.1541-0420.2008.01155.x. [PubMed]
- Mehrotra DV, Li X, Gilbert PB. A comparison of eight methods for the dual-endpoint evaluation of efficacy in a proof-of-concept HIV vaccine trial. Biometrics. 2006;62:893–900. [PubMed]
- Mogg R, Mehrotra DV. Analysis of antiretroviral immunotherapy trials with potentially non-normal and incomplete longitudinal data. Statistics in Medicine. 2007;26:484–497. [PubMed]
- Pitisuttithum P, Gilbert PB, Gurwith M, Heyward W, Martin M, Van Griensven F, Hu D, Tappero JW, Choopanya K. The Bangkok Vaccine Evaluation Group. Randomized, double-blind, placebo-controlled efficacy trial of a bivalent recombinant glycoprotein 120 HIV-1 vaccine among injection drug users in Bangkok, Thailand. The Journal of Infectious Diseases. 2006;194:1661–1671. [PubMed]
- Robins JM, Rotnitzky A, Zhao LP. Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association. 1995;90:106–121.
- Rosenbaum PR, Rubin DB. Assessing sensitivity to an unobserved binary covariate in an observational study with binary outcome. The Journal of the Royal Statistical Society, Series B. 1983;45:212–218.
- Rubin DB. Bayesian inference for causal effects. Annals of Statistics. 1978;6:34–58.
- Rubin DB. Comment on “Causal inference without counterfactuals” by A.P. Dawid. Journal of the American Statistical Association. 2000;95:435–437.
- Shepherd B, Gilbert PB, Jemiai Y, Rotnitzky A. Sensitivity analyses comparing outcomes only existing in a subset selected post-randomization, conditional on covariates, with application to HIV vaccine trials. Biometrics. 2006;62:332–342. [PubMed]
- Shepherd B, Gilbert PB, Lumley T. Sensitivity analyses comparing time-to-event outcomes only existing in a subset selected post-randomization, conditional on covariates, with application to HIV vaccine trials. Journal of the American Statistical Association. 2007;102:573–582. [PMC free article] [PubMed]
- Shepherd B, Gilbert PB, Mehrotra DV. Eliciting a counterfactual sensitivity parameter. The American Statistician. 2006;61:56–63.
- Shepherd BE, Redman MW, Ankerst DP. Does finasteride affect the severity of prostate cancer? A causal sensitivity analysis. Journal of the American Statistical Association. 2008;103:1392–1404. [PMC free article] [PubMed]
- Zhang JL, Rubin DB. Estimation of causal effects via principal stratification when some outcomes are truncated by “death. Journal of Educational and Behavioral Statistics. 2003;28:353–368.

Articles from Biostatistics (Oxford, England) are provided here courtesy of **Oxford University Press**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |