Data-Adaptive Bias-Reduced Doubly Robust Estimation

1.1 Some relevant literature overview

The discussions on the Kang and Schafer article [5, 11–15] revealed that many different doubly robust estimators may be obtained for a given target parameter by varying the choice of nuisance working model estimators. Different choices may lead to doubly robust estimators which all have potentially very different behavior and properties under misspecification of at least one working model but are asymptotically equivalent under correct specification of both working models. As a result, many alternatives to maximum likelihood estimation (MLE) have now been proposed in the literature for the estimation of the nuisance parameters/nuisance working models that index certain doubly robust estimators. Rubin and van der Laan [16], Cao et al. [17] and Tsiatis et al. [18] developed doubly robust estimators in specific missing data models which target minimal asymptotic variance when the missingness model is correctly specified. The TMLE (targeted maximum likelihood estimation) procedure [9, 19] and the procedures of Tan [20] and Rotnitzky et al. [21] guarantee doubly robust estimators that fall within the parameter range, with the latter procedures also having desirable efficiency properties. With the exception of TMLE, all these proposals focus on improving the efficiency of doubly robust estimators under misspecification of the working model for the dependence of the observed outcome on covariates/confounders in missing data/causal inference models. The collaborative TMLE (C-TMLE) procedure of van der Laan and Gruber [22] additionally focuses on the estimation of the missingness/exposure probabilities in a way that aims to attain an optimal bias-variance trade-off. In contrast, Vermeulen and Vansteelandt [8] focus on bias reduction under the constraints of parametric working models. For a review of many of the aforementioned procedures, we refer to Rotnitzky and Vansteelandt [23]; for a comparison of the performance of many of these alternatives, we refer to Tan [20], Porter et al. [24] and Vermeulen and Vansteelandt [8].

3.1 Parametric nuisance working models and MLE

Suppose we parameterize the working model for the conditional mean outcome by a q-dimensional parameter γ: qγ={Qˉ(γ)(X)|γ∈Rq} with Qˉ(γ)(X)=Q{γTk(X)}, Q an appropriate inverse link function and k=(1,k1,…,kq−1); e.g., a linear regression model Qˉ(γ)(X)=γ1+γ2TX for a continuous outcome Y. If the model m(qγ) is correctly specified and thus includes P0, we let γ0 be such that Qˉ(γ0)=Qˉ0. Further, let the p-dimensional parameter α parameterize the working model for the missingness mechanism: gα={g(α)(X)|α∈Rp} where g(α)(X)=G{αTl(X)}, G is an appropriate inverse link function and l=(1,l1,…,lp−1); e.g., a logistic regression model g(α)(X)=expit(α1+α2TX), expit(x)=1/(1+e−x). If the model m(gα) is correctly specified and thus includes P0, we let α0 be such that g(α0)=g0.

Root-n consistent and asymptotically normal estimators ˆγn and αˆn for the nuisance parameters γ and α can be obtained as solutions to estimating equations PnDQˉ(ˆγn)=0 and PnDg(αˆn)=0 with the estimating functions DQˉ and Dg such that P0{DQˉ(γ0)}=0 if P0∈m(qγ) and P0{Dg(α0)}=0 if P0∈m(gα).

Throughout, we will assume that ˆγn→pγ∗ for some γ∗ (with γ∗=γ0 under model m(qγ)) and αˆn→pα∗ for some α∗ (with α∗=α0 under model m(gα)) and, moreover, that these estimators are asymptotically linear with influence function −P0{Dγ,Qˉ(γ∗)}−1DQˉ(γ∗) and −P0Dα,g(α∗)−1Dg(α∗), where Dγ,Qˉ(γ∗)=∂DQˉ(γ)/∂γ|γ=γ∗ and Dα,g(α∗)=∂Dg(α)/∂α|α=α∗.

In practice, maximum likelihood estimation (MLE) (or least squares) is routinely employed to estimate the nuisance parameters. The MLEs ˆγnMLE and αˆnMLE solve the estimating equations PnDQˉMLE(ˆγnMLE)=0 and PnDgMLE(αˆnMLE)=0 with

3.2 Bias-reduced doubly robust estimation

3.2.2 Illustration: estimation of a population mean with incomplete data

We illustrate the bias-reduced doubly robust estimation principle for the missing data problem discussed in Section 2. Consider the working models qγ and gα from Section 3.1 but with k and l of the same dimension (i.e., such that q=p). Estimators for the nuisance parameters (ˆγnBR,T,αˆnBR,T)T are then obtained by solving (7) and (8):

(9)n−1∑i=1nRi{Yi−Q¯(γ^nBR)(Xi)}G'{α^nMLE,Tl(Xi)}G2{αnMLE,Tl(Xi)}l(Xi)=0,

(10)n−1∑i=1n{1−Rig(α^nBR)(Xi)}Q'{γ^nBR,Tk(Xi)}k(Xi)=0.

Vermeulen and Vansteelandt [8] show how (10) can be solved by optimizing an integrated estimating equation; (9) can be solved via weighted least squares based on the complete cases. They moreover discuss various desirable properties of the doubly robust estimator based on these nuisance parameter estimators. Here, we develop visual insight using two simple examples.

Example 1

Consider a sample o of size n=105, which we take to be large so that we can ignore sampling variability. In this first example, for each individual i, we let Xi=dN(0,1), Ri|Xi=dBer{π0(Xi)} with π0(Xi)=expit(−1+Xi3) and Yi|Xi=dN(Xi2,1). The true mean outcome equals ψ0=1. One-dimensional working models of the form g(α)(X)=expit(αX) and Qˉ(γ)(X)=γX are misspecified. We obtain αˆnMLE=1.115 and γˆnMLE=1.606, leading to ψˆnDR,MLE=0.237 and bias(αˆnMLE,γˆnMLE;1)=−0.763. The bias-reduced doubly robust estimation strategy yields αˆnBR=2.254 and γˆnBR=−0.959, leading to ψˆnBR=0.663 and a smaller bias of bias(αˆnBR,γˆnBR;1)=−0.337, as theoretically expected.

Figure 1 shows a contour plot of the log of the squared asymptotic bias as a function of the nuisance parameters α and γ. It thus indicates for each combination of fixed nuisance parameter values γ and α the magnitude of the asymptotic bias of the doubly robust estimator for these choices of γ and α. Darker values indicate bias closer to zero. The cross “×” indicates the point (αˆnMLE,γˆnMLE)=(1.115,1.606), which approximates (αMLE∗,γMLE∗) and the bullet “∙” indicates the point (αˆnBR,γˆnBR)=(2.254,−0.959), which approximates (αBR∗,γBR∗). Figure 1 illustrates the defining property of the bias-reduced estimation principle: as one varies the nuisance parameter values allowed by the chosen working models, the squared asymptotic bias of the doubly robust estimator is locally minimized at (2.254,−0.959)≈(αBR∗,γBR∗) in certain directions of γ and α. Bias-reduction is indeed local, as we observe that there are two other regions where even smaller bias near zero is attained. The location of these regions depends strongly on the underlying data generating mechanism P0, so that they cannot be obtained by solving an estimating equation without further knowledge of the true data-generating mechanism.□

Figure 1:

Contour plot of the log of the squared asymptotic bias log{bias2(α,γ;1)} as a function of the nuisance parameters α and γ for example 1 with ×=(1.115,1.606)≈(αMLE∗,γMLE∗), ∙=(2.254,−0.959)≈(αBR∗,γBR∗) and ▲▲=(1.115,0.169)≈(αMLE∗,γMLE−BR∗).

Example 2

To illustrate that the location of the aforementioned regions depends on the underlying data generating mechanism, we show a similar plot for a second example. With n=105, for each individual i, we let Xi=dN(1,1), Ri|Xi=dBer{π0(Xi)} with π0(Xi)=expit(−1+Xi2) and Yi|Xi=dN(Xi2,1). Under this data generating mechanism, the true mean outcome equals ψ0=2.

In Figure 2, the cross “×” now indicates the point (αˆnMLE,γˆnMLE)=(0.737,2.157), which approximates (αMLE∗,γMLE∗), leading to ψˆnDR,MLE=2.393 and bias(αˆnMLE,γˆnMLE; 2)=0.393. The bullet “∙” indicates the point (αˆnBR,γˆnBR)=(0.609,1.220), which approximates the probability limits (αBR∗,γBR∗), leading to ψˆnBR=2.316 and a smaller bias given by bias(αˆnBR,γˆnBR;2)=0.316. The behavior in this second example (see Figure 2) is similar to the behavior in the first example (see Figure 1). However, do note that the regions where smaller bias near zero is attained are located in a different position, not trackable by solving an estimating equation without further knowledge of the true data-generating mechanism.□

Figure 2:

Contour plot of the log of the squared asymptotic bias log{bias2(α,γ;2)} as a function of the nuisance parameters α and γ for example 2 with ×=(0.737,2.157)≈(αMLE∗,γMLE∗), ∙=(0.609,1.220)≈(αBR∗,γBR∗) and ▲▲=(0.737,0.859)≈(αMLE∗,γMLE−BR∗).

As previously noted, it follows from eq. (2) that small perturbations (of the order one over root-n) in the nuisance parameters γ and α do not affect the first-order asymptotic behavior of the doubly robust estimator when the intersection model m(qγ)∩m(gα) holds. Besides targeting bias reduction, the bias-reduced estimation procedure extends this property to working model misspecification. This is visually suggested by both Figures 1 and 2 which show that the estimators αˆnBR and γˆnBR are situated in a region where small perturbations of the nuisance parameters do not heavily affect the doubly robust estimator (in contrast to the regions that deliver near-zero bias, which are very unstable). This is unlike the MLE of the nuisance parameters. More formally, the doubly robust estimator ψˆnBR is first-order ancillary [28] under misspecification of the working models [8], Corollary 1). This means that even under possible model misspecification, ψˆnBR still admits the expansion ψˆnBR−ψ0=(P0−Pn)[D∗{Qˉ(γBR∗),g(αBR∗); ψ0}]+op(1/n) and thus ψˆnBR and ψˆn{Qˉ(γBR∗),g(αBR∗)} are asymptotically equivalent. This moreover implies that the asymptotic variance of the doubly robust estimator ψˆnBR can be straightforwardly estimated as one over n times the sample variance of the values D∗{Qˉ(ˆγnBR),g(αˆnBR);ψˆnBR}, without having to adjust for the estimation of the nuisance parameters.

A limitation of the bias-reduced estimation procedure is that it is confined to parametric working models of the same dimension. Equal dimensions are indeed required because the gradient of D∗ with respect to α is used as an estimating function for γ and vice versa. This can be overcome by targeting asymptotic bias reduction in the direction of a single nuisance parameter, for instance in the direction of α. The effect of this is best understood from Figures 1 or 2. For a given value α˜ of α, e.g., the MLE αˆnMLE, and for each value of γ, we evaluate how the bias of the doubly robust estimator changes as we move away from the chosen value α˜. We then choose the outcome regression parameter γ at which no change in the squared asymptotic bias is seen. For example 1, with MLE αˆnMLE=1.115, this leads to γˆnMLE−BR=0.169, ψˆnMLE−BR=0.518 and bias(αˆnMLE,γˆnMLE−BR;1)=−0.482. On Figure 1, this point (αˆnMLE,γˆnMLE−BR)=(1.115,0.169) is marked by means of the triangle “▲”. This is the point where no change of bias is seen in the direction of α. Note that the bias lies between the bias of the MLE and the bias-reduced doubly robust estimator. A qualitatively different result is seen for example 2. With MLE αˆnMLE=0.737, we find γˆnMLE−BR=0.859, ψˆnMLE−BR=2.302 and bias(αˆnMLE,γˆnMLE−BR; 2)=0.302, which is slightly smaller than the bias of the bias-reduced doubly robust estimator. This is an artefact of the true underlying data-generating mechanism P0; we are not aware of theoretical reasons for this bias to be smaller (see example 1 where it is higher). In the next section, we will explain how we can use this idea of targeting asymptotic bias reduction in one direction to be able to make use of data-adaptive learning algorithms.

4.1 Main idea

Our proposal is to start with a parametric model for the missingness mechanism g0(X), indexed by α, and an estimator αˆn of α, e.g., the MLE. Next, an estimator for the conditional mean outcome Qˉ0(X) is obtained, either by maximum likelihood estimation under a parametric model or by using data-adaptive learning algorithms such as super-learning [10]. With concern for misspecification of the missingness model, we next fluctuate this initial estimator of Qˉ0(X) through a parametric fluctuation model, indexed by a finite-dimensional parameter ε of at least the dimension of α. This model includes covariates that are chosen in such a way that the score of the fluctuation parameter ε equals the gradient Dα∗{Qˉ(γ),g(α)}. By doing so, we aim to reduce the asymptotic bias of the doubly robust estimator in the direction of α, as explained in Section 3.2. In Section 4.2, we will detail how this can be done. Note that we can thus allow for a missingness model of arbitrary dimension by exploiting the bias-reduced estimation principle in only a single dimension.

The idea of extending an initial fit of the conditional mean outcome is not new: this idea was already proposed in Bang and Robins [6] and it also underlies the TMLE procedure of van der Laan and Rubin [9], as well as other improved doubly robust estimation procedures [20, 21]. Our proposal is nonetheless different in that it explicitly targets bias reduction (in the direction of α). In contrast, the TMLE procedure aims at obtaining a doubly robust substitution estimator and the procedures by Tan [20] and Rotnitzky et al. [21] target a bounded doubly robust estimator with desirable efficiency properties.

4.2 Practical implementation of the procedure

The extension of the bias-reduced doubly robust estimation procedure can be implemented using the following four steps.

Step 1: Estimatorg^nMLEfor the missingness mechanismg0(X).

Postulate a parametric working model for g0(X): gα={g(α)(X)|α∈Rp} where g(α)(X)=G{αTl(X)} and G is an appropriate inverse link function and l=(1,l1,…,lp−1), e.g., the logistic regression model G(⋅)=expit(⋅). Obtain the MLE αˆnMLE solving

with DgMLE(α) given by eq. (6). Let gˆnMLE=g(αˆnMLE).

Step 2: Initial estimatorQ¯^n0for the conditional mean outcomeQ¯0.

The second step of the procedure is to obtain an initial estimator for Qˉ0. We describe two possibilities: (1) a parametric model and (2) a super-learner.

1. Parametric Working Model. The first option is to postulate a parametric working model for Qˉ0(X): qγ={Qˉ(γ)(X)|γ∈Rq} with Qˉ(γ)(X)=Q{γTk(X)}, Q an appropriate inverse link function and k=(1,k1,…,kq−1); e.g., a linear regression model with Q the identity link. Obtain the MLE ˆγnMLE solving

with DQˉMLE(γ) given by eq. (5). Let Q¯^nMLE=Q¯(γ^nMLE).

1. Super-Learner. Another option is to obtain an initial estimator based on data-adaptive learning algorithms, such as super-learning [10]. The super-learner is a machine learning algorithm which starts from a library of estimators {Q¯^j|j=1…,J}, which may consist of both parametric and nonparametric estimators. It then considers the family of all weighted averages of these estimators: Q¯^ω=∑j=1JωjQ¯^j, ω=(ω1,…,ωJ)T, ∑j=1Jωj=1 and ωj≥0 for j=1,…,J. Next, the optimal weight vector ωˆn is defined to be the choice of ω that minimizes the cross-validated risk with respect to some loss function lSL:(O,Qˉ)→lSL(Qˉ)(O)∈R which satisfies Qˉ0=argminQˉE{lSL(Qˉ)}, e.g., the squared error loss function l2(Qˉ)(O)=R{Y−Qˉ(X)}2. The super-learner of Qˉ0 is defined as the estimator Q¯^nSL=Q¯^ω^n.

   For further reference, we let the initial estimator Q¯^n0 denote either Q¯^nMLE or Q¯^nSL.

Step 3: FluctuationQ¯^n(c)of the initial estimatorQ¯^n0.

To construct an appropriate fluctuation model, we need to choose an appropriate loss-function. In this article, we consider the quasi-log-likelihood loss function with corresponding logistic fluctuation model. This choice is favored over the squared error loss function with linear fluctuation model [30] because it ensures predictions within the admissible range for the outcome. In particular, suppose Y is known to fall in the interval [a,b] with a<b. To be able to use the quasi-log-likelihood loss functions and the logistic fluctuation model, we rescale the outcome between zero and one as follows:

The procedure we describe below is based on the transformed outcome Y˜ and results in an estimator ψ˜n of ψ˜n=E(Y˜). Because E(Y)=(b−a)E(Y˜)+a, the final estimator is given by ψˆn=(b−a)ψ˜n+a. For notational convenience, without loss of generality, we will assume that a=0 and b=1 so that Y=Y˜∈[0,1] and we can drop the ∼-notation.

Having thus obtained an estimator for the missingness mechanism, such as gˆnMLE, and an initial estimator Q¯^n0 for the conditional mean outcome taking values in the interval [0,1], we can now fluctuate Q¯^n0 to target bias reduction in the direction of the finite-dimensional parameter α. Below, we consider three fluctuation models that will accomplish this goal with the second also guaranteeing the final estimator to be a substitution estimator, like the TMLE procedure.

1. Fluctuation model 1. Consider the fluctuation model {Q¯^n0(ε(1)):ε(1)∈ℝp} through the initial estimator (Q¯^n0(0)=Q¯^n0):

   (12)logit Q¯^n0(ε(1))(X)=logit Q¯^n0(X)+ε(1),TH(1)(g^nMLE)(X),

   with H(1)(g^nMLE)(X)=[G'{α^nMLE,Tl(X)}/G2{α^nMLE,Tl(X)}]l(X) and with logit(x)=log{x/(1−x)}. Then define εˆn(1) as the MLE under this parametric submodel, i.e., such that

   ε^n(1)=arg minε(1)Pn[ℒ(1){Q¯^n0(ε(1))}]

   where l(1) is the quasi-log-likelihood loss function;

   (13)l(1)(Qˉ)(O)=−R[Ylog{Qˉ(X)}+(1−Y)log{1−Qˉ(X)}].

   It is easily verified that εˆn(1) solves the estimating equation

   (14)∑i=1nRi{Yi−Q¯^n0(ε^n(1))(Xi)}H(1)(g^nMLE)(Xi)=0,

   which can be solved via standard logistic regression of the observed outcomes on the covariates H(1) using offset logit Q¯^n0. Because the score eq. (14) is like the estimation eq. (9), it targets bias-reduction (in the direction of α) in the sense that was previously described. Define the updated estimator Q¯^n(1)=Q¯^n0(ε^n(1)). Because the quasi-log-likelihood loss function is a valid loss-function for the conditional mean outcome Qˉ0 in the sense that Qˉ0=argminQˉE{l(1)(Qˉ)} (see Gruber and van der Laan [30], Lemma 1), Q¯^n(1) is an improved fit of Qˉ0 as compared to Q¯^n0 with respect to the loss function (13).

2. Fluctuation model 2. We extend eq. (12) so that the final doubly robust estimator of the target parameter will also be a substitution estimator (also known as a regression doubly robust estimator [12]). For this purpose, define the fluctuation model {Q¯^n0(ε(2)):ε(2)∈ℝp+1} through the initial estimator:

   logit Q¯^n0(ε(2))(X)

   (15)=logit Q¯^n0(X)+ε(2,1),TH(1)(g^nMLE)(X)+ε(2,2)H(2)(g^nMLE)(X),

   ε(2)=(ε(2,1),T,ε(2,2))T and H(2)(gˆnMLE)=1/gˆnMLE. Then define εˆn(2) as the MLE under this parametric submodel, i.e., such that

   ε^n(2)=arg minε(2)Pn[ℒ(1){Q¯^n0(ε(2))}]

   where l(1) is the quasi-log-likelihood loss function (13). It follows that εˆn(2) solves the estimating equation

   (16)∑i=1nRi{Yi−Q¯^n0(ε^n(2))(Xi)}×(H(1),T(g^nMLE)(Xi),H(2)(g^nMLE)(Xi))T=0,

   which can be easily solved via standard logistic regression of the observed outcomes on the covariates H(1) and H(2) using offset logit Q¯^n0. By adding H(2) to the fluctuation model, it follows from the second component of eq. (16) that the resulting doubly robust estimator is a substitution estimator (as in the TMLE procedure) because eq. (16) implies that ∑i=1nRi/g^nMLE{Yi−Q¯^n0(ε^n(2))(Xi)}=0 so that ψ^n{Q¯^n0(ε^n(2)),g^nMLE}=Pn{Q¯^n0(ε^n(2))}. Define the updated estimator Q¯^n(2)=Q¯^n0(ε^n(2)).

3. Fluctuation model 3. There is evidence in the literature that the use of a weighted loss function in the construction of a fluctuation model can improve the stability of the doubly robust estimator of interest; see for example Robins et al. [12], Díaz and Rosenblum [31] for a comparison of an unweighted versus a weighted loss function in the context of TMLE. We therefore propose the following alternative fluctuation model to eq. (12): {Q¯^n0(ε(3)):ε(3)∈ℝp} through the initial estimator:

   (17)logit Q¯^n0(ε(3))(X)=logit Q¯^n0(X)+ε(3),Tl(X).

   Then define εˆn(3) as the weighted MLE under this parametric submodel

   ε^n(3)=arg minε(3)Pn[G′{α^nMLE,Tl(X)}G2{α^nMLE,Tl(X)}ℒ(1){Q¯^n0(ε(3))}],

   with l(1) the quasi-log-likelihood loss function (13); thus εˆn(3) solves the estimating equations

   (18)0=∑i=1nRi{Yi−Q¯^n0(ε^n(2))(Xi)}G′{α^nMLE,Tl(Xi)}G2{α^nMLE,Tl(Xi)}l(Xi),

   which can be easily solved via standard weighted logistic regression of the observed outcomes on the covariates l using offset logit Q¯^n0. As before this targets bias reduction because eq. (18) is like the estimating eq. (9).Define the updated estimator Q¯^n(3)=Q¯^n0(ε^n(3)).