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

**|**HHS Author Manuscripts**|**PMC2700847

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Data Structure and Notation
- 3 Causal Estimand
- 4 Assumptions
- 5 Identification, Bias, and Sensitivity Analysis
- 6 Estimation
- 7 NSCOT Study
- 8 Simulation Study
- 9 Discussion
- 11 References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2010 June 1.

Published in final edited form as:

PMCID: PMC2700847

NIHMSID: NIHMS52102

Ph: 215-214-3917, e-mail: ude.cccf@notselgE.nairB

The publisher's final edited version of this article is available at Biometrics

See other articles in PMC that cite the published article.

We focus on estimation of the causal effect of treatment on the functional status of individuals at a fixed point in time *t** after they have experienced a catastrophic event, from observational data with the following features: (1) treatment is imposed shortly after the event and is non-randomized, (2) individuals who survive to *t** are scheduled to be interviewed, (3) there is interview non-response, (4) individuals who die prior to *t** are missing information on pre-event confounders, (5) medical records are abstracted on all individuals to obtain information on post-event, pre-treatment confounding factors. To address the issue of survivor bias, we seek to estimate the survivor average causal effect (SACE), the effect of treatment on functional status among the cohort of individuals who would survive to *t** regardless of whether or not assigned to treatment. To estimate this effect from observational data, we need to impose untestable assumptions, which depend on the collection of all confounding factors. Since pre-event information is missing on those who die prior to *t**, it is unlikely that these data are missing at random (MAR). We introduce a sensitivity analysis methodology to evaluate the robustness of *SACE* inferences to deviations from the MAR assumption. We apply our methodology to the evaluation of the effect of trauma center care on vitality outcomes using data from the National Study on Costs and Outcomes of Trauma Care.

Consider an observational study of individuals who have experienced a catastrophic event and are non-randomly assigned to a treatment or control condition shortly after the event. The study design calls for the collection of information on functional outcomes based on interviews among survivors at a fixed point in time *t** after the event, collection of information on pre-injury covariate data on all survivors, and abstraction of medical records on all individuals to obtain post-event, pre-injury covariate data. The scientific objective is to draw inferences about the effect of treatment on functional outcomes at time *t**.

Inference about functional outcomes is complicated due mortality. A traditional analysis which compares treatment-specific outcomes among observed survivors can be misleading, even after adjusting for measured covariates. This follows since the analysis is conditioned on an intermediate outcome (i.e., survival) which may be affected by treatment (Rosenbaum, 1984). To see this, suppose some individuals would survive regardless of their treatment assignment and possess moderate functional outcomes under the control condition and high functional outcomes under treatment; some who would die regardless of their treatment assignment; and the remaining would die under the control condition, but survive under the treatment condition with low functional outcomes. Thus, treatment has a positive effect on survival and a positive impact on functional outcomes among those who would survive regardless of treatment assignment. The untreated survivors will have moderate functional outcomes. In contrast, the treated survivors are a mixture of those who would live regardless of treatment assignment and those who live under treatment only and accordingly, the observed functional outcomes would be mixture of high and low scores. Depending on the mixture, the traditional analysis could suggest that treatment has a null or negative impact on functioning, even though treatment is obviously beneficial both in terms of survival and functioning.

To address this issue, Rubin (2000) and Robins and Greenland (2000) recommended that inference should focus on the effect of the treatment on functional outcomes among the cohort of individuals who would survive to time *t** regardless of whether or not they are assigned to treatment. The effect has been termed the “Survivor Average Causal Effect” (*SACE*; Hayden *et al.*, 2005).

Egleston *et al.*, (2007) proposed a sensitivity analysis methodology for identifying and estimating *SACE* with complete information on confounders and missing outcome data among survivors. Their assumptions do not identify *SACE* when there are confounding factors that are missing exclusively on those who die. This paper demonstrates how the method of sensitivity analysis can be used to identify and estimate *SACE* when confounders and outcomes are missing. Frangakis *et al.* (2007) also tackled the issue of missing factors due to death; however their approach does not identify *SACE*.

This paper is organized as follows. Section 2 introduces the data structure and notation. Section 3 defines the *SACE* estimand. Section 4 provides a review of the identification and estimation of *SACE* when there is complete information on confounders. Sections 5 and 6 discuss identification and estimation of *SACE* in the presence of missing confounder data. Section 7 provides an analysis of data from the National Study on the Costs and Outcomes of Trauma Care (NSCOT; MacKenzie *et al.*, 2006), in which moderately to severly injured individuals were assigned to received trauma or non-trauma center care. Section 8 discusses the results of a simulation study. The final section is devoted to a discussion.

Let *Z* denote treatment indicator, where *Z* = 1 (*Z* = 0) denotes that the individual has been assigned to treatment (control). Let *D*(*z*) be the indicator of death prior to time *t** had the individual been assigned to treatment *Z* = *z*. If *D*(*z*) = 0, then define *Y* (*z*) to be the functional outcome (continuous) at time *t** had the individual been assigned to treatment *Z* = *z*. The outcomes {*D*(0), (*Y*(0) : *D*(0) = 0), *D*(1), (*Y*(1) : *D*(1) = 0)} are referred to as potential outcomes, as they are only partially observable. The notation “:” represents belonging to the relevant subgroup. Let *D* = *D*(*Z*) be the observed mortality outcome. If *D* = 0, then define *Y* = *Y*(*Z*) as the functional outcome under the treatment actually received. Also, let *X* = {*X*_{1}, *X*_{2}} denote a set of confounding covariates, where *X*_{1} is observed and *X*_{2} is observed if *D* = 0. If *D* = 0, let *R* be an indicator of whether the functional outcome is observed. The observed data for an individual is *O* = {*Z, R, D, X*_{1}, (*X*_{2} : *D* = 0), (*Y* : *D* = 0, *R* = 1)}

For continuous outcomes, *SACE* has been defined as the mean of the difference in the potential functional outcomes under treatment and control, among the cohort of individuals who would survive regardless of whether or not assigned to treatment. That is,

$$\mathit{SACE}=E[Y\left(1\right)-Y\left(0\right)\mid D\left(0\right)=0,D\left(1\right)=0]$$

Egleston *et al.* (2007) proposed a number of non-verifiable assumptions, similar in spirit to Hayden *et al.*’s (2005) assumptions, to identify *SACE* when there were completely observed confounder data. These assumptions are as follows:

The stable unit treatment value assumption (SUTVA)

SUTVA (Rubin, 1980) states that an individual’s potential outcomes under treatment *z* (i.e., {*D*(*z*), (*Y*(*z*) : *D*(*z*) = 0)}) remain the same regardless of the mechanism used to assign the individual to treatment *z and* regardless of the treatment assignments of other individuals. This assumption would be violated if there is interference between individuals *or* if there exist unrepresented versions of treatments.

*D*(1) ≤ D(0) (Monotonicity)

Monotonicity states that treatment does not harm survival. If an individual were to live under the control condition then she would also live under the treatment condition.

*Z*{*D*(0), *Y*(0) : *D*(0) = 0, *D*(1), *Y*(1) : *D*(1) = 0} | *X* (No unmeasured confounders)

This assumption is the observational study equivalent of randomization in a randomized trial. It states that within levels of the confounding covariates *X*, treatment assignment is randomized. This implies that, within levels of *X*, the distribution of potential outcomes is the same for those assigned to treatment as those assigned to control.

*E*[*Y*(1)|*D*(0) = 0, *D*(1) = 0, *X*] = *E*[*Y*(1)|*D*(0) = 1, *D*(1) = 0, *X*]

This is consistent with one component of the “explainable non-random survival” assumption of Hayden *et al.* (2005). It states that, given survival under treatment and the full set of confounding factors, the mean functional outcome under treatment is unrelated to mortality experience under control. Hayden *et al.* (2005) and Egleston *et al.* (2007) propose sensitivity analyses related to this assumption.

*R**Y*|*D* = 0, *X*, *Z* (Missing at random)

This assumption states that missingness of the functional outcomes among survivors is not related to the outcome given the covariates and exposure.

We will discuss each of these assumptions in the context of NSCOT study in Section 7.

Egleston *et al.* (2007) demonstrated that, under the above assumptions,

$$E[Y\left(z\right)\mid D\left(0\right)=0,\phantom{\rule{thickmathspace}{0ex}}D\left(1\right)=0]=\frac{E\left[\mu (z,X)\rho (0,X)\right]}{E\left[\rho (0,X)\right]}.$$

(5.1)

where *μ*(*z, X*) = *E*[*Y*|*D* = 0, *R* = 1, *X*, *Z* = *z*] and *ρ*(*z, X*) = *P*[*D* = 0|*X*, *Z* = *z*]. When *X* is fully observed on all subjects, *μ*(*X, z*), *ρ*(*z, X*), and the marginal distribution of *X* are identifiable. When *X*_{2} is missing on subjects who die, *μ*(*z, X*) remains identifiable since it represents the expected outcomes of those who survive to one year and respond to the interview, and hence have observable *Y* and *X*. However, *ρ*(*z, X*) and the marginal distribution of *X* are not identified since they do not condition on survival and hence *X*_{2} is not observable on those who die.

The key to identification of *SACE* when *X*_{2} is missing is the following key identity:

$$E\left[k(Z,X)\right]=E\left[\frac{(1-D)}{\rho (Z,X)}k(Z,X)\right]=P[D=0]\int \frac{k(\stackrel{~}{z},\stackrel{~}{x})}{\rho (\stackrel{~}{z},\stackrel{~}{x})}dF(\stackrel{~}{z},\stackrel{~}{x}\mid D=0)$$

(5.2)

where *k*(*Z*, *X*) is any function of *Z* and *X*. In light of Equation (5.1), this identity implies that, if we could establish conditions for identifying *ρ*(*z*, *x*), then *SACE* would be identifiable since we would only need to take expectations with respect to the distribution of *Z* and *X* among survivors, which is identifiable. Specifically,

$$E[Y\left(z\right)\mid D\left(0\right)=0,D\left(1\right)=0]=\frac{E\left[\frac{(1-D)}{\rho (Z,X)}\mu (z,X),\rho (0,X)\right]}{E\left[\frac{(1-D)}{\rho (Z,X)}\rho (0,X)\right]}=\frac{\int \frac{\mu (z,\stackrel{~}{x})\rho (0,x)}{\rho (\stackrel{~}{z},\stackrel{~}{x})}dF(\stackrel{~}{z},\stackrel{~}{x}\mid D=0)}{\int \frac{\rho (0,x)}{\rho (\stackrel{~}{z},\stackrel{~}{x})}dF(\stackrel{~}{z},\stackrel{~}{x}\mid D=0)}$$

(5.3)

To identify *ρ*(*z*, *x*), we introduce a class of assumptions indexed by sensitivity analysis parameters.

logit {1-*ρ*(*z*, *x*)} = *h*(*z*, *x*_{1}) + *α*’*x*_{2} + *τ*’ *z*·*x*_{2}, where *h*(*z*, *x*_{1}) is an unknown function of *z* and *x*_{1}, *z* · *x*_{2} is a vector where each component of *x*_{2} is multiplied by *z*, and *α* and *τ* are fixed *non-identifiable* sensitivity analysis vectors of the same dimension as *x*_{2}.

Here, we assume no interaction between *x*_{1} and *x*_{2}. However, additional sensitivity analysis parameters could be added to the model to account for interactions between *x*_{1} and *x*_{2}. We also assume that *ρ*(*z*, *x*) > 0 for all *z* and *x*. Under assumption 6 with *α* and *τ* fixed, it can be shown that *ρ*(*z*, *x*) is identifiable from the observed data. Identification of *ρ*(*z*, *x*) is tantamount to showing that *h*(*z*, *x*_{1}) is identified from the observed data. It can be shown that

$$h(z,{x}_{1})=\mathrm{log}\left\{\frac{E[D\mid Z=z,{X}_{1}={x}_{1}]}{E[(1-D)\mathrm{exp}({\alpha}^{\prime}{X}_{2}+{\tau}^{\prime}Z\cdot {X}_{2})\mid Z=z,{X}_{1}={x}_{1}]}\right\},$$

where each of the conditional expectations on the right hand side above are identifiable from the distribution of the observed data. This result tells us that *SACE* is identifiable.

When ** α** =

What is the implication of falsely assuming that ** α** =

$${\mathit{SACE}}_{x1}=\sum _{x2=0}^{1}{E}_{x1}[Y\left(1\right)-Y\left(0\right)\mid D\left(0\right)=0,D\left(1\right)=0,{X}_{2}={x}_{2}]{P}_{x1}[{X}_{2}={x}_{2}]\frac{{P}_{x1}[D\left(0\right)=0\mid {X}_{2}={x}_{2}]}{{P}_{x1}[D\left(0\right)=0]}$$

(5.4)

When ** α** =

$${\mathit{SACE}}_{x1}=\sum _{x2=0}^{1}{E}_{x1}[Y\left(1\right)-Y\left(0\right)\mid D\left(0\right)=0,{X}_{2}={x}_{2}]{P}_{x1}[{X}_{2}={x}_{2}]$$

(5.5)

The right hand sides of (5.4) and (5.5) will be unequal since *P*_{x1}[*D*(0) = 0|*X*_{2} = 1] < *P*_{x1}[*D*(0) = 0|*X*_{2} = 0]. If 0 < *P*[*X*_{2} = 1] < 1 and *E*_{x1}[*Y*(1) - *Y*(0)|*D*(0) = 0, *X*_{2} = 1] > (<)*E*_{x1}[*Y*(1) - *Y*(0)|*D*(0) = 0, *X*_{2} = 0], then falsely assuming ** α** =

Since ** α** and

In order to estimate *SACE* from a dataset of size *N*, we need to estimate *μ*(*z, x*) and *h*(*z, x*_{1}). When *X*_{1} and *X* are high-dimensional (i.e., multiple discrete or continuous components), we cannot estimate these quantities non-parametrically at rates fast enough to guarantee that the resulting estimator of *SACE* will be consistent at -rates. This is called the curse of dimensionality (Robins and Ritov, 1997). Thus, we need to impose modeling restrictions on *μ*(*z*, *x*) and *h*(*z*, *x*_{1}). In particular, we assume that *h*(*z*, *x*_{1}) = *h*(*z, x*_{1}; *γ**) and *g*(*μ*(*z*, *x*)) = *b*(*z*, *x*; *β**). where *h*(*z*, *x*_{1}; *γ*) is a specified function of *z, x*_{1} and *γ*, *g*(·) is a specified link function, *b*(*z*, *x*; *β*) is a specified function of *z*, *x*, and *β*, and *γ** and *β** denote the true values of *γ* and *β*, respectively.

Throughout, we suppress dependence of our estimation procedure on ** α** and

To estimate *γ** and *β**, we develop unbiased estimating functions. Our estimating function for *γ** is motivated by the fact that

$${U}_{\gamma}(O;\gamma )={a}^{\prime}(Z,X;\gamma )(D-(1-D)\mathrm{exp}\left\{a(Z,X;\gamma )\right\})$$

has mean zero when evaluated at *γ**, where *a*’(*Z*, *X*; *γ*) is the derivative of *h*(*Z*, *X*_{1}; *γ*) with respect to *γ*. Our estimating function for *β** is motivated by the fact that

$${U}_{\beta}(O;\beta )={b}^{\prime}(Z,X;\beta )R(1-D)(Y-\mu (Z,X;\beta ))$$

has mean zero when evaluated at *β**, where *b*’ (*Z*, *X*; *β*) is the derivative of *b*(*Z*, *X*_{1}; *β*) with respect to *β*. Further, we can set up an unbiased estimating function for ${\nu}_{z}^{\ast}=E[Y\left(z\right)\mid D\left(0\right)=0,D\left(1\right)=0]$. Notice that

$${U}_{{\nu}_{z}}(O;{\nu}_{z},\gamma ,\beta )=\frac{(1-D)\rho (0,X;\gamma )}{\rho (Z,X;\gamma )}(\mu (z,X;\beta )-{\nu}_{z})$$

has mean zero at ${\nu}_{z}^{\ast}$, *γ**, and *β**.

In some observational studies (e.g., NSCOT), individuals are probabilistically sampled. Let *S* denote the sampling indicator. Then, *O* is only observed if *S* = 1. An individual can be sampled with probability depending on a subset of *O*. Let *p*(*O*) denote the known, sampling probability. Then the estimating functions above will remain unbiased when they are multiplied by *S/p*(*O*). We will refer to *O*^{†} as the sampled data (i.e., *O* and *p*(*O*) for those individuals with *S* = 1) plus *S*. Let *U*_{γ} (*O*^{†}; ** γ**) =

Let ** ψ** = (

Using M-estimation theory (Huber, 1964), under mild regularity conditions, we can show

$$\sqrt{N}(\widehat{\psi}-{\psi}^{\ast})\stackrel{D}{\to}MVN(0,{\Sigma}^{\ast})$$

where ${\Sigma}^{\ast}=E{\left[\frac{\partial \mathit{U}({O}^{\u2020};\psi )}{\partial \psi}\right]}^{-1}E\left[\mathit{U}({O}^{\u2020};\psi )\mathit{U}{({O}^{\u2020};\psi )}^{\prime}\right]E{\left[\frac{\partial \mathit{U}({O}^{\u2020};\psi )}{\partial \psi}\right]}^{-{1}^{\prime}}$. We can use the empirical operator and substitute $\widehat{\psi}$ for ψ* to obtain an estimate of Σ*. Specifically, $\widehat{\psi}\approx MVN({\psi}^{\ast},{\widehat{\Sigma}}_{N})$ where

$${\widehat{\Sigma}}_{N}={\left[\sum _{i=1}^{N}\frac{\partial U({O}_{i}^{\u2020};\psi )}{\partial \psi}\right]}^{-1}\sum _{i=1}^{N}\left[U({O}_{i}^{\u2020};\psi )U{({O}_{i}^{\u2020};\psi )}^{\prime}\right]{\left[\sum _{i=1}^{N}\frac{\partial U({O}_{i}^{\u2020};\psi )}{\partial \psi}\right]}^{-{1}^{\prime}}$$

Since *SACE* is a differentiable function of *ψ**, we can use the delta method to obtain an approximation to the distribution of $\widehat{\mathit{SACE}}$.

Eligible patients in the NSCOT study were all trauma patients ages 18-84 who were treated at one of 18 trauma centers and 51 non-trauma centers for a moderately severe to severe injury between July 2001 and November 2002, and met certain criteria. Medical records were scheduled to be collected on all hospital deaths and a stratified sample of those discharged alive, with stratum-specific sampling probabilities. The strata were defined by hospital, age, injury severity score and principal body region injured. For those identified, not all medical records could be ascertained and for those that were ascertained, some patients were deemed ineligible. The sampling weights reflect the initial sampling scheme, the lack of ascertainment of records, and the ineligibility of ascertained records.

Medical records were ascertained on 1,391 of 1,438 hospital deaths. Of these 1,104 were identified as being eligible. Of the 16,760 individuals discharged alive from the hospital and initially identified as eligible, 8,021 were selected, and medical records were abstracted for 4,866 of these individuals. Of these, 4,087 were considered eligible. Of 5,191 total eligible with abstracted medical records, 148 were further excluded because they were either transferred to an NSCOT hospital greater than 24 hours post injury or had a length of stay prior to transfer from an NSCOT center of less than 24 hours, leaving a total of 5,043 eligibles. Three month post-injury interview data was collected on all 3,794 eligible individuals who were discharged alive and alive at the time of interview. A further 86 individuals died between 3 and 12 months, at which time a follow-up interview was planned. Of 3,708 alive at that time, 3,031 responded to the interview and had SF-36 vitality data.

In our analysis, *t** is 12 months, *Z* is indicator of trauma center care, *Y* is vitality status, and *D* is indicator of death by 12 months. Vitality, as defined by the vitality perceptions subscale of the SF-36 (Ware *et al.*, 1993), is a mental health measure of energy, pep, and tiredness that instrument developers have found to be correlated with both physical and mental function (McHorney *et al.*, 1993). High vitality scores indicate better outcomes.

The medical records include extensive information on pre-treatment potential confounders (*X*_{1}) such as age, race/ethnicity, gender, insurance type, number of pre-existing medical conditions, injury type and injury severity. These factors, discussed in detail in MacKenzie *et al.* (2006), include, for example, the continuous New Injury Severity Score (NISS) and injury subgroup (categorized into three groups: head, lower extremity, or abdominal/thoracic, spinal cord, upper extremity). The three month interview contains information on additional pre-injury confounders (*X*_{2}) including education, poverty status, self reported health status, difficulty with activities of daily living (ADL) and difficulty with instrumental activities of daily living (IADL). ADL measures an individual’s physical ability to perform tasks such as eating, dressing, and bathing. IADL measures an individual’s ability to care for oneself, such as being able to manage finances, prepare meals, and shop for oneself. Since these pre-injury factors are missing on those who died prior to a three month interview and are likely related to treatment, survival, and post-injury functioning, study investigators felt that, to prevent biased inferences, it was important to properly control for these factors. They were particularly concerned about pre-injury self rated health and difficulty with ADL and IADL; they were relatively less concerned about education and poverty status, as these factors were partially captured by race/ethnicity and insurance status from the medical records. With this in mind and with an eye toward a low-dimensional sensitivity analysis, we defined *X*_{2}* _{,Dep}* to be the indicator of fair or poor self-reported pre-injury health status

Before proceeding further, it is important to discuss the reasonableness of the identification assumptions. If we define *Z* = 1 (*Z* = 0) to be treatment at the closest, available trauma (non-trauma) center, then SUTVA (Assumption 1) is reasonable. One would not expect interference between individuals and the emphasis on closest, available facility-type rules out multiple versions of treatment. Monotonicity (Assumption 2) is reasonable since those treated in trauma centers have been shown to have similar or better mortality outcomes compared to those treated in non-trauma centers, as discussed by MacKenzie *et al.* (2006). This assumption could be violated if travel time and mode of transport to the nearest, available trauma and non-trauma centers are very different. However, investigators believe that there are likely few people for whom specialized care is harmful.

After controlling for medical and interview information, investigators are confident that they have controlled for the major confounding factors for evaluating the treatment effect on functioning. While it always possible that there are unmeasured confounders, they feel that the residual bias will be relatively negligible. While Assumption 4 (explainable non-random survival) seemed relatively plausible to the investigators, we will evaluate the robustness of inferences to deviations from this assumption below. Missing at random (Assumption 5) also seemed relatively reasonable to investigators, although one can imagine that, even after adjusting for *X*, those with poorer functioning might be less likely to respond to the interview. One can conduct sensitivity analysis to this assumption using inverse-weighting techniques described in Rotnitzky *et al.* (1998).

In our model for death, *h*(*Z*, *X*_{1}; *γ*) was taken to be additive in *X*_{1} with smooth functions (restricted cubic splines with 2 interior knots) of age and NISS, interactions of the main effects with *Z*, and interactions of the splines terms for age with injury subgroup. In our model for vitality status among survivors, *b*(*Z*, *X*; *β*) was taken to have similar functional dependence on *X*_{1} and *Z* as in the death model, with the exception that NISS was not included (since NISS has never been validated as an instrument for functional outcomes). Further, *b*(*Z*, *X*; *β*) was taken to be additive in *X*_{2} and include interactions with *Z*. To account for hospital-level clustering in the estimation of *SACE*, we used the robust standard error technique as described by Williams (2000). For data missing for reasons other than death, multiple imputation was used as described in MacKenzie *et al.* (2006, 2008).

In our analysis, we hypothesize that the sensitivity parameters *α _{Dep}* and

For comparison, we fit a naive linear model of SF-36 Vitality outcomes using *X*_{1} and *X*_{2} as covariates to obtain a naive estimate of the treatment effect among observed survivors.

In the sample, 60% of individuals were treated in trauma centers. From here on, we discuss weighted results. 28% and 16% of those interviewed at 3 months and treated in non-trauma and trauma centers, respectively, reported having fair or poor pre-injury health/functioning (*X*_{2,Dep} = 1). 12% and 10% of those treated in non-trauma and trauma centers, respectively, died by one year post-injury. Among those alive at one year who responded to the interview, the average SF-36 Vitality score was 50.02 and 52.08 for those treated in non-trauma and trauma centers, respectively.

Figure 1 presents the *SACE* analysis. The contour lines represent estimates of *SACE* as a function of *α _{Dep}* and

Estimates of *SACE* for the SF-36 Vitality subscale as a difference in expectations over ranges of the sensitivity parameters *α*_{Dep} and *τ*_{Dep} which describe the model main and trauma center interaction effect of the pre-injury health/functional **...**

In comparison, the estimate of the treatment effect from the naive regression analysis is 2.04 (95% CI -0.44, 4.52; p=0.11). It is not surprising that the naive estimate is attenuated. Such attenuation would occur when trauma centers tend to save the lives of individuals who would otherwise die in non-trauma centers and these individuals are destined to have lower vitality scores. *SACE* removes these individuals from the analysis by focusing only on the people who would live to one year regardless of the care received.

In this subsection, we consider a sensitivity analysis with respect to Assumption 4. Specially, we assume that

$$\frac{E[Y\left(1\right)\mid D\left(0\right)=1,D\left(1\right)=0,X]}{E[Y\left(1\right)\mid D\left(0\right)=0,D\left(1\right)=0,X]}=\kappa \left(X\right)$$

(7.1)

where *κ*(*x*) is a specified function of *x*. Here, we assume that this function is constant in *x*, i.e., *κ*(*x*) = κ, so that Assumption 4 is equivalent to *κ* = 1; *κ* > (<)1 implies that individuals who would die under non-trauma care but live under trauma care have higher (lower) mean functional outcomes under treatment than those who would survive regardless of whether or not assigned to treatment. In the NSCOT setting, it is likely that *κ* ≤ 1.

The identification formula for *E*[*Y* (0)|*D*(0) = 0, *D*(1) = 0] remains the same. For given *κ*, we have the following new identification formula for *E*[*Y* (1)|*D*(0) = 0, *D*(1) = 0]:

$$E[Y\left(1\right)\mid D\left(0\right)=0,D\left(1\right)=0]=\frac{E\left[\frac{\rho (0,X)\rho (1,X)}{\rho (0,X)+\kappa \{\rho (1,X)-\rho (0,X)\}}\mu (1,X)\right]}{E\left[\rho (0,X)\right]}$$

(7.2)

The estimating functions remain the same, except we re-define

$${U}_{{\nu}_{1}}(O;{\nu}_{1},\gamma ,\beta )=\frac{(1-D)\rho (0,X;\gamma )}{\rho (Z,X;\gamma )}\left(\frac{\rho (1,X;\gamma )}{\rho (0,X;\gamma )+\kappa \{\rho (1,X;\gamma )-\rho (0,X;\gamma )\}}\mu (1,X;\beta )-{\nu}_{1}\right)$$

When we performed sensitivity analysis with respect to *κ* (*κ* = 0.25, 0.50, 0.75) for selected combinations of *α _{Dep}* and

Overall, our analysis indicates that trauma centers have a statistically significant impact on vitality as measured by the SF-36 score. The magnitude of the effect is modest as changes in SF-36 scores of 5 or more have been cited as clinically relevant (Kahn *et al.*, 2005).

We evaluated the finite sample behavior of our estimation procedure when a key confounder was missing among those who die. We considered a case with two confounding covariates, *X*_{1} (continuous; fully observed) and *X*_{2} (binary; missing on those who die). For our simulation, we generated data (under Assumptions 1-6) for an individual according to the following scheme: (1) generate *X*_{1} as Normal with mean 0 and variance 1; (2) generate *X*_{2} as Bernoulli with *P*[*X*_{2} = 1|*X*_{1}] = expit(0.5 + 0.25*X*_{1}); (3) generate *Z* as Bernoulli with *P*[*Z* = 1|*X*_{1},*X*_{2}] = expit(-2 + 2.5*X*_{1} + 3.5*X*_{2}); (4) generate *D*(0) as Bernoulli with *P*[*D*(0) = 1|*X*_{1},*X*_{2}] = expit(-0.5 + 2.5*X*_{1} + 3*X*_{2}); (5) if *D*(0) = 0 set *D*(1) = 0; otherwise, generate *D*(1) as Bernoulli with *P*[*D*(1) = 1|*D*(0) = 1, *X*_{1}, *X*_{2}] = expit(-2.25 + *X*_{1} + 1.75*X*_{2}); (6) if *D*(0) = 0, generate *Y* (0) as Normal with *E*[*Y* (0)|*D*(0) = 0, *X*_{1},*X*_{2}] = 40 + 5*X*_{1} + 40*X*_{2} and variance 1; (7) if *D*(1) = 0 generate *Y* (1) as Normal with *E*[*Y*(1)|*D*(1) = 0,*X*_{1},*X*_{2}] = 33 + *X*_{1} + 20*X*_{2}, and variance 1; (8) set *D* = *D*(*Z*) and *Y* = *Y* (*Z*) if *D* = 0.

Using the analysis model logit*P*[*D* = 1|*Z*,*X*_{1},*X*_{2}] = *β*_{0} + *β*_{1}*Z* + *β*_{2}*X*_{1} + *β*_{3}*X*_{1}*Z* + *αX*_{2} + *τX*_{2}*Z*, the true values of *α* and *τ* are approximately 3.0 and -1.0. Since these values are unknown when *X*_{2} is missing on those who die, we conducted a sensitivity analysis, where we considered combinations of these parameters in a 4 × 4 square, centered at the truth. We used a linear (analysis) model for the conditional (on *Z*, *X*_{1}, and *X*_{2}) mean of *Y* among observed survivors, where the functional form of the linear predictor was the same as that used in death model above. For fixed *α* and *τ*, we estimated *SACE* and derived confidence intervals using the large sample theory approximation for the variance as described in Section 6. We evaluated the proportion of confidence intervals that contained the true value of *SACE* = -10.4.

In our data generation scheme, higher values of *Y* (0) and *Y* (1) were considered worse outcomes and higher values of *X*_{1} and *X*_{2} were indicative of poorer health and worse outcomes. Those with poorer health were more likely to be assigned to the treatment group. The probability of being assigned to treatment was approximately 53%. The parameters resulted in approximately 53% of the sample being assigned to *Z* = 1, and 62% having *X*_{2} = 1. In Table 1, we present treatment-specific distributional summaries of key variables.

The unadjusted naive effect of the intervention indicates that the intervention is harmful (*E*[*Y*|*D* = 0, *Z* = 1] - *E*[*Y*|*D* = 0, *Z* = 0] = 2.8). Adjusting for *X*_{1} and *X*_{2} using an analysis model for *E*[*Y*|*D* = 0,*Z*,*X*_{1},*X*_{2}] which is linear in *Z*, *X*_{1}, and *X*_{2} and includes no interactions between *Z* and (*X*_{1}, *X*_{2}), yields an overly optimistic treatment effect of -16.7. The effect among the survivors under each treatment (*E*[*Y* (1)|*D*(1) = 0] - *E*[*Y* (0)|*D*(0) = 0]) is -5.99, a more modest effect. This example demonstrates the importance of estimating *SACE* even when observed death rates between the treatment arms are relatively similar (38% vs. 43%).

For each of three sample sizes, 500, 1000, and 2000, we generated 1000 simulated datasets. Figure 2 displays the results. Coverage of the 95% confidence intervals was good when *α* and *τ* are fixed at their true values: 94.9%, 93.4%, and 93.7% for the three respective sample sizes. The coverage worsened when *α* and *τ* were fixed at values different than the truth. Coverage was particularly poor when *α* was assumed to be smaller than the truth. In terms of bias, the trends were very similar: there was no evidence of bias when the sensitivity analysis parameters were correctly specified and absolute bias increases substantially when *α* was assumed to be smaller than the truth.

Results of simulations with the 95% coverage indicated using contour lines and the absolute bias of the point estimates represented by shading. The figures are centered at the true values of the sensitivity parameters, denoted *α** and *τ* **...**

The top left corner of the figures depicts the region in which *α* and *τ* are incorrectly assumed to equal zero. It is in this region where the bias and empirical coverage are worst. Specifically, the results are biased towards a more beneficial treatment benefit. The reason that the treatment estimates are more beneficial as *α* decreases is that such an assumption reduces the probability that those with high values of *X*_{2} will die. Since higher values of *X*_{2} are associated with a greater treatment effect, the point estimates of the intervention become stronger as more individuals with higher values of *X*_{2} are assumed to survive.

Our analysis model for the conditional probability of death, logit*P*[*D* = 1|*Z*,*X*_{1},*X*_{2}], might be incorrectly specified in this simulation. The true values of ** α** and

Overall, the simulations suggest the utility of the sensitivity analysis approach. There is little bias and good coverage when the sensitivity parameters are correctly specified. Since the true value of the sensitivity parameters are not identifiable without additional assumptions and/or data, it is important to present inferences over a range of these parameters.

This paper has presented a sensitivity analysis methodology for estimating *SACE* in observational studies, where key confounders are missing due to death. Our analysis of the NSCOT study demonstrated that our estimate of *SACE* is larger than a naive regression analysis that conditions on observed survivors. Further, our analysis was insensitive to the choices of the sensitivity analysis parameters. As our bias analysis and simulation study demonstrated, such insensitivity is by no means guaranteed.

MacKenzie *et al.* (2008) utilized our methodology to investigate the effect of trauma center care on health status outcomes of NSCOT patients with major lower limb trauma.

In our analyses of the NSCOT data, we did not utilize 3 month interview information (i.e., *X*_{2}) on the 86 individuals who died between the 3 and 12 month interviews. Since the number of individuals is small, the influence on our analysis is likely to be negligible. It is straightforward to generalize our methodology to utilize this extra information by appropriately enriching the data structure, models, and estimating functions.

Research supported in part by NIH grants P30 CA 06927, DA 009897, MH 47447, grants R49/CCR316840 and R49CE000198 from the National Center for Injury Prevention and Control of the Centers for Disease Control and Prevention, grant P30 MH066247 from the Center for Prevention and Early Intervention, and an appropriation from the Common-wealth of Pennsylvania

- Egleston BL, Scharfstein DO, Freeman EE, West SK. Causal inference for non-mortality outcomes in the presence of death. Biostatistics. 2007;8(3):526–545. [PubMed]
- Frangakis CE, Rubin DB, An M-W, MacKenzie E. Principal stratification designs to estimate input data missing due to death. Biometrics. 2007;63(3):641–649. [PubMed]
- Hayden D, Pauler DK, Schoenfeld D. An estimator for treatment comparisons among survivors in randomized trials. Biometrics. 2005;61(1):305–310. [PubMed]
- Huber PJ. Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics. 1964;35(1):73–101.
- Kahn SR, Ducruet T, Lamping DL, Arsenault L, Miron MJ, Roussin A, Desarais S, Joyal F, Kassis J, Solymoss S, Desjardins L, Jori M, Shrier I. Prospective evualation of health-related quality of life in patients with deep venous thrombosis. Archives of Internal Medicine. 2005;165:1173–1178. [PubMed]
- MacKenzie EJ, Rivara FP, Jurkovich GJ, Nathens AB, Frey KP, Egleston BL, Salkever DS, Scharfstein DO. A national evaluation of the effect of traumacenter care on mortality. New England Journal of Medicine. 2006;354:366–378. [PubMed]
- MacKenzie EJ, Rivara FP, Jurkovich GJ, Nathens AB, Frey KP, Egleston BL, Salkever DS, Scharfstein DO. The impact of trauma center care on one-year functional outcomes following major lower limb trauma. Journal of Bone and Joint Surgery. 2008;90:101–109. [PubMed]
- McHorney CA, Ware JE, Raczek AE. The MOS 36-item Short Form Health Survey (SF-36): Psychometric and clinical tests of validity in measuring physical and mental health constructs. Medical Care. 1993;31(3):247–263. [PubMed]
- Robins JM, Greenland S. Comment to causal inference without counterfactuals. Journal of the American Statistical Association. 2000;95:431–435.
- Robins JM, Ritov Y. Toward a curse of dimensionality appropriate (CODA) asymptotic theory for semi-parametric models. Statistics in Medicine. 1997;16:285–319. [PubMed]
- Rosenbaum PR. The consequences of adjustment for a concomitant variable that has been affected by the treatment. Journal of the Royal Statistical Society, Series A. 1984;147:656–666.
- Rotnitzky A, Robins JM, Scharfstein DO. Semiparametric regression for repeated outcomes with non-ignorable non-response. Journal of the American Statistical Association. 1998;93:1321–1339.
- Rubin DB. Comment on Randomization Analysis of Experimental Data: The Fisher Randomization Test, by D. Basu. Journal of the American Statistical Association. 1980;75:591–593.
- Rubin DB. Comment to causal inference without counterfactuals. Journal of the American Statistical Association. 2000;95:435–438.
- Ware JE, Snow KK, Kosinski M, Gandek B. SF-36® Health Survey Manual and Interpretation Guide. New England Medical Center, The Health Institute; Boston, MA: 1993.
- Williams RL. A note on robust variance estimation for cluster-correlated data. Biometrics. 2000;56:645–646. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library 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. |