Assessing the sensitivity of regression results to unmeasured confounders in observational studies

D. Y. Lin^{1}, B. M. Psaty^{2}, and R. A. Kronmal^{1}

1Department of Biostatistics, Box 357232, University of Washington, Seattle, WA 98195, USA 2Departments of Medicine, Epidemiology and Health Services, University of Washington, Seattle,

WA 98195, USA

SUMMARY

This paper presents a general approach for assessing the sensitivity of the point and interval estimates of the primary treatment or exposure effect in an observational study to the residual confounding effects of unmeasured variables after adjusting for measured covariates. The proposed method assumes that the true treatment effect can be represented in a regression model which includes the treatment indicator as well as the measured and unmeasured confounders. Under this formulation, we can use the corresponding reduced model which omits the unmeasured confounder to make statistical inferences about the true treatment effect by specifying the distributions of the unmeasured confounder in the two treatment groups along with the possibly treatment-dependent effects of the unmeasured confounder on the outcome variable. If the unmeasured and measured confounders are independent within each treatment group, there exists a simple algebraic relationship between the true treatment effect under the full model and the apparent treatment effect under the reduced model. We can then estimate the true treatment effect by making a simple adjustment to the point and interval estimates of the apparent treatment effect obtained from a software package or reported in a published study. The proposed method handles both binary response and censored survival time data, accommodates any study design, and allows the unmeasured confounder to be discrete or normally distributed. We describe applications to two major medical studies.

Key words: Case-control studies; Causal inference; Cohort studies; Covariate adjustment; Cox regression; Epidemiology; Logistic regression; Model misspecification; Non-randomized studies; Omitting covariate; Sensitivity analysis.

1. INTRODUCTION

Well-conducted randomized controlled experiments have been widely recognized as the decisive method for assessing treatments, exposures or policies. Due to ethical, economic and other constraints, randomized experiments are not feasible in many circumstances. Instead, observational studies, including cohort and case-control studies, are commonly undertaken. Because patients in observational studies are not assigned to treatments at random, the apparent association between treatment and outcome may result from some unknown or unmeasured covariates that are associated with both treatment and outcome. Thus, it would be desirable to investigate how the findings of an observational study may be affected by certain variations in the assumptions about the unmeasured confounders. If the conclusions are insensitive over a wide range of plausible assumptions, then the number of interpretations of the observational data is reduced, and the causal conclusions become more defensible.

This idea of sensitivity analysis can be traced back to Cornfield et al. (1959), who responded to Fisher (1958)’s constitutional hypothesis that the apparent association between cigarette smoking and lung cancer could be explained by the confounding effects of a genotype that predisposed to both smoking and lung cancer. A centerpiece of their argument was a theorem due to Cornfield: for a confounder to explain a relative risk of a given magnitude, say r, this confounder has to be r times more prevalent among the exposed than among the unexposed. Because the lung cancer relative risk was approximately 10 for cigarette smokers versus nonsmokers and increased to 20 fold for heavy smokers, the existence of such a confounder seemed implausible. This line of statistical reasoning played a critical role in the 1964 Surgeon General’s Report concluding that “Cigarette smoking is causally related to lung cancer in men”.

Bross (1966) and Schlesselman (1978) obtained results similar to Cornfield’s inequality. Rosenbaum (1995, Ch. 4) provided the bounds for the significance levels of several common two-sample tests and for Hodges-Lehmann estimates. In addition, Rosenbaum and Rubin (1983a) discussed the sensitivity analysis on a point estimate for the difference between two success probabilities in a 2 × 2 × S contingency table. The method of Rosenbaum and Rubin (1983a) specifies stratumspecific sensitivity parameters, consisting of the overall prevalence of the binary confounder, the confounder’s effect on the selection of treatment and the confounder’s effect on the binary outcome, and provides a range of point estimates corresponding to various combinations of the sensitivity parameters. By contrast, the methods of Cornfield et al. (1959) and Rosenbaum (1995) specify only a single sensitivity parameter and give the corresponding bound for the treatment effect which caters to the least favorable or most conservative values of the remaining parameters.

In this paper, we develop a general yet simple method of sensitivity analysis, which shares the spirit of the published methods, especially that of Rosenbaum and Rubin (1983a). Formulating the problem more broadly than Rosenbaum and Rubin, we provide an alternative solution that applies to a variety of statistical models and data structures. Specifically, we assume that the treatment effect and the effects of measured and unmeasured confounders can be formulated through a regression model. This formulation implies a second model which excludes the unmeasured confounder. It is possible to estimate the true treatment effect under the full model by fitting the reduced model to the observable data provided that one specifies the distribution of the unmeasured confounder in each treatment arm as well as the possibly treatment-dependent effects of the unmeasured confounder on the outcome variable. As a matter of fact, given these two sets of sensitivity parameters, there exists a simple algebraic relationship between the true treatment effect under the full model and the apparent treatment effect under the corresponding reduced model if the unmeasured and measured confounders are independent within each treatment group. In these situations, we can make inferences about the true treatment effect by making simple adjustments to the point and interval estimates of the apparent treatment effect, which are obtained either from the standard software package or from the published paper.

The proposed method offers several advantages over the existing ones. First, it applies to both binary response and censored survival time data, and to the corresponding logistic and proportional hazards models. Survival time or failure time is an important outcome variable in many observational cohort studies, but this endpoint has not been studied in the sensitivity analysis literature except for a brief mention of some matched-pair tests by Rosenbaum (1995, p. 115). Secondly, the proposed method allows adjustment for any type of measured confounders, discrete or continuous, whereas the published methods can at most accommodate stratification on categorical variables. Thirdly, the unmeasured confounders may be binary or normally distributed; the existing literature is confined to a bounded confounder. Finally, we provide simple and explicit formulas to adjust both the point and interval estimates of the apparent treatment effect in the reduced model. These formulas hold for both binary and survival endpoints and are invariant across study designs.

2. BINARY OUTCOME

Let Y be a binary response variable which takes the value 1 if the event of interest, such as disease or death, occurs, and let X be the exposure indicator taking the value 1 for the exposed. Also, let Z be a set of measured covariates, and let U be an unmeasured confounder. We first assume that the response probability is related to X, Z and U through the log-linear model:

Pr(Y = 1|X, Z, U) = exp (α + βX + γ_{X}U + θ′Z) ,  (2.1)

where α, β, γ_{X} (X = 0, 1) and θ are unknown regression parameters. Note that the parameterization βX + γ_{X}U is the same as βX + γ_{0}U + (γ_{1} − γ_{0})XU. Since U is unmeasured, one is forced to fit the reduced model

Pr(Y = 1|X, Z) = exp (α^{∗} + β^{∗}X + θ^{∗′}Z) ,  (2.2)

where α^{∗}, β^{∗} and θ^{∗} are potentially different from α, β and θ of model (2.1). We refer to β and β^{∗} as the true and apparent exposure effects, respectively. Since β^{∗} can be directly estimated from the observable data whereas β can not, it is of great interest to ascertain the relationship between β^{∗} and β.

Model (2.1) is slightly different from the logistic regression model

Pr(Y = 1|X, Z, U) = \frac{exp (α + βX + γ_{X}U + θ′Z)}{1 + exp (α + βX + γ_{X}U + θ′Z)},  (2.3)

i.e.,

logit{Pr(Y = 1|X, Z, U)} = α + βX + γ_{X}U + θ′Z.

For rare events, the right side of (2.1) is a good approximation to that of (2.3). If Z is excluded from models (2.1) and (2.3) and U is binary, then these two models are saturated and do not impose any restrictions on the effects of X and U on Y. In this case, models (2.1) and (2.3) are equivalent. In the presence of Z, model (2.1) is mathematically more tractable than (2.3). Thus, it is instructive to consider model (2.1) first.

2.1. Log-Linear Regression

Let F (u|X, Z) be the distribution function of U given X and Z. By the law of conditional expectation,

Pr(Y = 1|X, Z) = ∫_{−∞}^{∞} Pr(Y = 1|X, Z, u)dF (u|X, Z),  (2.4)

which becomes

Pr(Y = 1|X, Z) = exp (α + βX + θ′Z) ∫_{−∞}^{∞} e^{γ_{X}u}dF (u|X, Z)  (2.5)

under model (2.1). We will evaluate the right side of the above equation for the cases of binary and normally distributed U.

Binary confounder. Suppose that U is binary such that F (u|X, Z) is a Bernoulli distribution with success probability P_{X,Z}. Then equation (2.5) becomes

Pr(Y = 1|X, Z) = exp (α + βX + θ′Z) {e^{γ_{X}} P_{X,Z} + (1 − P_{X,Z})} .  (2.6)

For simplicity, we assume that U is independent of Z conditional on X, which implies that P_{X,Z} does not depend on Z. Under P_{X,Z} = P_{X}, equation (2.6) may be written as

Pr(Y = 1|X, Z) = exp [α + log \{e^{γ_{0}}P_{0} + (1 − P_{0})\} + \{β + log \frac{e^{γ_{1}}P_{1} + (1 − P_{1})}{e^{γ_{0}}P_{0} + (1 − P_{0})}\} X + θ′Z] .

(2.7)

Comparing the right sides of (2.2) and (2.7), we find that

β = β^{∗} − log \frac{e^{γ_{1}}P_{1} + (1 − P_{1})}{e^{γ_{0}}P_{0} + (1 − P_{0})}  (2.8)

or

R = R^{∗}/ \{\frac{Γ_{1}P_{1} + (1 − P_{1})}{Γ_{0}P_{0} + (1 − P_{0})}\} ,  (2.9)

where R = e^{β}, R^{∗} = e^{β^{∗}}, Γ_{1} = e^{γ_{1}} and Γ_{0} = e^{γ_{0}}. Note that R and R^{∗} are, respectively, the true and apparent relative risks of disease associated with X, while Γ_{1} and Γ_{0} are the relative risks of disease associated with U among the exposed and unexposed, respectively. Formula (2.9) shows succinctly how the prevalences of the unmeasured confounder, P_{1} and P_{0}, and its relative risks of disease, Γ_{1} and Γ_{0}, among the exposed and unexposed determine the relationship between the true and apparent relative risks of disease associated with the exposure.

Define

A = \frac{Γ_{1}P_{1} + (1 − P_{1})}{Γ_{0}P_{0} + (1 − P_{0})}.

It is of most interest to examine A when P_{1} > P_{0} and Γ_{1} ≥ Γ_{0} > 1. In such a setting, the adjustment factor A is bounded above by Γ_{1}; this bound is achieved when P_{0} = 0 and P_{1} = 1. Under Γ_{0} = Γ_{1} = Γ, A is also bounded above by P_{1}/P_{0}; it approaches this bound as Γ → ∞. In applications, it is more informative to evaluate A for specific values of P_{0}, P_{1}, Γ_{0} and Γ_{1} rather than using the upper bounds. Table 1 displays the actual values of A under Γ_{0} = Γ_{1} = Γ = 2, 3, 4, 6, 8, 10 and for P_{0} and P_{1} taking values from 0 to 1 with 0.1 increments.

Normal confounder. Suppose that, conditional on X and Z, the confounder U is normally distributed with mean μ_{X,Z} and unit variance. The standardization of U entails that γ_{0} and γ_{1} are standardized regression coefficients. A simple integral calculation yields

∫_{−∞}^{∞} e^{γ_{X}u} \frac{1}{\sqrt{2π}} exp \{−\frac{(u − μ_{X,Z})^{2}}{2} \} du = exp (γ_{X}μ_{X,Z} + 0.5γ_{X}^{2}) ,

which entails that equation (2.5) is

Pr(Y = 1|X, Z) = exp (α + βX + θ′Z + γ_{X}μ_{X,Z} + 0.5γ_{X}^{2}) .  (2.10)

As with the binary confounder, considerable simplification arises when U is conditionally independent of Z given X. Under this condition with μ_{X} = μ_{X,Z}, equation (2.10) can be written as

Pr(Y = 1|X, Z) = exp [α + γ_{0}μ_{0} + 0.5γ_{0}^{2} + {β + (γ_{1}μ_{1} − γ_{0}μ_{0}) + 0.5(γ_{1}^{2} − γ_{0}^{2})} X + θ′Z] .

Thus,

β = β^{∗} − {(γ_{1}μ_{1} − γ_{0}μ_{0}) + 0.5(γ_{1}^{2} − γ_{0}^{2})} ,  (2.11)

or

R = R^{∗}/ exp{(γ_{1}μ_{1} − γ_{0}μ_{0}) + 0.5(γ_{1}^{2} − γ_{0}^{2})}.  (2.12)

Further simplification arises when γ_{0} = γ_{1} = γ, in which case (2.11)–(2.12) reduce to

β = β^{∗} − γδ,  (2.13)
R = R^{∗}/Γ^{δ},  (2.14)

where δ = μ_{1} − μ_{0}. Under γ_{0} = γ_{1}, equations (2.13)–(2.14) hold even if U is not conditionally independent of Z, but the effects of X and Z on μ_{X,Z} must be additive, i.e., μ_{X,Z} = μ_{X} + q(Z), where q is some arbitrary function of Z. Note that the adjustment factor Γ^{δ} in (2.14) depends on the difference between μ_{1} and μ_{0}, not on their actual values.

Formulas (2.8)–(2.9) and (2.11)–(2.14) enable one to conduct the sensitivity analysis in an extremely simple fashion. The sensitivity parameters consist of (P_{0}, P_{1}, Γ_{0}, Γ_{1}) and (μ_{0}, μ_{1}, Γ_{0}, Γ_{1}) in the binary and normal cases, respectively. Since β^{∗} and consequently R^{∗} can be directly estimated by fitting model (2.2) to the observed data on (Y, X, Z), we can use formulas (2.8)–(2.9) and (2.11)– (2.14) to estimate β and R upon specifying the sensitivity parameters. Denote the estimators of β^{∗}, β, R^{∗} and R by β̂^{∗}, β̂, R̂^{∗} and R̂. Equations (2.8) and (2.11) show that, when the sensitivity parameters are fixed at specific values, the variance of β̂ is the same as that of β̂^{∗}. The latter can be estimated by the standard estimation theory. The fact that, for fixed values of the sensitivity parameters, the estimators β̂ and β̂^{∗} have the same variance implies that the adjustment factors for the limits of the confidence intervals for β are the same as that of the point estimate. In other words, formulas (2.8)–(2.9) and (2.11)–(2.14) apply to both the point and interval estimates of β. Furthermore, all these results hold for any study design: prospective or retrospective, matched or unmatched.

If the conditional independence of U and Z given X is not plausible, then there is generally no explicit relationship between β and β^{∗}. Nevertheless, when the sensitivity parameters (γ_{X}, P_{X,Z}) and (γ_{X}, μ_{X,Z}) are specified in the binary and normal cases, respectively, the right sides of models (2.6) and (2.10) involve only the measured variables X and Z and the corresponding regression parameters α, β and θ. One may then obtain β̂ and its variance estimator by applying the maximum likelihood method or other requisite methods to models (2.6) and (2.10).

2.2. Logistic Regression

Under the logistic model given in (2.3), equation (2.4) becomes

Pr(Y = 1|X, Z) = exp (α + βX + θ′Z) ∫_{−∞}^{∞} \frac{e^{γ_{X}u}}{1 + exp (α + βX + γ_{X}u + θ′Z)}dF (u|X, Z). (2.15)

In the presence of covariates Z, equation (2.15) generally does not reduce to a logistic model even if F (u|X, Z) has a simple form. In other words, the following logistic model usually fails,

Pr(Y = 1|X, Z) = \frac{exp (α^{∗} + β^{∗}X + θ^{∗′}Z)}{1 + exp (α^{∗} + β^{∗}X + θ^{∗′}Z)},

i.e.,

logit{Pr(Y = 1|X, Z)} = α^{∗} + β^{∗}X + θ^{∗′}Z.  (2.16)

Nonetheless, once γ_{X} and F (u|X, Z) are specified, the right side of equation (2.15) involves only the measured variables X and Z and the corresponding regression parameters α, β and θ. Thus, the maximum likelihood type method can be applied to model (2.15) to make inferences about β.

It is of practical importance to ascertain under which circumstances model (2.16) provides a reasonable approximation to model (2.15) and whether simple relationships such as (2.8) and (2.11) exist for the logistic regression. Straightforward algebraic manipulations show that model (2.15) is

logit{Pr(Y = 1|X, Z)} = α + βX + θ′Z + g(X, Z),  (2.17)

where

g(X, Z) = log \frac{e^{γ_{X}} \{1 + exp(α + βX + θ′Z)\} P_{X,Z} + \{1 + exp(α + βX + γ_{X} + θ′Z)\} (1 − P_{X,Z})}{\{1 + exp(α + βX + θ′Z)\} P_{X,Z} + \{1 + exp(α + βX + γ_{X} + θ′Z)\} (1 − P_{X,Z})}

for the binary confounder, and

g(X, Z) = γ_{X}μ_{X,Z}+0.5γ_{X}^{2} +log \frac{∫_{−∞}^{∞} \{1 + exp(α + γ_{X}^{2} + βX + γ_{X}u + θ′Z)\}^{−1} exp \{−\frac{(u−μ_{X,Z})^{2}}{2} \} du}{∫_{−∞}^{∞} \{1 + exp(α + βX + γ_{X}u + θ′Z)\}^{−1} exp \{−\frac{(u−μ_{X,Z})^{2}}{2} \} du}

for the normal confounder. If the event is rare or |γ_{X}| is small, g(X, Z) ≈ log {e^{γ_{X}} P_{X,Z} + (1 − P_{X,Z})} in the binary case and g(X, Z) ≈ γ_{X}μ_{X,Z} + 0.5γ_{X}^{2} in the normal case. It then follows that, under P_{X,Z} = P_{X} and μ_{X,Z} = μ_{X}, equation (2.17) becomes

logit\{Pr(Y = 1|X, Z)\} ≈ [α + log \{e^{γ_{0}}P_{0} + (1 − P_{0})\}] + \{β + log \frac{e^{γ_{1}}P_{1} + (1 − P_{1})}{e^{γ_{0}}P_{0} + (1 − P_{0})}\} X + θ′Z

(2.18)

and

logit{Pr(Y = 1|X, Z)} ≈ (α + γ_{0}μ_{0} + 0.5γ_{0}^{2})+{β + (γ_{1}μ_{1} − γ_{0}μ_{0}) + 0.5(γ_{1}^{2} − γ_{0}^{2})} X+θ′Z (2.19)

in the binary and normal cases, respectively. Models (2.18) and (2.19) take the form of (2.16) and imply that formulas (2.8)–(2.9) and (2.11)–(2..14) hold approximately for the logistic regression. Note that equations (2.13)–(2.14) hold approximately as long as the effects of X and Z on μ_{X,Z} are additive.

When model (2.3) and consequently model (2.15) hold whereas model (2.16) does not, the maximum likelihood estimator β̂^{∗} under model (2.16) converges to some well-defined limit, say β̃^{∗}, which minimizes the Kullback-Leibler distance between the true and assumed models (White, 1982). It is worthwhile to assess how close β̃^{∗} is to the value of β^{∗} determined by equation (2.8), (2.11) or (2.13), i.e., whether the substitution of β̂^{∗} under model (2.16) for β^{∗} in formula (2.8), (2.11) or (2.13) would lead to an approximately unbiased estimator of β. Unfortunately, there does not exist an explicit analytic expression for β̃^{∗}. We therefore resort to Monte Carlo simulation, as described below.

In our simulation studies, the binary responses were generated from model (2.3) with β = θ = 1, γ_{0} = γ_{1}, X taking the value 1 for 20% of the subjects and Z being a standard normal variable independent of X. The intercept term α was set at different values to yield certain overall event rates. There were 1000 simulation samples, each consisting of 1000 subjects. For each simulation sample, we calculated the maximum likelihood estimate β̂^{∗} under model (2.16) and replaced β^{∗} in formula (2.8) or (2.13) by β̂^{∗} to obtain β̂ corresponding to specific values of the sensitivity parameters. In the case of binary U, we first set P_{0} = 0.2 and P_{1} = 0.9. Under Γ = (2, 3, 4), the Monte Carlo estimates for the sampling means of β̂ were approximately (1.01, 1.02, 1.03) when the overall event rate was 10% and approximately (1.02, 1.05, 1.08) when the overall event rate was 50%. For P_{0} = 0.4 and P_{1} = 0.6, the sampling means of β̂ under Γ = (2, 3, 4) were estimated at (1.00, 0.98, 0.96) and (0.98, 0.96, 0.94) when the overall event rates were 10% and 50%, respectively. Clearly, all these sampling means are very close to 1, the true value of β. In the case of normal U, we set the overall event rate at 10%. When μ_{0} = 0 and μ_{1} = 1, the Monte Carlo estimates for the sampling means of β̂ were approximately (0.99, 0.96, 0.90 0.80) under γ = (0.25, 0.5, 0.75, 1). The corresponding estimates were approximately (0.99, 0.94, 0.86, 0.75) when μ_{0} = 0 and μ_{1} = 2. The bias of β̂ became smaller as the event rate decreased.

These numerical results suggest that, for binary U, one may obtain β̂^{∗} and its variance estimate from an existing logistic regression software and then use formulas (2.8)–(2.9) to perform the sensitivity analysis, while for normal U, formulas (2.11)–(2.14) are good approximations only when the overall event rate is below 10% or |γ_{X}| < 0.75, i.e., 0..5 < Γ_{X} < 2.

2.3. An Alternative Parameterization for the Association Between U and X

In Sections 2.1 and 2.2, we parameterized the association between the unmeasured confounder U and the exposure indicator X in terms of the conditional distribution of U given X. An alternative parameterization, as used by Rosenbaum and Rubin (1983a) and Rosenbaum (1995, Ch. 4), specifies the effect of U on X. There is a simple relationship between these two types of parameterizations when U is binary. Specifically, let ϕ and P be the overall prevalences of the exposure and unmeasured confounder, respectively, and let Ψ be the relative risk of being exposed associated with the unmeasured confounder, i.e., ϕ = Pr(X = 1), P = Pr(U = 1), and Ψ = Pr(X = 1|U = 1)/Pr(X = 1|U = 0). Note that Ψ measures directly the degree of departure from randomization due to the unmeasured confounder. By Bayes’ theorem,

Ψ = \frac{P_{1}/P}{(1 − P_{1})/(1 − P)},

which implies that

P_{1} = \frac{ΨP}{1 − P + ΨP} .  (2.20)

On the other hand, because P = P_{1}ϕ + P_{0}(1 − ϕ), we have

P_{0} = \frac{P − ϕP_{1}}{1 − ϕ} .  (2.21)

By replacing P_{1} and P_{0} in (2.9) by (2.20) and (2.21) and also replacing P_{1} in (2.21) by expression (2.20), we obtain an expression that relates R^{∗} to R through Γ_{0}, Γ_{1}, Ψ, P and ϕ.

We may also measure the effect of U on X by the odds ratio, say Q. Clearly,

Q = \frac{P_{1}(1 − P_{0})}{P_{0}(1 − P_{1})}.  (2.22)

Replacing P_{0} in (2.22) by expression (2.21), we can express P_{1} as a function of P, ϕ and Q by solving a quadratic equation. The resulting expression for P_{1}, together with (2.21), enables us to relate R^{∗} and R through Γ_{0}, Γ_{1}, Q, P and ϕ. This formulation is essentially the same as that of Rosenbaum and Rubin (1983a). In view of this connection, the proposed methodology may be regarded as a generalization of the pioneering work by Rosenbaum and Rubin. It seems preferable to use the relative risk formulation rather than the odds ratio formulation for the effect of U on X not only because the former leads to simpler formulas but also because the relative risk is more relevant than the odds ratio especially when the exposure is not rare enough to ensure adequate approximation of the relative risk by the odds ratio.

The parameter ϕ may be estimated from the current study or from an external data source, the resulting estimator being denoted by ϕ̂. With the restriction that Γ_{0} = Γ_{1} = Ψ or Γ_{0} = Γ_{1} = Q, the results of the sensitivity analysis can be summarized in a concise manner, especially if the results are insensitive to the variation of P over a broad range of plausible values. The estimation of the variance of β̂ or R̂ is more complicated under the specification of Ψ or Q than under the direct specification of (P_{0}, P_{1}) if one wishes to adjust for the extra random variation due to the estimation of ϕ. Nonetheless, it is straightforward to establish the asymptotic joint distribution of (β̂^{∗}, ϕ̂) by the multivariate central limit theorem and the law of large numbers; the asymptotic normality and variance of β̂ then follow from the delta-method. The variation of ϕ̂ is negligible for large studies.

3. SURVIVAL TIME

Let T denote the survival time or failure time of interest, and let X, Z and U be as defined at the start of Section 2. The proportional hazards regression (Cox, 1972) specifies that the hazard functions of T conditional on the sets of covariates (X, Z, U) and (X, Z) are, respectively,

λ(t|X, Z, U) = λ_{0}(t) exp (βX + γ_{X}U + θ′Z) ,  (3.1)

and

λ(t|X, Z) = λ_{0}^{∗}(t) exp (β^{∗}X + θ^{∗′}Z) ,  (3.2)

where λ_{0}(·) and λ_{0}^{∗}(·) are arbitrary baseline hazard functions, and (β, γ_{0}, γ_{1}, θ) and (β^{∗}, θ^{∗}) are unknown regression parameters. Models (3.1) and (3.2) parallel (2.1) and (2.2) as well as (2.3) and (2.16). As in Section 2, we wish to ascertain the relationship between β^{∗} and β.

Again, let F (u|X, Z) be the distribution function of U given X and Z. Also, let f(t|·) and S(t|·) denote the conditional density and survival functions of T. By elementary probability arguments,

λ(t|X, Z) = \frac{f(t|X, Z)}{S(t|X, Z)} = \frac{∫_{−∞}^{∞} f(t|X, Z, u)dF (u|X, Z)}{∫_{−∞}^{∞} S(t|X, Z, u)dF (u|X, Z)}.  (3.3)

Under model (3.1),

∫_{−∞}^{∞} f(t|X, Z, u)dF (u|X, Z) = ∫_{−∞}^{∞} λ_{0}(t)e^{βX+γ_{X}u+θ′Z} exp {−Λ_{0}(t)e^{βX+γ_{X}u+θ′Z}} dF (u|X, Z),

∫_{−∞}^{∞} S(t|X, Z, u)dF (u|X, Z) = ∫_{−∞}^{∞} exp {−Λ_{0}(t)e^{βX+γ_{X}u+θ′Z}} dF (u|X, Z),

where Λ_{0}(t) = ∫_{0}^{t} λ_{0}(s)ds. Thus, (3.3) becomes

λ(t|X, Z) = λ_{0}(t) exp (βX + θ′Z) h(t; X, Z),  (3.4)

where

h(t; X, Z) = \frac{∫_{−∞}^{∞} e^{γ_{X}u} exp \{−Λ_{0}(t)e^{βX+γ_{X}u+θ′Z}\} dF (u|X, Z)}{∫_{−∞}^{∞} exp \{−Λ_{0}(t)e^{βX+γ_{X}u+θ′Z}\} dF (u|X, Z)}.

We will evaluate the right side of (3.4) for the binary and normal U.

Binary confounder. If F (u|X, Z) is a Bernoulli distribution with success probability P_{X,Z}, then

h(t; X, Z) = \frac{e^{γ_{X}} exp \{−Λ_{0}(t)e^{βX+γ_{X}+θ′Z}\} P_{X,Z} + exp \{−Λ_{0}(t)e^{βX+θ′Z}\} (1 − P_{X,Z})}{exp \{−Λ_{0}(t)e^{βX+γ_{X}+θ′Z}\} P_{X,Z} + exp \{−Λ_{0}(t)e^{βX+θ′Z}\} (1 − P_{X,Z})}.  (3.5)

To make the right side of (3.5) more tractable, we assume that the event is rare or the effect of U on T is small, which implies that Λ_{0}(t) or |γ_{X}| is small. Then (3.5) simplifies to

h(t; X, Z) ≈ e^{γ_{X}} P_{X,Z} + (1 − P_{X,Z}),

which entails that (3.4) is

λ(t|X, Z) ≈ λ_{0}(t) exp (βX + θ′Z) {e^{γ_{X}} P_{X,Z} + (1 − P_{X,Z})} .  (3.6)

If U is conditionally independent of Z given X, then (3.6) is equivalent to

λ(t|X, Z) ≈ λ_{0}(t) \{e^{γ_{0}}P_{0} + (1 − P_{0})\} exp [\{β + log \frac{e^{γ_{1}}P_{1} + (1 − P_{1})}{e^{γ_{0}}P_{0} + (1 − P_{0})}\} X + θ′Z] .  (3.7)

In view of (3.2) and (3.7),

β ≈ β^{∗} − log \frac{e^{γ_{1}}P_{1} + (1 − P_{1})}{e^{γ_{0}}P_{0} + (1 − P_{0})},  (3.8)

which takes the same form as formula (2.8).

Normal confounder. For the normally distributed U, it can be shown that

h(t; X, Z) = exp (γ_{X}μ_{X,Z} + 0.5γ_{X}^{2}) \frac{∫_{−∞}^{∞} exp \{−Λ_{0}(t)e^{γ_{X}^{2}+βX+γ_{X}u+θ′Z}\} exp \{−\frac{(u−μ_{X,Z})^{2}}{2} \} du}{∫_{−∞}^{∞} exp \{−Λ_{0}(t)e^{βX+γ_{X}u+θ′Z}\} exp \{−\frac{(u−μ_{X,Z})^{2}}{2} \} du},

(3.9)

which can be approximated by exp(γ_{X}μ_{X,Z} + 0.5γ_{X}^{2}) if Λ_{0}(t) or |γ_{X}| is small.. Then

λ(t|X, Z) ≈ λ_{0}(t) exp (βX + θ′Z + γ_{X}μ_{X,Z} + 0.5γ_{X}^{2}) ,  (3.10)

which is similar to (2.10). Under μ_{X,Z} = μ_{X}, equation (3.10) entails that

β ≈ β^{∗} − {(γ_{1}μ_{1} − γ_{0}μ_{0}) + 0.5(γ_{1}^{2} − γ_{0}^{2})} ,  (3.11)

If γ_{0} = γ_{1} = γ and μ_{X,Z} = μ_{X} + q(Z), then

β ≈ β^{∗} − γδ,  (3.12)

where again δ = μ_{1} − μ_{0}. Formulas (3.11) and (3.12) parallel (2.11) and (2.13).

The above results indicate that, at least for rare events, the sensitivity analysis for survival time data can be performed in the same manner as that of binary outcome discussed in Section 2. Given the sensitivity parameters (γ_{X}, P_{X,Z}) and (γ_{X}, μ_{X,Z}), the partial likelihood theory can be applied to models (3.6) and (3.10), respectively, to make inferences about β. If U is conditionally independent of Z given X, one may estimate β^{∗} by fitting model (3.2) with any existing Cox regression software and then use formula (3.8) and (3.11) to conduct the sensitivity analysis. Likewise, if γ_{0} = γ_{1} and μ_{X,Z} = μ_{X} + q(Z), then formula (3.12) may be used.

As in the setting of logistic regression, it is important to investigate how accurate the approximations (3.8), (3.11) and (3.12) are. In general, model (3.2) does not hold exactly because, as evident from (3.5) and (3.9), the factor h(t; X, Z) in (3.4) is a complicated function of t that depends on X and Z, which implies that the right side of (3.4) does not reduce to that of (3.2). In particular, the hazard functions for the exposed and unexposed are not proportional over time. If one fits model (3.2) to the observable data, then the maximum partial likelihood estimator of β^{∗}, denoted by β̂^{∗}, will converge to some well-defined limit, say β̃^{∗}, which represents an average over time of the log hazard ratios between the exposed and unexposed (Struthers and Kalbfleisch, 1986; Lin and Wei, 1989). The question then arises as to how different β̃^{∗} is from β^{∗} used in formulas (3.8), (3.11) and (3.12).

To answer the above question, we conducted a series of Monte Carlo studies similar to that of Section 2.2. The survival times were generated from model (3.1) with λ_{0}(·) = 1, β = θ = 1, γ_{0} = γ_{1}, X taking the value 1 for 20% of the subjects and Z being a standard normal variable independent of X. The censoring times were generated from the uniform(0, τ) distribution, where τ was chosen such that a desired percentage of the survival times were censored. Again, there were 1000 simulation samples, each with 1000 subjects. For each sample, we calculated the maximum partial likelihood estimate β̂^{∗} under model (3.2) and replaced β^{∗} in formula (3.8) or (3.12) by β̂^{∗} to yield β̂. In the case of binary U, we set P_{0} = 0.2 and P_{1} = 0.9. Under Γ = (2, 3, 4), the Monte Carlo estimates for the sampling means of β̂ were approximately (1.01, 1.01, 1.02) with 90% censoring and (1.00, 0.98, 0.97) with 50% censoring. Similar results were obtained for other values of P_{0} and P_{1}. In the case of normal U, we set μ_{0} = 0 and μ_{1} = 1 and created 90% censoring. The sampling means of β̂ were found to be approximately (0.99, 0.95, 0.88, 0.79) under γ = (0.25, 0.5, 0.75, 1). The bias of β̂ increased slightly when μ_{1} was changed from 1 to 2 and decreased as censoring became heavier.

These numerical results suggest that, for binary U, the maximum partial likelihood estimator β̂^{∗} under model (3.2) may be used in formula (3.8) to yield an approximately unbiased estimator of β given the sensitivity parameters, while for normal U, formulas (3.11) and (3.12) are good approximations only when the censoring percentage is greater than 90% or |γ_{X}| < 0.75, i.e., 0.5 < Γ_{X} < 2. Note that the variance of β̂^{∗} may be estimated by the partial likelihood method, or more appropriately by the sandwich estimator of Lin and Wei (1989), which is available in most software packages. For correlated survival times, the robust variance estimator can be found in Lin (1994).

4. REAL EXAMPLES

We now apply the methods developed in Sections 2 and 3 to two recently published medical studies. The first one is a case-control study investigating the potential role of appetite-suppressant drugs in the development of primary pulmonary hypertension (Abenhaim et al., 1996), and the second one is a prospective cohort study examining the association between the use of right heart catheterization during the first 24 hours of care in the intensive care unit and subsequent survival (Connors et al., 1996). Both studies aroused considerable controversy because they were observational studies that raised concerns about the safety of commonly used interventions.

4.1. Primary Pulmonary Hypertension Study

Primary pulmonary hypertension (PPH) is a rare, often fatal disease with unknown causes. In the early 1990s, PPH developed in a cluster of French patients exposed to derivatives of fenfluramine in appetite suppressants used for weight control. This phenomenon prompted a major European case-control study which evaluated 95 patients with PPH and 355 controls matched individually on gender and age. The investigators used conditional logistic regression to study the association between the exposure to appetite suppressants and the development of PPH after adjustment for weight-related confounding variables and other suspected risk factors. With the adjustment, the odds ratio of PPH among those who had taken appetite suppressants for a total of more than 3 months was estimated at 23.1, with a 95% confidence interval of (6.9, 77.7).

Because the factors leading to the development of PPH remain enigmatic, it is worthwhile to assess the sensitivity of the findings of this case-control study to the assumptions about unmeasured confounders. To this end, we display in Table 2 the point and interval estimates for the odds ratio of interest after adjustment for an unmeasured binary confounder according to formula (2.9) with Γ_{0} = Γ_{1} = Γ. Here, P_{1} and P_{0} are, respectively, the prevalences of the unmeasured confounder among those who had used appetite suppressants for more than 3 months and among those who had never used them, and Γ is the odds ratio of PPH associated with the unmeasured confounder. We only display the estimates for P_{1} ≥ P_{0}. The estimates will be inflated upward if P_{1} < P_{0}. In Table 2, the point estimates are always much higher than 1 while the lower limits of the 95% confidence intervals are greater than 1 except for two entries under Γ = 10. Recall that the adjustment factor A is bounded above by Γ. Thus, regardless of the values of P_{0} and P_{1}, the point estimate will remain above 1 as long as Γ is less than 23.1, and the lower confidence limit will be greater than 1 if Γ is below 6.9. When Γ is larger than 6.9, the lower confidence limit may drop below 1 if P_{0} is close to 0 and P_{1} is close to 1.

It is also informative to assess how the adjustment for an unmeasured normal confounder might affect the estimates. Formula (2.14) implies that, under δ = 1, the point or interval estimate adjusting for the normal confounder will remain above 1 as long as Γ is less than the unadjusted estimate. Under δ = 2, the point estimate and the lower limit of the 95% confidence interval will remain above 1 if Γ is less than 4.81 and 2.63, respectively.

The results of this sensitivity analysis strongly support the investigators’ conclusion that the long-term use of appetite-suppressant drugs increased the risk of PPH. It should be pointed out that sensitivity analyses do not and can not account for potential detection and recall biases, which the investigators judged to be minimal in this study.

4.2. Right Heart Catheterization Study

Many cardiologists and intensive-care physicians believe that the direct measurement of cardiac function provided by right heart catheterization (RHC) is necessary to guide therapy for certain critically ill patients and that such management leads to better outcomes. This popular belief is so strong that it has prevented the conduct of a randomized clinical trial even though the benefit of RHC has never been demonstrated in a clinical trial. Connors et al. (1996) reported a prospective cohort study designed to assess the effectiveness of RHC in the initial care of critically ill patients. The study involved 5,735 patients receiving care in the intensive care unit for one of 9 prespecified disease categories. The investigators did an exceptionally thorough analysis of their data. Not only did they attempt to identify all the major risk factors, but they also constructed a propensity score (Rosenbaum and Rubin, 1983b) for RHC to adjust for treatment selection bias. Using a proportional hazards model that included the propensity score as well as 9 risk factors, they estimated the hazard ratio or relative hazard of death for patients managed with RHC at 1.21 with a 95% confidence interval of (1.09, 1.25). It is interesting to assess the sensitivity of this unexpected finding to potential unknown confounders.

Table 3 shows the point and interval estimates for the relative hazard of death associated with RHC after adjustment for an unmeasured binary confounder according to formula (3.8) with γ_{0} = γ_{1} ≈ 1.1, i.e., Γ_{0} = Γ_{1} = Γ = 3. Here, P_{1} and P_{0} are, respectively, the prevalences of the unmeasured confounder among those who were managed with RHC and among those who were not, and Γ is the relative hazard of death for the unmeasured confounder. As shown in the table, the point estimates for the relative hazard of death are no larger than 1.01 if P_{1} −P_{0} ≥ 0.2, and the lower limits of the 95% confidence intervals drop below 1 if P_{1} − P_{0} ≥ 0.1. Furthermore, the point estimates are no larger than 0.81 if (P_{0}, P_{1}) = (0, 0.3), (0.1, 0.4), (0.2, 0.6), (0.3, 0.7), (0.4, 0.9) or (0.5, 1.0), and the upper confidence limits are about 0.8 or smaller if P_{1} − P_{0} > 0.5.

Table 4 displays the estimates for the relative hazard of death associated with RHC under the alternative parameterization described in Section 2.3. Here, Ψ is the relative risk of RHC use associated with the unmeasured confounder, and Γ is again the relative hazard of death for the unmeasured confounder. Under Ψ = Γ = 2, the point estimates for the relative hazard of death are very close to the null value of 1 if 0.2 ≤ P ≤ 0.5; under Ψ = Γ = 3, the estimates lie between 0.73 and 0.82 if 0.1 ≤ P ≤ 0.5. In other words, over a broad range of plausible values for the prevalence of the unmeasured confounder, ranging from 0.1 to 0.5, a true relative hazard of 0.8 could be misrepresented as 1.21 if an unmeasured confounder increased both the hazard of death and probability of RHC 3-fold, and a true relative hazard of 1.0 could be misrepresented as 1.21 if an unmeasured confounder increased both the hazard of death and probability of RHC 2-fold.

Connors et al. (1995) also performed a sensitivity analysis. They stated that “the covariate missing from the propensity score would have to increase the risk of death 6-fold and increase the probability of RHC 6-fold for a true relative hazard of 0.80 to be misrepresented as a relative hazard of 1.21. In addition, a true relative hazard of 1.0 could be misrepresented as 1.21 if a missing covariate that simultaneously increased the risk of death 3-fold and the probability of RHC 3-fold was not accounted for in the propensity score.” Apparently, these results differ from those just described in the preceding paragraph. There are several reasons for these discrepancies. First of all, due to the lack of sensitivity analysis methods for censored survival time data, Connors et al. obtained their results by applying the method of Rosenbaum and Rubin (1983a) to the binary outcome of 30-day mortality rather than patient survival time. The point estimate for the odds ratio of death as actually used in their sensitivity analysis was 1.29, which is slightly higher than the reported relative hazard estimate of 1.21. More important, the method of Rosenbaum and Rubin formulates the effect of the unmeasured confounder on treatment selection, which in this case is the use of RHC, in terms of the odds ratio, but Connors et al. described this odds ratio as the relative risk or probability ratio in their statement. Because RHC is a common procedure, an increase in the probability of RHC translates into a much larger increase in the odds of RHC. Thus, the misreporting of the odds ratio of RHC as the probability ratio of RHC exaggerated the effects of an unmeasured confounder that would be required to misrepresent a neutral or beneficial effect of RHC as harmful. Dr. Connors confirmed that “risk of death”, “probability of RHC” and “relative hazard of 1.21” in their statement should read “odds of death”, “odds of RHC” and “odds ratio of 1.29” and has submitted a correction letter to Journal of the American Medical Association. Given this clarification, we performed an additional sensitivity analysis on the odds ratio estimate of 1.29 using the sensitivity parameters that were actually used by Conners et al. For P = 0.5, which is the overall prevalence implicitly assumed by Connors et al., we estimated the true odds ratios of death associated with RHC at 0.99 and 0.73 under Q = Γ = 3 and Q = Γ = 6, respectively, which are close to the values of 1 and 0.8 obtained by Connors et al. Here, Q is the odds ratio of RHC associated with the unmeasured confounder.

It should be pointed out that, because of the extraordinary efforts made by the investigators to identify and adjust for all the known risk factors, it is unlikely that there exists an unmeasured confounder with Q = Γ = 3 or 6 even though these were the values considered by the investigators in their sensitivity analysis. On the other hand, the lower limit of the observed 95% confidence interval was 1.09, which is not far above the null value of 1.0. Thus, the adjustment for an unmeasured confounder that has only small to modest effects on the use of RHC and death would drop the lower confidence limit below 1 or even 0.8.

The results of our sensitivity analysis provide additional insights into this important study, and imply perhaps greater uncertainty about the effect of RHC than those stated in the original report. The observed relative hazard of 1.21 with a 95% confidence interval of (1.09, 1.25) suggests a fairly small increase in the hazard of death due to RHC, and could represent a true relative hazard of 1 or 0.8 if there existed an unmeasured confounder of moderate independent effects. Thus, this study did not provide strong evidence for either a harmful or beneficial effect of RHC. It would be desirable to carry out a randomized controlled clinical trial, as recommended by Connors et al.

Incidentally, the censoring percentage in this study was approximately 65%. According to the simulation results reported at the end of Section 3, formula (3.8) is a good approximation under such a censoring level whereas formulas (3.11)–(3.12) may not be. Thus, we did not perform a sensitivity analysis with an unmeasured normal confounder for this study.

5. DISCUSSION AND FURTHER RESULTS

Observational studies provide a valuable source of information for investigating the effects of treatments, exposures, etc., though these studies are subject to potential biases induced by unmeasured or unknown confounders. The methods of sensitivity analysis as presented in this paper enable one to assess quantitatively how the findings of a study might be altered by unknown confounders of various magnitudes. These methods, however, cannot be used to test whether or not unknown confounders exist, nor can they be used to assess or eliminate the potential effects of other forms of bias such as selection bias.

If Γ_{0} = Γ_{1}, then formula (2.9) takes the same form as formula (1) of Schlesselman (1978), which becomes formula (1) of Cornfield et al. (1959) under R = 1. Thus, Section 2 establishes that the basic relationships between the true and apparent associations derived by Cornfield and Schlesselman for the simple 2×2 table with an unmeasured binary confounder continue to hold when the log-linear/logistic regression models are used to adjust for the measured confounders that are conditionally independent of the unmeasured binary confounder. Section 3 further demonstrates that similar relationships hold for the Cox regression with censored survival data. The development of simple and explicit formulas, such as (2.8)–(2.9), (2.11)–(2.14), (3.8) and (3.11)–(3.12), will facilitate the use of the powerful tool of sensitivity analysis by practitioners.

Compared with the published methods, the proposed technique has the advantage of accommodating both discrete and continuous covariates. In observational studies, it is desirable and customary to collect information on all previously identified confounders and to adjust for them in the estimation of the exposure effect. The method of sensitivity analysis presented in this paper enables one to assess the residual confounding effects due to unmeasured or unknown covariates after adjusting for the effects of the measured covariates. Without the adjustment for the known and measured confounders, the sensitivity analysis would need to encompass much broader ranges of values for the sensitivity parameters, and therefore become less informative.

One complication in assessing the sensitivity of the treatment estimate which adjusts for covariates through nonlinear regression is that the reduced model omitting the unmeasured confounder does not hold exactly. However, the analytical and numerical investigations in Sections 2 and 3 demonstrated that the reduced logistic model (2.16) and proportional hazards model (3.2) provide reasonable approximations for rare events. The requirement of rare events does not seem to be a major limitation of the proposed methodology because the rarity of the event is often one of the main reasons for conducting observational studies, especially case-control studies. Furthermore, the numerical studies in Sections 2 and 3 indicated that, for unmeasured binary confounders, the approximations are adequate even if the event is not rare.

Sections 2 and 3 considered regression models which are sufficiently general to be useful in a variety of applications and which at the same time are simple enough to permit straightforward analysis. The ideas presented in these sections can be applied to other, possibly more complicated, models. In particular, it is straightforward to show that formulas (2.8)–(2.9), (2.11)–(2.14), (3.8) and (3.11)–(3.12) hold when the corresponding regression models are further stratified on some categorical variables provided that β, θ and the sensitivity parameters remain the same among the strata.

In this paper, we focused on binary response and censored survival data because they are the two most common endpoints in observational studies, especially in biomedical applications. The approach taken here may also be applied to other outcome measures. It is worthwhile to discuss briefly the linear regression for continuous response variables. Let Y be a continuous response variable, and let X, Z and U be as defined in Sections 2 and 3. Suppose that Y is related to X, Z and U through the following linear model

Y = βX + γ_{X}U + θ′Z + ϵ,

and the fitted model is

Y = β^{∗}X + θ^{∗′}Z + ϵ^{∗},

where (β, γ_{0}, γ_{1}, θ) and (β^{∗}, θ^{∗}) are unknown regression parameters, and ϵ and ϵ^{∗} are error terms with unspecified distributions. For survival data, it is common to set Y = log T, which corresponds to the accelerated failure time model (Kalbfleisch and Prentice, 1980, pp. 32–34). It is easy to show that, under the conditional independence of U and Z given X,

β = β^{∗} − (γ_{1}m_{1} − γ_{0}m_{0}),  (5.1)

where m_{1} and m_{0} are the means of U for the exposed and unexposed, respectively. If γ_{0} = γ_{1} = γ, then equation (5.1) reduces to

β = β^{∗} − γd,  (5.2)

where d = m_{1} −m_{0}. Note that formula (5.2) holds as long as the effects of X and Z on the mean of U are additive. Incidentally, Rosenbaum (1986) used a formula similar to (5.2) to study the effect of dropping out of high school on cognitive achievement test scores. Of course, for noncensored data, β^{∗} can be estimated by the familiar least-squared method. The estimation of β^{∗} for censored survival data has been studied by Tsiatis (1990), Wei, Ying and Lin (1990) and Lai and Ying (1991).

The new method as well as all the published methods only deal with a single unmeasured confounder. Nevertheless, the approach taken in Sections 2 and 3 can be generalized to the setting of multiple unmeasured confounders. One would then need to specify the joint distribution of the confounders as well as their joint effects on the outcome variable. A more practical alternative is to simply treat U as a composite score constructed from a set of unmeasured confounders. In this regard, a normal U may represent a linear combination of several confounders, while a binary U may represent the dichotomy of “high risk” versus “low risk” determined by multiple risk factors. Note that the use of a single unmeasured confounder may not be a severe restriction if all the known confounders are adjusted for through regression modelling.

Finally, this paper complements the current literature on misspecified nonlinear regression models. Gail, Wieand and Piantadosi (1984), Lagakos and Schoenfeld (1984) and Struthers and Kalbfleisch (1986) investigated the bias for the estimator of treatment effect in randomized experiments when relevant covariates are omitted from logistic and proportional hazards models. Our paper extends their results to non-randomized studies, in which the omitted covariates are not independent of treatment.

ACKNOWLEDGEMENTS

The authors are grateful to Frank Harrell for helpful discussions. This research was supported by the U.S. National Institutes of Health Grants GM47845, CA39929, AI29168, HL40628, HL43201 and AG09556, and Contract HC85079.

REFERENCES

Abenhaim, L., Moride, Y., Brenot, F., et al. (1996). Appetite-suppressant drugs and the risk of primary pulmonary hypertension. The New England Journal of Medicine, 335, 609–616.

Bross, I. D. J. (1966). Spurious effects from an extraneous variable. Journal of Chronic Diseases, 19, 637–647.

Connors, A. F., Speroff, T., Dawson, N. V., et al. (1996). The effectiveness of right heart catheterization in the initial care of critically ill patients. Journal of the American Medical Association, 276, 889–897.

Cornfield, J., Haenszel, W., Hammond, E. C., Lilienfeld, A. M., Shimkin, M. B. and Wynder, E. L. (1959). Smoking and lung cancer: Recent evidence and a discussion of some questions. Journal of the National Cancer Institute, 22, 173-203.

Cox, D. R. (1972). Regression models and life-tables (with discussion). Journal of the Royal Statistical Society B, 34, 187–220.

Fisher, R. A. (1958). Lung cancer and cigarettes? Nature, 182, 108.

Gail, M. H., Wieand, S. and Piantadosi, S. (1984). Biased estimates of treatment effect in randomized experiments with nonlinear regressions and omitted covariates. Biometrika, 71, 431–444.

Kalbfleisch, J. D. and Prentice, R. L. (1980). The Statistical Analysis of Failure Time Data. New York: Wiley.

Lagakos, S. W. and Schoenfeld, D. A. (1984). Properties of proportional-hazards score tests under misspecified regression models. Biometrics, 40, 1037–1048.

Lai, T. L. and Ying, Z. (1991). Large sample theory of a modified Buckley-James estimator for regression analysis with censored data. The Annals of Statistics, 19, 1370–1402.

Lin, D. Y. (1994). Cox regression analysis of multivariate failure time data: The marginal approach. Statistics in Medicine, 13, 2233–2247.

Lin, D. Y. and Wei, L. J. (1989). The robust inference for the Cox proportional hazards model. Journal of the American Statistical Association, 84, 1074–1078.

Rosenbaum, P. R. (1986). Dropping out of high school in the United States: An observational study. Journal of Educational Statistics, 11, 207–224.

Rosenbaum, P. R. (1995). Observational Studies. New York: Springer-Verlag.

Rosenbaum, P. R. and Rubin, D. B. (1983a). Assessing sensitivity to an unobserved binary covariate in an observational study with binary outcome. Journal of the Royal Statistical Society, B 11, 212–218.

Rosenbaum, P. R. and Rubin, D. B. (1983b). The central role of the propensity score in observational studies for causal effects. Biometrika, 70, 41–55.

Schlesselman, J. J. (1978). Assessing effects of confounding variables. American Journal of Epidemiology, 108, 3–8.

Struthers, C. A. and Kalbfleisch, J. D. (1986). Misspecified proportional hazard models. Biometrika, 73, 363–369.

Tsiatis, A. A. (1990), Estimating regression parameters using linear rank tests for censored data. The Annals of Statistics, 18, 354–372.

Wei, L. J., Ying, Z. and Lin, D. Y. (1990). Linear regression analysis of censored survival data based on rank tests. Biometrika, 77, 845–851.

White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica, 50, 1–25.

Table 1. The Adjustment Factor A = {ΓP_{1} + (1 − P_{1})}/{ΓP_{0} + (1 − P_{0})} for a Unmeasured

Binary Confounder

P_{0} = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Γ = 2 P_{1} = 0.0  1.00 0.91 0.83 0.77 0.71 0.67 0.62 0.59 0.56 0.53 0.50
P_{1} = 0.1  1.10 1.00 0.92 0.85 0.79 0.73 0.69 0.65 0.61 0.58 0.55
P_{1} = 0.2  1.20 1.09 1.00 0.92 0.86 0.80 0.75 0.71 0.67 0.63 0.60
P_{1} = 0.3  1.30 1.18 1.08 1.00 0.93 0.87 0.81 0.76 0.72 0.68 0.65
P_{1} = 0.4  1.40 1.27 1.17 1.08 1.00 0.93 0.87 0.82 0.78 0.74 0.70
P_{1} = 0.5  1.50 1.36 1.25 1.15 1.07 1.00 0.94 0.88 0.83 0.79 0.75
P_{1} = 0.6  1.60 1.45 1.33 1.23 1.14 1.07 1.00 0.94 0.89 0.84 0.80
P_{1} = 0.7  1.70 1.55 1.42 1.31 1.21 1.13 1.06 1.00 0.94 0.89 0.85
P_{1} = 0.8  1.80 1.64 1.50 1.38 1.29 1.20 1.12 1.06 1.00 0.95 0.90
P_{1} = 0.9  1.90 1.73 1.58 1.46 1.36 1.27 1.19 1.12 1.06 1.00 0.95
P_{1} = 1.0  2.00 1.82 1.67 1.54 1.43 1.33 1.25 1.18 1.11 1.05 1.00
Γ = 3 P_{1} = 0.0  1.00 0.83 0.71 0.62 0.56 0.50 0.45 0.42 0.38 0.36 0.33
P_{1} = 0.1  1.20 1.00 0.86 0.75 0.67 0.60 0.55 0.50 0.46 0.43 0.40
P_{1} = 0.2  1.40 1.17 1.00 0.87 0.78 0.70 0.64 0.58 0.54 0.50 0.47
P_{1} = 0.3  1.60 1.33 1.14 1.00 0.89 0.80 0.73 0.67 0.62 0.57 0.53
P_{1} = 0.4  1.80 1.50 1.29 1.12 1.00 0.90 0.82 0.75 0.69 0.64 0.60
P_{1} = 0.5  2.00 1.67 1.43 1.25 1.11 1.00 0.91 0.83 0.77 0.71 0.67
P_{1} = 0.6  2.20 1.83 1.57 1.38 1.22 1.10 1.00 0.92 0.85 0.79 0.73
P_{1} = 0.7  2.40 2.00 1.71 1.50 1.33 1.20 1.09 1.00 0.92 0.86 0.80
P_{1} = 0.8  2.60 2.17 1.86 1.62 1.44 1.30 1.18 1.08 1.00 0.93 0.87
P_{1} = 0.9  2.80 2.33 2.00 1.75 1.56 1.40 1.27 1.17 1.08 1.00 0.93
P_{1} = 1.0  3.00 2.50 2.14 1.88 1.67 1.50 1.36 1.25 1.15 1.07 1.00
Γ = 4 P_{1} = 0.0  1.00 0.77 0.62 0.53 0.45 0.40 0.36 0.32 0.29 0.27 0.25
P_{1} = 0.1  1.30 1.00 0.81 0.68 0.59 0.52 0.46 0.42 0.38 0.35 0.33
P_{1} = 0.2  1.60 1.23 1.00 0.84 0.73 0.64 0.57 0.52 0.47 0.43 0.40
P_{1} = 0.3  1.90 1.46 1.19 1.00 0.86 0.76 0.68 0.61 0.56 0.51 0.47
P_{1} = 0.4  2.20 1.69 1.38 1.16 1.00 0.88 0.79 0.71 0.65 0.59 0.55
P_{1} = 0.5  2.50 1.92 1.56 1.32 1.14 1.00 0.89 0.81 0.74 0.68 0.62
P_{1} = 0.6  2.80 2.15 1.75 1.47 1.27 1.12 1.00 0.90 0.82 0.76 0.70
P_{1} = 0.7  3.10 2.38 1.94 1.63 1.41 1.24 1.11 1.00 0.91 0.84 0.77
P_{1} = 0.8  3.40 2.62 2.12 1.79 1.55 1.36 1.21 1.10 1.00 0.92 0.85
P_{1} = 0.9  3.70 2.85 2.31 1.95 1.68 1.48 1.32 1.19 1.09 1.00 0.93
P_{1} = 1.0  4.00 3.08 2.50 2.11 1.82 1.60 1.43 1.29 1.18 1.08 1.00

Table 1. The Adjustment Factor A = {ΓP_{1} + (1 − P_{1})}/{ΓP_{0} + (1 − P_{0})} for a Unmeasured

Binary Confounder (Continued)

P_{0} = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Γ = 6  P_{1} = 0.0  1.00 0.67 0.50 0.40 0.33 0.29 0.25 0.22 0.20 0.18 0.17
P_{1} = 0.1  1.50 1.00 0.75 0.60 0.50 0.43 0.38 0.33 0.30 0.27 0.25
P_{1} = 0.2  2.00 1.33 1.00 0.80 0.67 0.57 0.50 0.44 0.40 0.36 0.33
P_{1} = 0.3  2.50 1.67 1.25 1.00 0.83 0.71 0.62 0.56 0.50 0.45 0.42
P_{1} = 0.4  3.00 2.00 1.50 1.20 1.00 0.86 0.75 0.67 0.60 0.55 0.50
P_{1} = 0.5  3.50 2.33 1.75 1.40 1.17 1.00 0.88 0.78 0.70 0.64 0.58
P_{1} = 0.6  4.00 2.67 2.00 1.60 1.33 1.14 1.00 0.89 0.80 0.73 0.67
P_{1} = 0.7  4.50 3.00 2.25 1.80 1.50 1.29 1.12 1.00 0.90 0.82 0.75
P_{1} = 0.8  5.00 3.33 2.50 2.00 1.67 1.43 1.25 1.11 1.00 0.91 0.83
P_{1} = 0.9  5.50 3.67 2.75 2.20 1.83 1.57 1.38 1.22 1.10 1.00 0.92
P_{1} = 1.0  6.00 4.00 3.00 2.40 2.00 1.71 1.50 1.33 1.20 1.09 1.00
Γ = 8  P_{1} = 0.0  1.00 0.59 0.42 0.32 0.26 0.22 0.19 0.17 0.15 0.14 0.12
P_{1} = 0.1  1.70 1.00 0.71 0.55 0.45 0.38 0.33 0.29 0.26 0.23 0.21
P_{1} = 0.2  2.40 1.41 1.00 0.77 0.63 0.53 0.46 0.41 0.36 0.33 0.30
P_{1} = 0.3  3.10 1.82 1.29 1.00 0.82 0.69 0.60 0.53 0.47 0.42 0.39
P_{1} = 0.4  3.80 2.24 1.58 1.23 1.00 0.84 0.73 0.64 0.58 0.52 0.48
P_{1} = 0.5  4.50 2.65 1.87 1.45 1.18 1.00 0.87 0.76 0.68 0.62 0.56
P_{1} = 0.6  5.20 3.06 2.17 1.68 1.37 1.16 1.00 0.88 0.79 0.71 0.65
P_{1} = 0.7  5.90 3.47 2.46 1.90 1.55 1.31 1.13 1.00 0.89 0.81 0.74
P_{1} = 0.8  6.60 3.88 2.75 2.13 1.74 1.47 1.27 1.12 1.00 0.90 0.83
P_{1} = 0.9  7.30 4.29 3.04 2.35 1.92 1.62 1.40 1.24 1.11 1.00 0.91
P_{1} = 1.0  8.00 4.71 3.33 2.58 2.11 1.78 1.54 1.36 1.21 1.10 1.00
Γ = 10 P_{1} = 0.0  1.00 0.53 0.36 0.27 0.22 0.18 0.16 0.14 0.12 0.11 0.10
P_{1} = 0.1  1.90 1.00 0.68 0.51 0.41 0.35 0.30 0.26 0.23 0.21 0.19
P_{1} = 0.2  2.80 1.47 1.00 0.76 0.61 0.51 0.44 0.38 0.34 0.31 0.28
P_{1} = 0.3  3.70 1.95 1.32 1.00 0.80 0.67 0.58 0.51 0.45 0.41 0.37
P_{1} = 0.4  4.60 2.42 1.64 1.24 1.00 0.84 0.72 0.63 0.56 0.51 0.46
P_{1} = 0.5  5.50 2.89 1.96 1.49 1.20 1.00 0.86 0.75 0.67 0.60 0.55
P_{1} = 0.6  6.40 3.37 2.29 1.73 1.39 1.16 1.00 0.88 0.78 0.70 0.64
P_{1} = 0.7  7.30 3.84 2.61 1.97 1.59 1.33 1.14 1.00 0.89 0.80 0.73
P_{1} = 0.8  8.20 4.32 2.93 2.22 1.78 1.49 1.28 1.12 1.00 0.90 0.82
P_{1} = 0.9  9.10 4.79 3.25 2.46 1.98 1.65 1.42 1.25 1.11 1.00 0.91
P_{1} = 1.0 10.00 5.26 3.57 2.70 2.17 1.82 1.56 1.37 1.22 1.10 1.00

Table 2. The Point Estimates and 95% Confidence Intervals for the Odds Ratio of PPH Associated With Long-Term Use of Appetite Suppressant Drugs With Adjustment for an

Unmeasured Binary Confounder of Specified Properties.

P_{0} = 0.0  0.1  0.2  0.3  0.4  0.5
Γ = 2 P_{1} = 0.0  23.1

(6.9, 77.7)

P_{1} = 0.2 19.3 21.2 23.1
(5.8, 64.8) (6.3, 71.2) (6.9, 77.7)
P_{1} = 0.4 16.5 18.2 19.8 21.5 23.1
(4.9, 55.5) (5.4, 61.1) (5.9, 66.6) (6.4, 72.2) (6.9, 77.7)
P_{1} = 0.6 14.4 15.9 17.3 18.8 20.2 21.7
(4.3, 48.6) (4.7, 53.4) (5.2, 58.3) (5.6, 63.1) (6.0, 68.0) (6.5, 72.8)
P_{1} = 0.8 12.8 14.1 15.4 16.7 18.0 19.3
(3.8, 43.2) (4.2, 47.5) (4.6, 51.8) (5.0, 56.1) (5.4, 60.4) (5.8, 64.8)
P_{1} = 1.0 11.6 12.7 13.9 15.0 16.2 17.3
(3.5, 38.9) (3.8, 42.7) (4.1, 46.6) (4.5, 50.5) (4.8, 54.4) (5.2, 58.3)
Γ = 6 P_{1} = 0.0  23.1

(6.9, 77.7)

P_{1} = 0.2 11.6 17.3 23.1
(3.5, 38.9) (5.2, 58.3) (6.9, 77.7)
P_{1} = 0.4 7.7 11.6 15.4 19.3 23.1
(2.3, 25.9) (3.5, 38.9) (4.6, 51.8) (5.8, 64.8) (6.9, 77.7)
P_{1} = 0.6 5.8 8.7 11.6 14.4 17.3 20.2
(1.7, 19.4) (2.6, 29.1) (3.5, 38.9) (4.3, 48.6) (5.2, 58.3) (6.0, 68.0)
P_{1} = 0.8 4.6 6.9 9.2 11.6 13.9 16.2
(1.4, 15.5) (2.1, 23.3) (2.8, 31.1) (3.5, 38.9) (4.1, 46.6) (4.8, 54.4)
P_{1} = 1.0 3.9 5.8 7.7 9.6 11.6 13.5
(1.2, 13.0) (1.7, 19.4) (2.3, 25.9) (2.9, 32.4) (3.5, 38.9) (4.0, 45.3)
Γ = 10 P_{1} = 0.0  23.1

(6.9, 77.7)

P_{1} = 0.2 8.3 15.7 23.1
(2.5, 27.8) (4.7, 52.7) (6.9, 77.7)
P_{1} = 0.4 5.0 9.5 14.1 18.6 23.1
(1.5, 16.9) (2.9, 32.1) (4.2, 47.3) (5.6, 62.5) (6.9, 77.7)
P_{1} = 0.6 3.6 6.9 10.1 13.4 16.6 19.9
(1.1, 12.1) (2.0, 23.1) (3.0, 34.0) (4.0, 44.9) (5.0, 55.8) (5.9, 66.8)
P_{1} = 0.8 2.8 5.4 7.9 10.4 13.0 15.5
(0.8, 9.5) (1.6, 18.0) (2.4, 26.5) (3.1, 35.1) (3.9, 43.6) (4.6, 52.1)
P_{1} = 1.0 2.3 4.4 6.5 8.5 10.6 12.7
(0.7, 7.8) (1.3, 14.8) (1.9, 21.8) (2.6, 28.7) (3.2, 35.7) (3.8, 42.7)

Note: P_{1} and P_{0} are, respectively, the prevalences of the unmeasured confounder among those had used appetite suppressants for more than 3 months and among those had never used them, and Γ is the odds ratio of PPH associated with the unmeasured confounder.

Table 3. The Point Estimates and 95% Confidence Intervals for the Relative Hazard of Death Associated With RHC Use With Adjustment for an Unmeasured Binary Confounder Having

Relative Hazard of Death of 3.

P_{0} = 0.0  0.1  0.2  0.3  0.4  0.5
P_{1} = 0.0  1.21
(1.09, 1.25)
P_{1} = 0.1 1.01 1.21
(0.91, 1.04) (1.09, 1.25)
P_{1} = 0.2 0.86 1.04 1.21
(0.78, 0.89) (0.93, 1.07) (1.09, 1.25)
P_{1} = 0.3 0.76 0.91 1.06 1.21
(0.68, 0.78) (0.82, 0.94) (0.95, 1.09) (1.09, 1.25)
P_{1} = 0.4 0.67 0.81 0.94 1.08 1.21
(0.61, 0.69) (0.73, 0.83) (0.85, 0.97) (0.97, 1.11) (1.09, 1.25)
P_{1} = 0.5 0.60 0.73 0.85 0.97 1.09 1.21
(0.55, 0.62) (0.65, 0.75) (0.76, 0.88) (0.87, 1.00) (0.98, 1.12) (1.09, 1.25)
P_{1} = 0.6 0.55 0.66 0.77 0.88 0.99 1.10
(0.50, 0.57) (0.59, 0.68) (0.69, 0.80) (0.79, 0.91) (0.89, 1.02) (0.99, 1.14)
P_{1} = 0.7 0.50 0.60 0.71 0.81 0.91 1.01
(0.45, 0.52) (0.55, 0.62) (0.64, 0.73) (0.73, 0.83) (0.82, 0.94) (0.91, 1.04)
P_{1} = 0.8 0.47 0..56 0.65 0.74 0.84 0.93
(0.42, 0.48) (0.50, 0.58) (0.59, 0.67) (0.67, 0.77) (0.75, 0.87) (0.84, 0.96)
P_{1} = 0.9 0.43 0.52 0.60 0.69 0.78 0.86
(0.39, 0.45) (0.47, 0.54) (0.55, 0.62) (0.62, 0.71) (0.70, 0.80) (0.78, 0.89)
P_{1} = 1.0 0.40 0.48 0.56 0.65 0.73 0.81
(0.36, 0.42) (0.44, 0.50) (0.51, 0.58) (0.58, 0.67) (0.65, 0.75) (0.73, 0.83)

Note: P_{1} and P_{0} are the prevalences of the unmeasured confounder among the RHC users and non-users, respectively.

Table 4. The Prevalences of the Unmeasured Confounder Among the Users and Non-users of RHC, P_{1} and P_{0}, and the Corresponding Point Estimates and 95% Confidence Intervals for the Relative Hazard of Death Associated With RHC, R̂ and 95% CI, Under a Direct Parameterization

for the Effect of an Unmeasured Binary Confounder on the Use of RHC.

Ψ = Γ = 2  Ψ = Γ = 3
P_{1} P_{0}  R̂  95% CI  P_{1} P_{0}  R̂  95% CI
P = 0.1  0.18 0.05 1.07 (0.97, 1.11)  0.25 0.01 0.82 (0.74, 0.85)
P = 0.2  0.33 0.12 1.01 (0.91, 1.05)  0.43 0.06 0.73 (0.66, 0.75)
P = 0.3  0.46 0.20 0.99 (0.90, 1.03)  0.56 0.14 0.73 (0, 66, 0.75)
P = 0.4  0.57 0.29 1.00 (0.90, 1.03)  0.67 0.24 0.76 (0.69, 0.79)
P = 0.5  0.67 0.40 1.01 (0.91, 1.05)  0.75 0.35 0.82 (0.74, 0.85)
P = 0.6  0.75 0.51 1.04 (0.94, 1.08)  0.82 0.47 0.89 (0.80, 0.92)
P = 0.7  0.82 0.62 1.08 (0.97, 1.11)  0.88 0.59 0.96 (0.87, 0.99)
P = 0.8  0.89 0.75 1.12 (1.01, 1.16)  0.92 0.72 1.04 (0.94, 1.08)
P = 0.9  0.95 0.87 1.16 (1.05, 1.20)  0.96 0.86 1.12 (1.01, 1.16)

Note: P is the overall prevalence of the unmeasured confounder; Ψ and Γ are, respectively, the relative risk of RHC and relative hazard of death associated with the unmeasured confounder. The observed prevalence of RHC, ϕ̂, was 0.38, with an estimated standard error of 0.006. We did not adjust for the variation of ϕ̂ in constructing the 95% confidence intervals, but the variation is too small to make an appreciable difference.