Role of placebo samples in observational studies

2.1 Nonparametric identification

Consider a binary exposure A ∈ { 0 , 1 } and potential outcomes { Y ( 0 ) , Y ( 1 ) } . Throughout the article, we assume the consistency assumption and stable unit treatment value assumption (SUTVA) so that the observed outcome Y satisfies Y = A Y ( 1 ) + ( 1 − A ) Y ( 0 ) [21]. Suppose that we have a binary baseline covariate S such that S = 0 represents the placebo sample that is unaffected by the exposure. We call subjects with S = 1 the primary sample.

Denote the other baseline covariates as X . Suppose that a random sample of n subjects is obtained, which is written as { ( Y i , A i , X i , S i ) : i = 1 , … , n } and assumed to be independent and identically distributed according to the joint law of ( Y ( A ) , A , X , S ) . Our target parameter is the average treatment effect of the treated (ATT) in the primary sample

Assumption 2 is our key identification assumption that links the subjects with S = 0 and S = 1 . Figure 1 illustrates the idea.

Figure 1

Illustration of Assumption 2.

The left-hand side of Assumption 2 encodes the level of confounding bias. Because the observed covariates, X may fail to render the exposure and control groups in the S = 1 stratum comparable, E { Y ( 0 ) ∣ S = 1 , A = 1 , X } ≠ E { Y ( 0 ) ∣ S = 1 , A = 0 , X } in general. However, Assumption 2 states that the extent of residual confounding bias in the S = 1 stratum is precisely equal to that in S = 0 , making it possible to debias using the placebo sample.

A similar assumption also appears in the DID and NCO literature. With repeated cross-sectional data, DID assumes a parallel trends assumption, which can be thought of as a special case of Assumption 2 that uses the pre- and post-exposure time indicator as S . As such, DID with repeated cross-sectional data can be recast as a placebo sample approach that uses a sample collected before the onset of the treatment as a placebo sample. For the NCO literature, the additive equi-confouding assumption says that the residual confounding bias for an NCO and the outcome is the same [12], compared to which Assumption 2 may be more reasonable as it is about the same outcome and is invariant to scaling of the outcome variable. Additional discussion of the assumptions can be found in Supplement §S1.2.

Consistent with the literature on DID, our focus is on the ATT, which can be identified without imposing restrictions on treatment effects for the primary sample. The average treatment effect is not identifiable without additional assumptions;  therefore, it is not considered in this article.

Finally, we make a standard positivity assumption.

We give some intuition before formally stating the identification results. For s = 0 , 1 , a = 0 , 1 , and every x , let μ Y ( s , a , x ) = E { Y ∣ S = s , A = a , X = x } , and Δ s ( x ) = μ Y ( s , 1 , x ) − μ Y ( s , 0 , x ) be observed data functions. The target causal parameter of interest can be decomposed into a contrast term C ( X ) and a bias term B ( X ) as follows:

By the consistency assumption, the contrast function C ( X ) equals Δ 1 ( X ) and is identified from observed data. By Assumption 2, the bias function B ( X ) equals E { Y ( 0 ) ∣ S = 0 , A = 1 , X } − E { Y ( 0 ) ∣ S = 0 , A = 0 , X } , which then equals Δ 0 ( X ) by Assumption 1. Finally, the target parameter E { Y ( 1 ) − Y ( 0 ) ∣ S = 1 , A = 1 , X } equals Δ 1 ( X ) − Δ 0 ( X ) and is identified from observed data. Alternatively, identification can be obtained from inverse probability weighting (IPW) that avoids modeling the outcome distribution. These two identification results are stated in Proposition 1.

All proofs in this article can be found in Supplement §S2. There are interesting connections among the DID, NCO, and placebo sample methods. Both NCO and placebo sample methods exploit the known null effect: the former leverages a placebo outcome known not to be affected by the treatment, whereas the latter utilizes a placebo population known not to be affected by the treatment. We view these two methods as complementary devices that can even be used in the same study. For instance, in the study of EITC on deaths of despair, Dow et al. [20] used both devices – a cancer outcome as an NCO and college graduates as a placebo sample. In another example, to study the causal effect of ivacaftor on lung functions in cystic fibrosis patients, Newsome et al. [22] also utilized both devices – outcomes measured in the pre-ivacaftor period as NCOs and patients ineligible for ivacaftor due to their genotype as a placebo sample. In addition, the DID with longitudinal data can be interpreted as an NCO method that uses the pre-exposure outcome as an NCO, and the DID with repeated cross-sectional data can be interpreted as a placebo sample method that uses a sample drawn before the exposure is in place as the placebo sample.

We note that Proposition 1 (b) slightly generalizes the IPW method proposed by Abadie [23], which is applicable only when S is independent of all the other variables.

2.2 Estimation and semiparametric inference

Proposition 1 suggests two estimators. The first one is a regression-based estimator

where n 11 = ∑ i = 1 n I ( S i = 1 , A i = 1 ) , Δ s ( x ; β ˆ ) = μ Y ( s , 1 , x ; β ˆ ) − μ Y ( s , 0 , x ; β ˆ ) , μ Y ( s , a , x ; β ) is a model of μ Y ( s , a , x ) parameterized by the finite-dimensional parameter β , and β ˆ is an estimator of β . The regression-based estimator θ ˆ reg is consistent for the target parameter when outcome regression models μ Y ( s , a , x ) are correctly specified by μ Y ( s , a , x ; β ) . In practice, researchers often specify generalized linear regression models for these outcome models.

A second estimator is an IPW estimator

where π S ( x ; ψ ) is a model of π S ( x ) parametrized by finite-dimensional parameter ψ , π A ( x , s ; α ) a model of π A ( x ) parametrized by α , and ψ ˆ and α ˆ estimators of ψ and α . In practice, researchers often fit logistic regression models for π S ( x ; ψ ) and π A ( x ; ψ ) . IPW-type estimators could be unstable when the (estimated) π A ( ⋅ ) and 1 − π S ( ⋅ ) are close to zero, and a common way to stabilize the weights is by normalization [24].

The reliability of the regression-based estimator and the IPW estimator depends on correct specification of different parts of the likelihood. When we are uncertain about which models are correctly specified, it is of interest to develop a doubly robust estimator that is guaranteed to be consistent and deliver valid inference about θ 0 , provided that either { μ Y } or { π S , π A } , but not necessarily both, are correctly specified. The next theorem derives the efficient influence function for θ [25,26] in the nonparametric model, where Assumptions 1–3 are assumed for identification, but the model is otherwise unrestricted. Similar to the discussion in Park and Tchetgen Tchetgen [27] regarding DID, our Assumptions 1 and 2 do not have any testable implications on the observed data (except for a falsification condition if the outcome is known to satisfy certain support conditions). Theorem 1 also provides the basis for constructing a doubly robust estimator.

The efficient influence function gives an estimator θ ˆ dr defined as the solution to ∑ i = 1 n EIF ( O i ; θ , η ˆ ) = 0 , where η ˆ = ( μ ˆ Y , π ˆ A , π ˆ S , λ ˆ ) denotes the collection of nuisance parameters, and λ ˆ = n 11 ⁄ n estimates E { S A } . We prove in Supplement §S2 that θ ˆ dr is doubly robust.

Next, we derive the asymptotic properties of θ ˆ dr . Let ‖ f ‖ 2 = { ∫ f 2 ( o ) d P ( o ) } 1 ⁄ 2 denote the L 2 ( P ) norm of any real-valued function f , and ‖ f ‖ 2 = ∑ j = 1 ℓ ‖ f j ‖ 2 for any collection of real-valued functions f = ( f 1 , … , f ℓ ) , where P denotes the distribution of O . Moreover, let η 0 = ( μ Y 0 , π A 0 , π S 0 , λ 0 ) denote the true values of the nuisance parameters.

Assumption 4 (a) describes the double robustness of our proposed estimator. Assumption 4 (b) is standard for M-estimators [26, Chapter 5.4].

Theorem 2 summarizes the doubly robust and locally efficient property of θ ˆ dr .

The first part of Theorem 2 characterizes the convergence rate of θ ˆ dr ;  it also implies that θ ˆ dr is doubly robust in the sense that it is consistent when either { μ ˆ Y } or { π ˆ S , π ˆ A } is consistent. The second part of Theorem 2 says that if the nuisance parameters are consistently estimated with fast rate, e.g., if they are estimated using parametric methods, then their variance contributions are negligible, and θ ˆ dr achieves the semiparametric efficiency bound. The results in Theorems 1 and 2 generalize the results in the study of Sant-Anna and Zhao [28] to allow for S being correlated with A and X .

When (2) holds, a plug-in variance estimator for n θ ˆ dr can be easily constructed as n − 1 ∑ i = 1 n EIF ( O i ; θ ˆ dr , η ˆ ) 2 . Even if (2) does not hold, e.g., when only the model for { μ Y } or the models for { π S , π A } are correctly specified, but all the nuisance parameters are finite-dimensional and in the form of M-estimators, n θ ˆ dr is still consistent and asymptotically normal from standard M-estimation theory [29, Chapter 6]. Thus, a consistent variance estimator for n θ ˆ dr can be constructed under the M-estimation framework, allowing for doubly robust inference (see details in Supplement §S3). Alternatively, the nonparametric bootstrap is commonly used in practice. However, when nonparametric, data-adaptive estimators of these nuisance parameters are used, conducting doubly robust inference becomes much more challenging if one of the nuisance parameters is inconsistently estimated. In such cases, θ ˆ dr may become irregular and exhibit a slower convergence rate than the typical root- n . Recent advancements in addressing these challenges can be found in previous works [30,31].

Finally, we note that the Donsker condition in Assumption 4 can be relaxed by using sample splitting [32], which enables estimating the nuisance parameters using complex machine learning methods that do not satisfy the Donsker condition. For example, even without the Donsker condition in Assumption 4, the result in (2) holds when the nuisance parameters are estimated at faster than n − 1 ⁄ 4 -rates.

3.1 Sensitivity analysis under linear sensitivity models

We propose a sensitivity analysis method of Assumptions 1–2 under linear sensitivity models. Recall in Section 2.1, we decompose the target parameter as θ 0 = E { C ( X ) − B ( X ) ∣ S = 1 , A = 1 } , where C ( X ) = Δ 1 ( X ) is identified from the observed data. The key is then to bound B ( X ) when Assumptions 1–2 are violated. Consider the following sensitivity model.

Sensitivity Model 1.A E { Y ( 1 ) − Y ( 0 ) ∣ S = 0 , A = 1 , X } = λ 0 + λ X T X almost surely.

Sensitivity Model 1.B E { Y ( 0 ) ∣ S = 1 , A = 1 , X } − E { Y ( 0 ) ∣ S = 1 , A = 0 , X } − [ E { Y ( 0 ) ∣ S = 0 , A = 1 , X } − E { Y ( 0 ) ∣ S = 0 , A = 0 , X } ] = δ 0 + δ X T X almost surely.

Here, H ≔ { λ 0 , λ X , δ 0 , δ X } are sensitivity parameters. When all these parameters equal zero, Sensitivity Models 1.A and 1.B degenerate to Assumptions 1–2 in the primary analysis. Under Sensitivity Models 1.A and 1.B, we have that θ 0 = θ 0 ( H ) , where

Fix parameters Γ L , Γ U , Λ L , and Λ U and consider the set of sensitivity parameters

where p = dim ( X ) . Then, θ 0 ( H ) is in the set

The set Θ 0 is often referred to as the partially identified set under Sensitivity Models 1.A and 1.B and parameters Γ L , Γ U , Λ L , and Λ U .

A general method to construct confidence intervals for the set Θ 0 is the union method. Taking an estimator θ ˆ ∈ { θ ˆ reg , θ ˆ ipw , θ ˆ dr } of E { Δ 1 ( X ) − Δ 0 ( X ) ∣ S = 1 , A = 1 } , an estimator of θ 0 ( H ) is

Suppose that for any fixed H , we construct a confidence interval [ L ( H ) , U ( H ) ] for θ 0 ( H ) , e.g., by nonparametric bootstrap, the interval [ inf H ∈ ℋ L ( H ) , sup H ∈ ℋ U ( H ) ] is an asymptotic confidence interval of θ 0 with at least ( 1 − α ) -coverage under the collection of sensitivity models ℋ . However, as discussed in the study of Zhao et al. [33], L ( H ) and U ( H ) may be complicated functions of H , and optimization over H may be complex. Therefore, we use

as our confidence interval, where inf H ∈ ℋ θ ˆ * ( H ) is the analog of inf H ∈ ℋ θ ˆ ( H ) computed using each bootstrap sample, and Q α ⁄ 2 ( inf H ∈ ℋ θ ˆ * ( H ) ) is the α ⁄ 2 -percentile of inf H ∈ ℋ θ ˆ * ( H ) in the bootstrap distribution; the quantities in the upper bound are defined similarly. The interval in (3) is shown by Qingyuan Zhao et al. [33] to cover θ 0 with probability at least 1 − α asymptotically. Note that as θ ˆ * ( H ) is monotone in every element in H , the confidence interval (3) is computationally efficient.

3.2 Sensitivity analysis under marginal sensitivity models

We consider a second sensitivity model. With a slight abuse of notation, we still denote the sensitivity parameters as Λ and Γ , although they have distinct meanings compared to those in Section 3.1:

Sensitivity Model 2.A − Λ ≤ E { Y ( 1 ) − Y ( 0 ) ∣ S = 0 , A = 1 , X } ≤ Λ almost surely.

Sensitivity Model 2.B Let U denote the set of unmeasured confounders such that A ⊥ ( Y ( 0 ) , Y ( 1 ) ) ∣ U , X , S . Suppose that (i)

where π S ( U , A , X ) = P ( S = 1 ∣ U , A , X ) , π S ( A , X ) = P ( S = 1 ∣ A , X ) , and OR ( p 1 , p 2 ) = p 1 ⁄ ( 1 − p 1 ) p 2 ⁄ ( 1 − p 2 ) ; and (ii)

is a function of X alone.

Sensitivity models 2.A and 2.B allow the treatment to have a small effect on the placebo sample and that the distribution of U being dependent on S after conditioning on A and X . In particular, the model in (4) is termed the marginal sensitivity model in the literature [33]. When Γ = 1 , Sensitivity Model 2.B reduces to the case where U ⊥ S ∣ A , X . Finally, we suppose that E { Y ( 0 ) ∣ S = 0 , A , X } ≥ 0 almost surely, which is not a restriction because we can always add a constant to every Y i to make Y i ≥ 0 without affecting the parameter of interest.

Let r ( A , U , X ) = P ( U ∣ S = 1 , A , X ) ⁄ P ( U ∣ S = 0 , A , X ) , V ( X ) = E { Y ( 0 ) ∣ S = 1 , A = 1 , X } − E { Y ( 0 ) ∣ S = 1 , A = 0 , X } − [ E { Y ( 0 ) ∣ S = 0 , A = 1 , X } − E { Y ( 0 ) ∣ S = 0 , A = 0 , X } ] , and T ( X ) = E { Y ( 1 ) − Y ( 0 ) ∣ S = 0 , A = 1 , X } . Under the Sensitivity Model 2.A, we have T ( X ) ∈ [ − Λ , Λ ] . From the Bayes formula, we have r ( A , U , X ) = OR { π S ( U , A , X ) , π S ( A , X ) } , and thus, Γ − 1 ≤ r ( A , U , X ) ≤ Γ under the Sensitivity Model 2.B. We show in Supplement §2 that

Finally, with θ 0 written as θ 0 = E [ Δ 1 ( X ) − Δ 0 ( X ) ∣ S = 1 , A = 1 ] + E [ T ( X ) − V ( X ) ∣ S = 1 , A = 1 ] , we can construct the bounds for θ 0 as [ θ L , θ U ] , where θ U = E [ Δ 1 ( X ) − Δ 0 ( X ) ∣ S = 1 , A = 1 ] + Λ − E [ E { Y ∣ S = 0 , A = 1 , X } ∣ S = 1 , A = 1 ] ( Γ − 1 − 1 ) + E [ E { Y ∣ S = 0 , A = 0 , X } ∣ S = 1 , A = 1 ] ( Γ − 1 ) and θ L = E [ Δ 1 ( X ) − Δ 0 ( X ) ∣ S = 1 , A = 1 ] − Λ − E [ E { Y ∣ S = 0 , A = 1 , X } ∣ S = 1 , A = 1 ] ( Γ − 1 ) + E [ E { Y ∣ S = 0 , A = 0 , X } ∣ S = 1 , A = 1 ] ( Γ − 1 − 1 ) . In correspondence to the three estimation approaches, the bounds in the sensitivity analysis can also be estimated in three ways (see Supplement §1 for their expressions).