On Estimation of the Survivor Average Causal Effect in Observational Studies when Important Confounders are Missing Due to Death

doi:10.1111/j.1541-0420.2008.01111.x

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

Published in final edited form as:

Biometrics. 2009 June; 65(2): 497–504.

doi: 10.1111/j.1541-0420.2008.01111.x

PMCID: PMC2700847

NIHMSID: NIHMS52102

On Estimation of the Survivor Average Causal Effect in Observational Studies when Important Confounders are Missing Due to Death

Brian L. Egleston,¹ Daniel O. Scharfstein,² and Ellen MacKenzie³

¹Biostatistics Facility, Fox Chase Cancer Center, 333 Cottman Avenue, Philadelphia, PA 19111-2497

²Department of Biostatistics, Johns Hopkins School of Public Health

³Department of Health Policy and Management, Johns Hopkins School of Public Health

Ph: 215-214-3917, e-mail: Brian.Egleston/at/fccc.edu

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

See other articles in PMC that cite the published article.

Abstract

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.

1 Introduction

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.

2 Data Structure and Notation

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₁, X₂} denote a set of confounding covariates, where X₁ is observed and X₂ 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₁, (X₂ : D = 0), (Y : D = 0, R = 1)}

3 Causal Estimand

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,

4 Assumptions

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:

Assumption 1

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.

Assumption 2

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.

Assumption 3

{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.

Assumption 4

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.

Assumption 5

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.

5 Identification, Bias, and Sensitivity Analysis

5.1 Identifiability

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

(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₂ 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₂ is not observable on those who die.

The key to identification of SACE when X₂ is missing is the following key identity:

(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,

(5.3)

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

Assumption 6

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

Here, we assume no interaction between x₁ and x₂. However, additional sensitivity analysis parameters could be added to the model to account for interactions between x₁ and x₂. 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₁) is identified from the observed data. It can be shown that

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 α = τ = 0, Assumption 6 implies that X₂ is independent of D(z) given X₁, for z = 0, 1. In some circumstances, this assumption might be reasonable. For example, when X₂ records sociodemographic factors, then one might be able to argue that these factors are not related to the potential mortality outcomes, after conditioning on health status factors recorded in X₁.

5.2 Bias

What is the implication of falsely assuming that α = τ = 0? To answer this question, consider a stratum of X₁, say x₁. Define P_x₁ [A|B] and E_x₁ [A|B] as the conditional probability and expectation of event A given X₁ = x₁ and event B, respectively. Now, suppose that X₂ is dichotomous, where subjects with X₂ = 1 are sicker than those with X₂ = 0. Thus, P_x₁[D(0) = 0|X₂ = 1] < P_x₁[D(0) = 0|X₂ = 0]. In general, we can write the x₁-specific SACE estimand, SACE_x₁, as

(5.4)

When α = τ = 0, the above expression reduces to

(5.5)

The right hand sides of (5.4) and (5.5) will be unequal since P_x₁[D(0) = 0|X₂ = 1] < P_x₁[D(0) = 0|X₂ = 0]. If 0 < P[X₂ = 1] < 1 and E_x₁[Y(1) - Y(0)|D(0) = 0, X₂ = 1] > (<)E_x₁[Y(1) - Y(0)|D(0) = 0, X₂ = 0], then falsely assuming α = τ = 0 will lead to an over(under)-estimate of the stratum-specific estimand. If P_x₁[D(0) = 0|X₂ = 1] < P_x₁[D(0) = 0|X₂ = 0] for all x₁ (i.e., α > 0) and E_x₁[Y(1) - Y(0)|D(0) = 0, X₂ = 1] is larger (smaller) than E_x₁[Y(1) - Y(0)|D(0) = 0, X₂ = 0] for all x₁, then there will be over(under)-estimation of SACE when α = τ = 0 is falsely assumed.

5.3 Sensitivity Analysis

Since α and τ are not identifiable, we recommend conducting inference about SACE over plausible ranges of these sensitivity analysis vectors. A major limitation of this approach is that when α and τ have collectively greater than, say 4, components, it becomes difficult to display and synthesize the results of the sensitivity analysis. Thus, we recommend that (1) components of X₂ be replaced by a low-dimensional summary and/or (2) additional conditional independence assumptions be imposed by setting some components of α and τ equal to zero. Specifically, we partition X₂, α and τ so that equation M8

, where the Dep and Ind subscripts refer to the components of X₂ that are conditionally dependent and independent, respectively.

6 Estimation

In order to estimate SACE from a dataset of size N, we need to estimate μ(z, x) and h(z, x₁). When X₁ 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 An external file that holds a picture, illustration, etc.
Object name is nihms-52102-ig0001.jpg

An external file that holds a picture, illustration, etc.
Object name is nihms-52102-ig0001.jpg

-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₁). In particular, we assume that h(z, x₁) = h(z, x₁; γ*) and g(μ(z, x)) = b(z, x; β*). where h(z, x₁; γ) is a specified function of z, x₁ 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 τ. We define a(Z, X; γ) = h(Z, X₁; γ) + α’X₂ + τ’Z · X₂, ρ(Z, X; γ) = 1 - expit{a(Z, X; γ)}, and μ(Z, X; β) = g^-1{b(Z, X; β)}.

6.1 Estimating Functions

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

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

has mean zero when evaluated at β*, where b’ (Z, X; β) is the derivative of b(Z, X₁; β) with respect to β. Further, we can set up an unbiased estimating function for equation M13

. Notice that

has mean zero at equation M15

, γ*, 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^†; γ) = SU_γ(O; γ)/p(O), U_β(O^†; β) = SU_β(O; β)/p(O), U_{ν_z}(O^†; ν_z, γ, β) = SU_{ν_z}(O; ν_z, γ, β)/p(O).

Let ψ = (γ’, β’, ν₀, ν₁)’ and ψ* denote the true value of ψ. Further, define U (O^†; ψ) = (U_γ(O^†; γ)’, U_β(O^†; β)’, U_ν₀(O^†; γ, β, ν₀), U_ν₁(O^†; γ, β, ν₁). Then, ψ* is estimated by equation M16

, as the solution to equation M17

, where N denotes the size of the population from which individuals were probabilistically drawn. Since SACE is a function of equation M18

and

, we can then estimate SACE, by equation M20

6.2 Large Sample Theory

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

where

. We can use the empirical operator and substitute equation M23

for ψ* to obtain an estimate of Σ*. Specifically, equation M24

where

Since SACE is a differentiable function of ψ*, we can use the delta method to obtain an approximation to the distribution of equation M26

7 NSCOT Study

7.1 Features

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₁) 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₂) 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₂_,Dep to be the indicator of fair or poor self-reported pre-injury health status or one or more pre-injury ADL or IADL difficulties, and X₂_,Ind to include income (low - reference, middle, high) and education (less than high school-reference, high school, some college or more).

7.2 Assumptions

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).

7.3 Models

In our model for death, h(Z, X₁; γ) was taken to be additive in X₁ 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₁ 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₂ 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 τ_Dep lie between 0 and 2.7 and between -2.7 and 2.7, respectively, subject to the constraint that 0 ≤ α_Dep + α_Dep ≤ 2.7. These constraints imply that, within levels of X₁, the odds of dying under trauma (non-trauma) center care is the same or at most 15 times greater for patients with pre-injury health/function problems that those who do not. This incorporates an extreme assumption about the magnitude of the effects.

For comparison, we fit a naive linear model of SF-36 Vitality outcomes using X₁ and X₂ as covariates to obtain a naive estimate of the treatment effect among observed survivors.

7.4 Results

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 τ_Dep. The shading in the figure represent the associated p-values. When α_Dep and τ_Dep are both equal to 0, mean vitality score among those who would survive to one-year regardless of whether or not assigned to treatment is 51.94 (95% CI 50.57, 53.33) under trauma center care and 48.90 (95% CI 47.25, 50.56) when under non-trauma center care. This gives an estimate of SACE of 3.04 (95% CI 0.92, 5.17; p=0.005), suggesting that trauma center care improves vitality outcomes. As the figure demonstrates, the results are highly insensitive to choice of α_Dep and τ_Dep.

Figure 1

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 (more ...)

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.

7.5 Sensitivity Analysis for Assumption 4

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

(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]:

(7.2)

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

When we performed sensitivity analysis with respect to κ (κ = 0.25, 0.50, 0.75) for selected combinations of α_Dep and τ_Dep, we found that inferences about SACE were generally similar to those reported in Figure 1. For example, when α_Dep = 0 and τ_Dep = 2, the estimated value of SACE is 2.97 at κ = 0.75 and 3.10 at κ = 0.50 and the results are statistically significant. Only when α_Dep is highly positive and τ_Dep highly negative did the results become marginally statistically significant. For example, when α_Dep = 2.7, τ_Dep = -2.7, and κ = 0.50, then the estimated value of SACE is 2.22 (p = 0.06).

7.6 Summary

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).

8 Simulation Study

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₁ (continuous; fully observed) and X₂ (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₁ as Normal with mean 0 and variance 1; (2) generate X₂ as Bernoulli with P[X₂ = 1|X₁] = expit(0.5 + 0.25X₁); (3) generate Z as Bernoulli with P[Z = 1|X₁,X₂] = expit(-2 + 2.5X₁ + 3.5X₂); (4) generate D(0) as Bernoulli with P[D(0) = 1|X₁,X₂] = expit(-0.5 + 2.5X₁ + 3X₂); (5) if D(0) = 0 set D(1) = 0; otherwise, generate D(1) as Bernoulli with P[D(1) = 1|D(0) = 1, X₁, X₂] = expit(-2.25 + X₁ + 1.75X₂); (6) if D(0) = 0, generate Y (0) as Normal with E[Y (0)|D(0) = 0, X₁,X₂] = 40 + 5X₁ + 40X₂ and variance 1; (7) if D(1) = 0 generate Y (1) as Normal with E[Y(1)|D(1) = 0,X₁,X₂] = 33 + X₁ + 20X₂, and variance 1; (8) set D = D(Z) and Y = Y (Z) if D = 0.

Using the analysis model logitP[D = 1|Z,X₁,X₂] = β₀ + β₁Z + β₂X₁ + β₃X₁Z + αX₂ + τX₂Z, the true values of α and τ are approximately 3.0 and -1.0. Since these values are unknown when X₂ 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₁, and X₂) 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₁ and X₂ 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₂ = 1. In Table 1, we present treatment-specific distributional summaries of key variables.

Table 1

Distributional summaries of simulated variables, stratified by treatment assignment

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₁ and X₂ using an analysis model for E[Y|D = 0,Z,X₁,X₂] which is linear in Z, X₁, and X₂ and includes no interactions between Z and (X₁, X₂), 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.

Figure 2

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 τ (more ...)

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₂ will die. Since higher values of X₂ are associated with a greater treatment effect, the point estimates of the intervention become stronger as more individuals with higher values of X₂ are assumed to survive.

Our analysis model for the conditional probability of death, logitP[D = 1|Z,X₁,X₂], might be incorrectly specified in this simulation. The true values of α and τ are those that minimize the Kullback-Leibler distance between the true and (possibly) misspecified model. They are found by simulating a massive dataset and estimating the value of these parameters via maximum likelihood. The coverage probabilities and bias at these values of α and τ suggest that our analysis model works well. Further, when we performed simulations using a more flexible analysis model that entered X₁ into the model using restricted cubic splines, our results did not substantively change. This suggests that misspecification of the analysis model did not affect our simulation inferences.

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.

9 Discussion

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₂) 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.

10 Acknowledgements

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

11 References

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]