An approach to nonparametric inference on the causal dose–response function

2.1 Identification of the dose–response function

Let O 1 , … , O n be i.i.d. random vectors drawn from a probability distribution P 0 , which resides in a model ℳ . We allow ℳ to be a rich nonparametric model that is essentially unrestricted and must only satisfy some mild regularity conditions. We express our data as O = ( W , A , Y ) , where Y is a bounded real-valued outcome with sample space Y ⊂ R , A is a bounded exposure variable with sample space A ⊂ R , and W is a q -dimensional vector of covariates, with sample space W ⊂ R q . Throughout, we write E P to denote the expectation under any probability distribution P ∈ ℳ , and we use the shorthand notation E 0 ≔ E P 0 .

Let Y ( a ) denote the counterfactual outcome under exposure level a , or the potential outcome that would have been observed if an individual had been observed with exposure level A = a . We define the counterfactual mean outcome as the mapping P ↦ θ P ( a ) ≔ E P [ Y ( a ) ] for a ∈ A . Our objective is to perform inference on θ 0 ≔ θ P 0 , which is commonly called the causal dose–response function.

The stochastic process of counterfactual outcomes for subject i , { Y i ( a ) : a ∈ A } is not observable, as we only observe the outcome under one single exposure level A i . However, under some standard causal assumptions, the counterfactual mean outcomes can nonetheless be estimated from the data. Let Q P ( w , a ) ≔ E P [ Y ∣ W = w , A = a ] denote the conditional mean of Y given the exposure and covariates, and let g P ( a ∣ w ) ≔ d d a P ( A ≤ a ∣ W = w ) be the conditional density of the exposure given the covariates. Similarly as above, we use the shorthand notation Q 0 ≔ Q P 0 and g 0 ≔ g P 0 . The dose–response function can be estimated from the data under the following assumptions:

• Assumption A1 (Consistency). A = a implies Y = Y ( a ) .

• Assumption A2 (No unmeasured confoundedness). Y ( a ) is conditionally dependent on A , given W .

• Assumption A3 (Positivity). There exists g min > 0 so that inf a ∈ A g 0 ( a ∣ w ) > g min for all w ∈ W .

The consistency assumption states that the observed outcome for each subject is the potential outcome under their observed exposure level. The no unmeasured confoundedness assumption says that all confounding variables are contained within W . The positivity assumption means that it is possible for subjects with any covariate measurement to receive any treatment a ∈ A . When Assumptions A1–A3 hold, the counterfactual mean outcome at any exposure level a can be expressed as

where E 0 , W denotes the expectation with respect to the marginal distribution of W under P 0 .

Assumptions A1–A3 typically hold in randomized trials, but they are generally unverifiable in observational studies. However, even when the causal assumptions fail, θ 0 retains an interpretation as a mean regression function, which can be used to study the conditional association between the exposure and the outcome, given the covariates.

2.2 Summary of proposed methodology

We first derive a test of the null hypothesis that the dose–response function, centered about its mean, is equal to a given candidate null function. Let θ ¯ P : a ↦ θ P ( a ) − E P [ θ P ( A ) ] be the mean-centered dose–response function. For an arbitrary function θ * from A to R , let θ ¯ P * = θ * − E P [ θ * ( A ) ] , and let θ ¯ 0 * ≔ θ ¯ P 0 * . We are interested in testing the null hypothesis that θ ¯ 0 is equal to the candidate null parameter θ ¯ 0 * ,

Of particular interest in many applications is the case where θ ¯ * = 0 , which corresponds to the null hypothesis that the dose–response function is flat, or that the counterfactual mean outcomes are the same at every level of the exposure.

We also propose an approach for constructing confidence sets for the centered dose–response function. That is, for any α ∈ ( 0 , 1 ) , we construct a set C n α that contain θ ¯ 0 with probably at least 1 − α as n tends to infinity, i.e.,

When it is possible to obtain an estimator for θ 0 that has negligible asymptotic bias and a characterizable limiting distribution, it is straightforward to derive hypothesis tests and construct confidence sets based on the estimator. When the dose–response function is pathwise differentiable, meaning that θ P ( a ) changes smoothly in P with respect to local fluctuations around P 0 , there are well-established theory and methodology for constructing n 1 ⁄ 2 -consistent and asymptotically efficient estimators [21]. However, θ P ( a ) is typically only pathwise differentiable in a nonparametric model when the event { A = a } occurs with probability greater than zero. While this condition can be satisfied when A has a discrete support under P 0 , it is generally not met when the exposure variable is continuous.

While the dose–response function is not pathwise differentiable when the exposure is continuous, one can perform inference by instead estimating a collection of pathwise differentiable functionals of the dose–response function, as suggested by, e.g., Westling [18] and Hudson et al. [19]. In this work, we present an approach based on the estimation of linear transformations of the causal dose–response function. In what follows, we briefly summarize our approach to inference.

For any function h from A to R , we define

as the L 2 ( P ) inner product of θ ¯ P − θ ¯ P * and h , and let ψ 0 , θ * ( h ) ≔ ψ P 0 , θ * ( h ) . Observe that if θ ¯ 0 is equal to θ ¯ 0 * almost everywhere, ψ 0 , θ * ( h ) = 0 for all functions h , and if θ ¯ 0 is not almost everywhere equal to θ ¯ 0 * , there exists a function h * such that ψ 0 , θ * ( h * ) ≠ 0 . Therefore, one could test the null hypothesis in (1) by assessing whether, for a large class ℋ of bounded functions from A to R ,

or in other words, whether there exists a linear transformation of θ ¯ 0 − θ ¯ 0 * that is nonzero.

The manner in which the supremum in (2) is calculated will depend on the construction of ℋ . For instance, if ℋ is countable, one may approximate the supremum as the maximum of a finite approximation of { ψ 0 , h : h ∈ ℋ } . We also consider the setting in which ℋ is constructed using a basis expansion, and the supremum corresponds to an optimal weighting of the coefficients for the basis functions.

We develop a test of the null hypothesis based on the estimation of Ψ 0 , θ * ( ℋ ) . We later show that ψ 0 , θ * ( h ) is pathwise differentiable for any fixed h and can therefore be estimated at the parametric rate of n 1 ⁄ 2 . If, in addition, ℋ is not too complex, one can construct an estimator ψ n , θ * such that the standardized process { n 1 ⁄ 2 [ ψ n , θ * ( h ) − ψ 0 , θ * ( h ) ] : h ∈ ℋ } converges weakly to a Gaussian process G ≔ { G ( h ) : h ∈ ℋ } . As a consequence, we have that Ψ n , θ * ( ℋ ) ≔ sup h ∈ ℋ ∣ ψ n , θ * ( h ) ∣ is a consistent estimator for Ψ 0 , θ * ( ℋ ) , and under the null, when ψ 0 , θ * ( h ) = 0 for all h , n 1 ⁄ 2 Ψ n , θ * ( ℋ ) converges weakly to sup h ∈ ℋ ∣ G ( h ) ∣ . Because Ψ n , θ * ( ℋ ) has a tractable limiting distribution, it can be used as a test statistic. A hypothesis test that rejects the null when n 1 ⁄ 2 Ψ n , θ * ( ℋ ) is larger than the 1 − α quantile of the null limiting distribution would achieve the type-1 error level α asymptotically.

We can invert our proposed hypothesis test to obtain a confidence set for θ ¯ 0 . Let Θ be a large nonparametric class of functions from A to R that serve as candidate values for the dose–response function. Consider the set

We can interpret C n α as the set of functions in Θ that are in accordance with the observed data. If θ 0 belongs to Θ , then θ ¯ 0 belongs to C n α with probability tending 1 − α , and C n α is an asymptotically valid confidence set for θ ¯ 0 . Even when θ 0 does not belong to Θ , C n α retains a nice interpretation as long as Θ contains a good approximation for θ 0 .

The confidence set C n α can be difficult to visualize as it is a complex set of infinite-dimensional objects. We propose to simply summarize C n α by displaying the smallest band that contains all functions in C n α , which can be defined point-wise as

As we discuss later in Section 3.2, this confidence band is generally conservative for correctly specified Θ .

The inferential procedure described above requires selection of the class of functions ℋ , the elements of which index the set of linear transformations of θ ¯ 0 − θ ¯ 0 * that we need to estimate. Asymptotic type-1 error control is preserved using any choice ℋ that satisfies some mild regularity conditions, to be discussed later. The choice of ℋ does, however, influence the test’s statistical power. We use the following intuition for constructing an ℋ so that our test is well-powered. Because the inner product ψ 0 , θ * ( h ) is a measure of orthogonality of θ ¯ 0 − θ ¯ 0 * to a given h , the amount of evidence ψ 0 , θ * ( h ) provided against the null depends on the shape of h , rather than the scale. It is therefore sensible to construct ℋ as a class of equally scaled functions that contains the function that is least orthogonal θ ¯ 0 . We scale each h to have unit variance; while other measures of scale could be chosen, we scale by the standard deviation for simplicity. It can be shown by an application of the Cauchy-Schwarz inequality that the maximizer of ∣ ψ 0 , θ * ( h ) ∣ , among all h with unit standard deviation, can be expressed as

where Var 0 is the population variance under P 0 . The maximizer h 0 is a scaled difference between the true dose–response function and the candidate null parameter, and the maximal inner product ψ 0 , θ * ( h 0 ) is equal to the standard deviation of θ 0 ( A ) − θ 0 * ( A ) . We construct ℋ as a class of smooth functions that contain a good approximation of h 0 so that our target estimand Ψ 0 , θ * ( ℋ ) can be interpreted as an approximation of the standard deviation of the difference between the truth and the null.

The hypothesis testing procedure proposed by Westling [18] can be viewed as a special case of our method with ℋ taken as the class of binary indicators of whether an input exposure level a is greater than any given cutoff a 0 ∈ A . While Westling [18] showed that this choice results in an omnibus test of the null hypothesis, we later show that in finite samples, performance can be improved by considering a class of functions that contain the function that is least orthogonal to the difference between the true dose–response function and the null.

Our proposal is closely related to the nonparametric score test presented in the study of Hudson et al. [19]. The authors use the representation of function-valued nonpathwise differentiable parameters as the minimizer of a population risk functional to derive a set of estimating equations, indexed by functions h in a large class ℋ , that the true population parameter must satisfy. They then construct hypothesis tests that assesses whether the candidate null parameter also satisfies these estimating equations. Their proposal rejects the null hypothesis when the data provide contrary evidence. Our proposal can be viewed as a special case of this approach, where we assess whether θ * satisfies the estimating equations ψ 0 , θ * ( h ) = 0 for all h ∈ ℋ .

3.1 Estimation of ψ 0 , θ * ( h )

Recall that our proposal requires us to have an estimator ψ n ( h ) of ψ 0 ( h ) and a class of functions ℋ such that the process { n 1 ⁄ 2 [ ψ n ( h ) − ψ 0 ( h ) ] : h ∈ ℋ } converges weakly to a Gaussian process. In this subsection, we specify conditions on ψ n , θ * and ℋ such that the weak convergence property is satisfied, and we discuss the construction of a weakly convergent estimator.

As noted in Section 2, for any estimand that is pathwise differentiable in the sense of Bickel et al. [21], one can construct an estimator that, when centered around the true target estimand, converges weakly to a Gaussian distribution at the parametric rate of n 1 ⁄ 2 . Constructing such an estimator and establishing efficiency typically require knowledge of the efficient influence function of the estimand of interest. The following lemma states that ψ P , θ * ( h ) is pathwise differentiable in a nonparametric model and provides the form of the efficient influence function.

As before, we use the shorthand notation ϕ 0 , θ * ≔ ϕ P 0 , θ * to denote the efficient influence function at P 0 . Lemma 1 generalizes a result presented in [18] that characterizes the efficient influence function for the special case in which h is an indicator of whether the observed treatment level is greater than a specified cutoff, and θ ¯ 0 * = 0 .

Because ψ 0 , θ * ( h ) is pathwise differentiable, it is possible to construct an estimator ψ n , θ * ( h ) that is asymptotically linear in the sense that

where r n ( h ) = o P ( n − 1 ⁄ 2 ) . In view of the central limit theorem and the fact that ϕ 0 , θ * has zero mean and finite variance, an asymptotically linear estimator ψ n , θ * ( h ) is asymptotically Gaussian for any fixed h . However, our proposal requires a stronger notion of uniform convergence of ψ n , θ * ( h ) for h in a large collection ℋ , so some additional conditions are needed. The following lemma, which is a consequence of Slutsky’s theorem, provides conditions under which the desired uniform convergence holds.

The first condition of Lemma 2 is a constraint on the complexity of ℋ and typically holds when ℋ is a P 0 -Donsker class. The second condition directly involves the estimator ψ n and requires that the remainder term is asymptotically negligible in a uniform sense.

In what follows, we present two strategies for constructing weakly convergent estimators for { ψ 0 , θ * ( h ) : h ∈ ℋ } . We begin by considering a naïve plug-in estimator of ψ 0 ( h ) . Suppose we have a consistent estimator P ˆ n for P 0 . In practice, we do not need to estimate the entire distribution P 0 , and we must only estimate nuisance parameters upon which ψ 0 , θ * and ϕ 0 , θ * depend. In our setting, there are four nuisance parameters: (i) F 0 , W , the marginal cumulative distribution for W ; (ii) F 0 , A , the marginal cumulative distribution function for A ; (iii) Q 0 , the conditional mean of Y given A and W ; and (iv) g 0 , the conditional density of A given W . One can obtain estimators F n , W and F n , A for F 0 , W and F 0 , A nonparametrically using the empirical distribution function, and one typically requires machine learning to construct consistent nonparametric estimators Q n and g n for Q 0 and g 0 . Given the initial estimator P ˆ n , one can obtain the naïve plug-in estimator ψ P ˆ n ( h ) , which can be expressed as

where θ n : a ↦ n − 1 ∑ i = 1 n Q n ( a , W i ) is the plug-in estimator for the dose–response function.

The plug-in estimator ψ P ˆ n , θ * ( h ) will typically retain non-negligible asymptotic bias for ψ 0 , θ * ( h ) and consequently will not be asymptotically linear. This bias is attributable to the fact that nonparametric estimators Q n of Q 0 are usually obtained by balancing a bias-variance tradeoff that is sub-optimal for the objective of estimating ψ 0 , θ * ( h ) . We discuss two widely used strategies for correcting the bias of the naïve plug-in: one-step estimation [22] and TML-based estimation [23,24].

The estimation strategies we discuss require the following assumptions:

Assumption B1. There exists a P 0 -Donsker class Φ that contains ϕ 0 , θ * ( ⋅ ; h ) and ϕ P ˆ n , θ * ( ⋅ ; h ) for each h ∈ ℋ with probability tending to one.

Assumption B2. The nuisance parameter estimators satisfy

Assumptions B1 and B2 impose conditions on the estimators for the nuisance parameters upon which the target estimand and the efficient influence function depend. Assumption B1 places a constraint on the complexity of the family of candidate estimators for the nuisance components. Many flexible nonparametric estimators, e.g., those constructed via the highly adaptive LASSO [25], satisfy this condition. Assumption B2 places a requirement on the rates of convergence that the conditional mean and conditional density estimators must achieve. This condition holds when both estimators are consistent, and the product of the convergence rates is greater than n 1 ⁄ 2 . In view of Assumption B2, the estimators of ψ 0 , θ * ( h ) we develop are doubly robust in the sense that consistency and asymptotic normality are achieved if one of the nuisance parameters is estimated at a slow rate as long as the other nuisance is estimated at a fast enough rate to compensate.

We first construct a one-step estimator for ψ 0 , θ * ( h ) . Given an estimator for the nuisance parameters upon which ϕ 0 , θ * ( h ) depends, we can obtain an estimator ϕ P ˆ n , θ * ( ⋅ ; h ) for the efficient influence function ϕ 0 , θ * ( ⋅ ; h ) . The empirical average of the estimator of the influence function ϕ P ˆ n , θ * can be shown to serve as a first-order approximation to the bias of the naïve plug-in [22]. This allows one to perform a, so-called, one-step bias correction. We define the one-step estimator as

The following theorem states that, under mild regularity conditions, the one-step estimator is asymptotically linear and hence asymptotically Gaussian.

While one-step estimators are asymptotically efficient, they are not guaranteed to be compatible in the sense that there exists a probability distribution P in the model ℳ such that ψ P , θ * ( h ) = ψ n , θ * I ( h ) for all h in ℋ . TML-based estimation is an appealing alternative strategy that can be used to construct an estimator for { ψ 0 , θ * ( h ) : h ∈ ℋ } that is compatible in this sense. TML estimators correct the bias of the naïve plug-in by updating the initial estimator P ˆ n of P 0 to obtain a new estimator P ˜ n such that the updated plug-in { ψ P ˜ n , θ * : h ∈ ℋ } that takes as input P ˜ n has reduced bias for the target estimand { ψ 0 , θ * ( h ) : h ∈ ℋ } . Therefore, as long as P ˜ n resides within ℳ , the TML estimator is compatible. In our presentation, we only briefly summarize some of the main principles of targeted learning, and we refer readers to previous studies [23,24] for a comprehensive discussion.

The main idea behind TML-based estimation is to construct an updated estimator P ˜ n based on the initial estimator P ˆ n so that the following efficient influence function estimating equations are satisfied:

and such that P ˜ n remains sufficiently close to P ˆ n , so as to remain a good estimator for P 0 . Because estimating a marginal distribution function using the empirical distribution function does not generate bias for the target estimand, we do not need to update the initial estimators F A , n and F W , n for F A , 0 and F W , 0 . We only need to obtain an updated estimator Q ˜ n for the conditional mean Q 0 , since the initial estimator makes a bias-variance trade-off that is suboptimal for the estimation of the target parameter.

It can be verified algebraically that for any choice Q ˜ n , the empirical average of the efficient influence function evaluated at P ˜ n can be expressed as

where we define Z n ( w , a ; h ) as

Thus, Q ˜ n satisfies the efficient influence function estimating equations at an adequate level if

We now discuss how to construct a Q ˜ n that satisfies (9). Let Q n , β be a parametric working model indexed by a scalar parameter β ∈ R , for which Q n , β = Q n when β = 0 . We construct the working model so that the derivative of the squared error loss is equal in magnitude to the supremum over ℋ of the empirical average of ϕ P ˜ n ( ⋅ ; h ) with Q ˜ n = Q n , β , i.e.,

We then take Q ˜ n = Q n , β n , where β n is a near minimizer of the squared error loss and satisfies

for a small positive sequence δ n ↓ 0 . Because β ˜ n is a near minimizer of the loss, we can see that this choice of Q ˜ n satisfies (9) for δ n sufficiently small. We note that a sub-model satisfying (10) is referred to as a universal least favorable submodel and provides the maximal reduction of the bias of { ψ P ˜ n , θ * ( h ) : h ∈ ℋ } as Q n , β moves away from Q n toward Q n , β n . Strategies based on a locally least favorable submodel, which would satisfy (10) only when β = 0 , could alternatively have been considered, but such approaches tend to perform worse in small samples, in particular when the target estimand is multidimensional or infinite-dimensional [26]. We recursively define the universal least favorable sub-model Q n , β point-wise as

where h n , β is a solution to

It can be verified using the fundamental theorem of calculus that the above working model satisfies (10).

After obtaining the updated estimator Q ˜ n , we can construct the TML estimator as

where θ P ˜ n : a ↦ n − 1 ∑ i = 1 n Q ˜ n ( W i , a ) is the updated TML estimator for the dose–response function. The following theorem states that the TML estimator satisfies the conditions of Lemma 2.

3.2 Inference on Ψ 0 , θ * ( ℋ )

We are at this point prepared to discuss inference on Ψ 0 , θ * ( ℋ ) = sup h ∈ ℋ ∣ ψ 0 , θ * ( h ) ∣ . Suppose that we have an estimator { ψ n , θ * ( h ) : h ∈ ℋ } for { ψ 0 , θ * ( h ) : h ∈ ℋ } and a function class ℋ that satisfy the conditions of Lemma 2. Such an estimator could be obtained using either of the strategies presented in Section 3.1. Consider the plug-in estimator Ψ n , θ * ( ℋ ) ≔ sup h ∈ ℋ ∣ ψ n , θ * ( h ) ∣ . The continuous mapping theorem implies that n 1 ⁄ 2 sup h ∈ ℋ ∣ ψ n , θ * ( h ) − ψ 0 , θ * ( h ) ∣ converges weakly to sup h ∈ ℋ ∣ G ( h ) ∣ where G is the Gaussian process in Lemma 2. This in combination with the reverse triangle inequality together imply that

Furthermore, because when the null hypothesis holds, ψ 0 , θ * ( h ) = 0 for all h , n 1 ⁄ 2 Ψ n , θ * ( ℋ ) converges weakly to sup h ∈ ℋ ∣ G ( h ) ∣ . Thus, Ψ n , θ * ( ℋ ) is a consistent estimator for Ψ 0 , θ * ( ℋ ) that has a fully characterizable null limiting distribution. This makes it possible to construct an asymptotically valid hypothesis test based on the estimator.

We make this explicit by providing the form of our proposed test below. Letting t 1 − α * denote the 1 − α quantile of the null limiting distribution of n 1 ⁄ 2 Ψ n , θ * ( ℋ ) , we consider the hypothesis test

This test achieves type-1 error control at the level α when the conditions of Lemma 2 are satisfied. Additionally, when ℋ contains a function that is nonorthogonal to θ ¯ 0 * (so that Ψ 0 , θ * ( ℋ ) > 0 ), the above test is consistent in the sense that the probability of rejecting the null tends to one. To see this, we observe that n 1 ⁄ 2 Ψ n , θ * ( ℋ ) tends to infinity because Ψ n , θ * ( ℋ ) converges to Ψ 0 , θ * ( ℋ ) in probability.

Moreover, the fact that our proposed test achieves nominal type-1 error control asymptotically implies that the confidence band in (4) is a uniform confidence band. This follows from the facts that

and that P 0 ( θ ¯ 0 ∈ C n α ) approaches 1 − α in the limit of large n when Θ contains θ 0 . This confidence band will in general achieve a coverage rate not lower than 1 − α , though this lower bound is not tight. This type of confidence band construction has been shown to result in over-coverage by Hudson et al. [19].

While the null limiting distribution of Ψ n , θ * ( ℋ ) can indeed be characterized, a closed form expression may not be available. It may therefore be necessary to use an approximation. We use the multiplier bootstrap method presented in the study of Hudson et al. [19], which makes use of the asymptotic linearity of ψ n , θ * ( h ) . We first note that due to the uniform asymptotic linearity of ψ n , θ * , under the null hypothesis, Ψ n , θ * ( ℋ ) can be expressed as the sum of the supremum of an empirical process and an asymptotically negligible remainder. That is,

We can approximate the null distribution of Ψ n , θ * ( ℋ ) as the supremum of a bootstrapped empirical process that attains the same limiting distribution as the empirical process in (13), conditional on the observed data. For m = 1 , … , M and M large, let ξ 1 m , … , ξ n m be i.i.d. random variables, independent of O 1 , … , O n , with mean zero, unit variance, and E 0 [ ∣ ξ i m ∣ 2 + u ] < ∞ for some u > 0 , and let ϕ n , θ * be an estimator for the efficient influence function. We define the m th bootstrap sample of Ψ n , θ * ( ℋ ) as