Improved Doubly Robust Estimation in Marginal Mean Models for Dynamic Regimes

2.1 Notation and Assumptions

Without loss of generality, we follow much the same notation given in [11]. Let A_j, j = 1, . . . , K, be the stochastic treatment decision at time t_j and Ā_j = (A₁, . . . , A_j) be the history of treatment decisions through time t_j. Similarly, define tailoring and auxiliary variables, S_j and X_j, respectively, for j = 0, . . . , K − 1, as well as their histories S̅_j = (S₀, . . . , S_j) and X̅_j = (X₀, . . . , X_j). Let 𝒢 denote subgroups of interest. The set of all possible potential outcomes is {Y(a¯K) | a¯K∈A¯K} , where A¯K is the collection of all possible treatment allocation vectors; similarly, the set of potential intermediate outcomes is {S¯K(a¯K) | a¯K∈A¯K} . Hence, we assume that the potential outcome does not depend on the treatment assignment mechanism nor is affected by others’ treatment [14, 18]. We also assume a sequential randomization assumption [14], that is, the treatment assignment A_j is conditionally independent of {𝒢, S₀, . . . , S_K(ā_K−1), Y(ā_K−1)} given {L₀, A₁, L₁, . . . , A_j−1, L_j−1}, with L_j = (S_j, X_j). In words, the sequential randomization implies (i) that at no point in time does treatment assignment depend on future outcomes, and (ii) that the history of auxiliary variables X̅_j−1 is sufficiently rich to ensure that treatment assignment A_j does not vary systematically by intermediate or potential outcomes within levels of (Ā_j−1, L̅_j−1). The observed data is O = {𝒢, L₀, A₁, L₁, A₂, . . . , L_K−1, A_K, Y}, where Y = Y(Ā_K), S_j = S_j(Ā_j), X_j = X_j(Ā_j) for all j = 1, . . . , K.

Let d_j(L̅_j−1) be a treatment decision or treatment assignment rule at time t_j, j = 1, . . . , K, and the sequence of decision rules for the entire treatment allocation period defines the dynamic treatment regimen, i.e. d̅_K = (d₁, . . . , d_K). Now, define the product

where, at time t_m, π_m0(· | Ā_m−1, L̅_m−1) is the treatment assignment probability in the observational data and treatment assignment in the dynamic regime is degenerate according to the treatment decision rule d_m(L̅_m−1). Without loss of generality, we define π¯j=(π1,…,πj) and the vector of true treatment assignment probabilities π¯j0=(π10,…,πj0) for j = 1, . . . , K. The product in (1) plays a critical role in our ability to estimate causal estimands from the observed data and is a version of the non-stabilized inverse PS weight. Let P_d̅K and E_d̅K be probability distribution and expectation under the dynamic regime, respectively, and P and E be the probability distribution and expectation in the observed data, respectively. Then, if there is non-trivial probability that a randomly selected subject can follow the regime d̅_K in the observational, then by Lemma 4.1 of [11], the distribution of (Y, S̅_K−1, Ā_K, 𝒢) under P_d̅K is absolutely continuous with respect to the distribution of (Y, S̅_K−1, Ā_K, 𝒢) under P and a version of the Radon-Nikodym derivative is

2.2 Ordinary doubly robust estimation

Suppose that we are interested in estimating the parameter vector β from the marginal model,

where μ(β, 𝒢) is a parameterization of the conditional mean of Y given subgroups 𝒢 under the dynamic regime d̅_K. Then, under regularity assumptions outlined in Section 2.1 as well as in Lemma 4.1 of [11], the statistic

where μ˙(β,G)=(∂/∂β)μ(β,G) , has mean zero at the true marginal parameter β = β₀. If the treatment assignment probabilities {π_j0} were known a priori, this statistic could be used to define an estimator for β. However, in observational studies, the treatment assignment probabilities are unknown a priori and instead one posits a statistical model. To this end, we model the treatment assignment probabilities through a finite-dimensional vector of parameters γ and estimate it through the estimating function,

where uj(aj | A¯j−1,L¯j−1;γ)=π˙j(aj | A¯j−1,L¯j−1;γ)Vj−1(aj | A¯j−1,L¯j−1;γ){aj−πj(aj | A¯j−1,L¯j−1;γ)} , π˙j(aj | A¯j−1,L¯j−1;γ)=(∂/∂γ)πj(aj | A¯j−1,L¯j−1;γ) and 𝒱_j(a_j | Ā_j−1, L̅_j−1; γ) is the conditional variance of a_j given Ā_j−1 and L̅_j−1. The estimating function U_γ(O; γ) may be regarded as proportional to the first-order partial derivative of the log likelihood with respect to γ for the treatment assignment mechanism. As such, EU_γ(O; γ₀) = 0 when the treatment assignment probabilities are correctly modeled. Let γ^ denote the maximum likelihood estimator of γ₀, i.e. the solution to the estimating equations, 0 = ℙ_nU_γ(O; γ). Then, the inverse probability weighted estimator (IPW) for β is β^IPW , the solution to the system of estimating equations,

where

A concern of IPW estimators, assuming {π_j} are correctly modeled, is their efficiency. Briefly, the likelihood may be written as the product of two expressions, ℒ(O) = ℒ_OR(O)ℒ_PS(O), where ℒPS(O)=∏j=1Kπj0(Aj | A¯j−1,L¯j−1) and ℒ_OR(O) is the product of conditional densities, ∏j=1KP(Lj*∈dLj | A¯j−1,L¯j−1) , where Y ≡ L_K and P(Lj*∈B | ⋅)=P(Lj∈B | ⋅) for every measureable rectangle B. The causal estimand E_d̅K (Y | 𝒢) is a function of the conditional densities in ℒ_OR but not {π_j0} in ℒ_PS; hence, the treatment selection probabilities are nuisance parameters for the estimand of interest. The efficient estimator for β will be orthogonal to the score function for the nuisance parameters and, hence, the IPW estimator can be improved by subtracting from ψ_IPW its projection onto the score function of the treatment selection probabilities. The projection of ψ_IPW onto the score function of the treatment selection probabilities [15] is

where the regression functions are nested within one another through the following relationship: g_K0(Ā_K, L̅_K−1) = E(Y | Ā_K, L̅_K−1) and for j = 1, . . . , K − 1,

Again, we adopt the notation g̅_j = (g₁, . . . , g_j) and the history of true regression models is g̅_j0 = (g₁₀, . . . , g_j0). In semi-parametric theory for missing data problems [22, 26], the augmentation term A(O;β,π¯K0,g¯K0) in (5) has mean zero by construction when evaluated at β = β₀, and the treatment assignment probabilities are modeled correctly so that {π_j0} ≡ {π_j(γ₀)}.

In order to use the augmentation in practice, the true regression functions {g_j0} must be known or modeled. Often, g_j0(Ā, L̅_j−1) is modeled parametrically or semi-parametrically as g_j(Ā_j, L̅_j−1; α) through the finite dimensional parameter α. For example, Murphy et al. [2001] suggested the semi-parametric estimator α^MLR , the solution to the estimating equations, 0 = ℙ_n[U_α(O; α)],

where g˙j(A¯j,L¯j−1;α)=(∂/∂α)gj(A¯j,L¯j−1;α) for j = 1, . . . , K. Using α^MLR defined through (7), then the usual augmented inverse probability weighted (AIPW) estimator, β^USUAL , is defined as the solution to the system of equations

where

and {g_j(α)} = {g_j(Ā_j, L̅_j−1; α)}. The AIPW estimator has the double robustness property: that is, β^USUAL is consistent if one of either the treatment selection probabilities {π_j} or the regression functions {g_j} are modeled correctly. If both functions are modeled correctly, the AIPW estimator is semi-parametric efficient [16, 22].

3.1 New Estimation with {π_j} as Known Functions

In this subsection, we consider constrained doubly robust estimation of β when {π_j} are known functions of (Ā_j−1, L̅_j−1), i.e., when γ = γ₀ is known. To begin, consider the coefficient estimator β˜ defined as the solution to

where α^* is the limiting value of the estimator α^ regardless of whether {g_j} are modeled correctly or not. The asymptotic variance of n1/2(β˜−β0) is

where Γ is the asymptotic slope matrix of the right-hand side of (8) with respect to β₀, say,

The asymptotic slope matrix is invariant to the choice of α_* either directly or in the limit. Therefore, when the {π_j0} are modeled correctly, minimizing the asymptotic variance V˜(α∗) reduces to a problem of minimizing var [ψAIPW(O;β0,π¯K0,g¯K(α*))] as a function in α_*. [For a related problem, see expression (4) on p.8 of [23]].

Now, because ψAIPW(O;β0,π¯K0,g¯K(α*)) has mean zero when evaluated at the true parameters, then minimizing var [ψAIPW(O;β0,π¯K0,g¯K(α*))] is equivalent to minimizing

where α_* is the probabilistic limit of some estimator α^ . Define α_opt as the minimizer of (9) as a function in α_*. Note that α_opt also satisfies the following system of equations,

where

When {g_j} are correctly specified, α_opt = α₀; otherwise, α_opt is some other value in the parameter space. To satisfy (P2), we wish to define an estimator α^opt→pαopt while, at the same time, α^opt→pαopt≡α0 when {g_j} are correctly specified, so that (P1) is satisfied.

To proceed with construction of the improved DR estimator via constrained augmentation, it is convenient to re-express (10). Using Lemmas 1–3 in Appendix A, we show that the estimating function on the right-hand side of (10) can be written as

We can then use (11) to propose an estimator for α. We propose to estimate α by solving the system of estimating equations,

where, for r = 1, . . . , K,

In Appendix A, we show that the estimating function on the right-hand side of (12) has mean zero when one of either {π_j} or {g_j} are correctly specified and, hence, the double robustness property (P1) is satisfied. Furthermore, when the PS model is correctly specified and (12) is evaluated with {π_j0}, we show that (12) has mean zero at α = α_opt, and therefore, the constrained augmentation improved doubly robust (IDR) estimator β^IDR has the smallest asymptotic variance among the class of AIPW estimators with correctly specified {π_j}, where β^IDR is the solution to the AIPW estimating equations with α^ solving (12). Thus, both properties (P1) and (P2) are satisfied.

3.2 New Estimation with {π_j} as Functions of Unknown γ

We now consider the case where the probabilities {π_j} are functions of an unknown parameter vector γ. In general, var [ψAIPW(O;β,π¯K(γ0),g¯K(α)] is not equal to var [ψAIPW(O;β,π¯K(γ^),g¯K(α)] in either finite or large samples. Exceptions include when {π_j} are known by design or when the estimators γ^ are constructed to be asymptotically orthogonal. The details given in this subsection are appropriate for a majority of cases where there is a non-negligible effect on var [ψAIPW(O;β,π¯K(γ0),g¯K(α)] when γ₀ is replaced with a consistent estimator γ^ . When the variance cost to the estimating function ψ_AIPW is zero or asymptotically negligible, methods described in Section 3.1 are appropriate exactly or in the limit.

For j = 1, . . . , K, let π_j(γ) = π_j(a_j | Ā_j−1, L̅_j−1; γ) and π¯j(γ)={π1(γ),…,πj(γ)} . Similar to the outline in Section 3.1, we aim to define an estimator β^ that solves the AIPW estimating equations,

where γ^ and α^ are estimators of γ and α, respectively, such that the double robustness property in (P1) and optimality property in (P2) are both satisfied. From Theorem 9.1 in [22], when the PS model is correctly specified, the AIPW estimating equations in (14) are asymptotically equivalent to

where α_* is the limiting value of some estimator α^ and

Let c_* be the value of c that minimizes

when α_* is substituted for α. Our objective is thus to find ζ_opt that minimizes (16) where ζ = (a, c), and the procedure is analogous to finding α_opt that solves (10) in Section 3.1. However, in the current scenario, we allow for the fact that {π_j0} are functions of the unknown parameter γ that is estimated from the data. We can then use this correspondence to Section 3.1 to propose an estimator of ζ_opt, say ζ^opt , that will in turn lead to the improved doubly robust (IDR) estimator β^IDR by solving (14) with smallest asymptotic variance when the PS models are correctly specified but the OR models may or may not be; see Appendix B for a detailed derivation of this section.

Building on the principles and logic laid out above, we therefore propose to estimate ζ by solving the following estimating equation,

where

and

and g˜jα(dj(L¯j−1),A¯j−1,L¯j−1;ζ,γ)=g˙j(dj(L¯j−1),A¯j−1,L¯j−1;α) , and u_j(a_j | Ā_j−1, L̅_j−1; γ) is defined in (3). We observe that γ^ converges to γ₀ when the PS models {π_j} are correctly specified regardless of whether the OR models {g_j} are correctly specified or not. Using similar arguments in Appendix A that show the results in Section 3.1, ζ^opt solving the system of estimating equations in (17), will converges to ζ_opt and thus satisfying (P2). On the other hand, when the PS models are misspecified but the OR models are correct, the expectation of (17) can be shown to be zero when ζ ≡ (α, c) = (α₀, 0) and therefore the double-robustness property (P1) is satisfied. Thus, our estimator β^IDR , defined as the solution to 0=ψAIPW(O;β0,π¯K(γ^),g¯K(α^opt)) where α^opt is defined through ζ^opt via (17) and γ^ is an M-estimator defined in (3), is DR and has smallest asymptotic variance among all DR estimators when the PS models are correct but OR models may not be.

We may use a theory of Z-estimation to describe the asymptotic properties of β^IDR . Specifically, we combine the estimating equations for θ = (β^⊤, γ^⊤, α^⊤)^⊤ and derive the asymptotic properties for β^ from the asymptotic properties of θ^=(β^⊤,γ^⊤,α^⊤)⊤ . Combining the estimating functions for β and γ with the proposed estimating function (17) with respect to ζ = (α, c), θ^ solves the system of estimating equations,

Under conditions used to derive the improved DR estimator [e.g., 11, Appendix], θ^ is a Z-estimator with an asymptotic normal distribution whose covariance follows standard formulae [e.g.3, Ch. 7] and can be estimated in various ways, including direct evaluation of the empirical covariance matrix and resampling methods.

3.3 Caveats

The technique described above is general and applicable for all non-random dynamic treatment regimes via the marginal model in (2). This development applies to the work by [4] and [23] who detailed the idea for a population mean in the presence of missing outcome and monotone coarsening, respectively.

Although the principles behind locally efficient doubly robust estimation is technically sound, putting the principles into practice can be a substantial challenge. Firstly, for K-stage trials with K even moderately large and non-trivial treatment effect, there are many regression functions to specify. Such bookkeeping can be an onerous task and encourages the use of a subclass of {g_j} that restrict regression coefficients to be common across stages. Second, the regression functions {g_j} are highly constrained through their nested definition in (6). This caveat led [11] to propose that {g_j} depend on auxiliary variables X̅_j−1 but not on treatment Ā_j−1, i.e. no treatment effect. Unfortunately, if there is a treatment effect, then one is relying entirely on a correct PS model via {π_j} for consistency of the AIPW estimator. As we demonstrated above, even if the PS model is correct, incorrectly specified {g_j} could have undesirable consequences on the precision of the AIPW estimator and downstream statistical inference. Thus, efficiency of the AIPW estimator in the current context is both germane and important.

4.1 The Setup: Adaptive Treatment Length Studies

Our numerical and real-data examples are illustrations of our continued work in dynamic regimes for infusion studies, where a clinical goal of interest is to estimate the mean outcome if an attending physician stops the infusion at a time of their choice or at the time of an unexpected random event (i.e. a time not of their choice), whichever time comes first [8]. In these applications, as with many dynamic treatment regimes, evaluating eligibility criteria for treatment assignment is a key feature of the analysis. Eligibility for treatment assignment is evaluated sequentially at each stage or clinic visit at time t_j and depends on treatment history, covariate history, and the history of eligibility. In the infusion study, treatment assignment at t_j is synonymous with a physician's decision to stop or continue treatment by choice at time t_j provided the patient satisfies eligibility criteria. In the observational study, treatment assignment depends on treatment history, covariate history, and the history of eligibility whereas treatment assignment depends on treatment and eligibility history in the non-random dynamic regime for our infusion study application.

To define the decision rule leading to our dynamic regime, let a target treatment length t_r be given, r ∈ {1, . . . , K}. Let S_j be a binary tailoring variable. We consider the decision rule

In the infusion study application, a subject initiates infusion treatment at t₀. According to the decision rule in (18), at t_j, j < r, the attending physician continues treatment as long as the subject satisfies the eligibility criteria, i.e. S̅_j = 0, where S_j is a random (in)eligibility event in the interval (t_j, t_j+1]. If Y is the clinical endpoint of interest in the infusion study, then our goal is to estimate the causal estimand, E_d̅K (Y | 𝒢) = β, using the data.

The observed data for the proposed application is outlined in Table 1. In contrast to a dynamic regime whose treatment length is stopped at t_r provided S̅_r−1 = 0, subjects may stop at any t_j, j < K, in the observational study. In addition, the concept of eligibility for treatment assignment at t_j is broader in the observational study. In order to be eligible for treatment assignment at t_j, a subject must have not stopped treatment at any prior time point t_j′, j′ < j, as well as avoided random adverse events {S_j′ = 1} for all j′ < j. In our notation, A_j = 1 implies a decision to stop or “assign” treatment at t_j by choice and the probability that a provider stops treatment at t_j given a subject is eligible at t_j is

Then, the treatment assignment mechanism is modeled as

Adhering to the non-stochastic dynamic regime that defines the adaptive treatment length policy is given by the degenerate treatment assignment mechanism,

Combining the degenerate treatment assignment mechanism above with the mechanism in (19) from the observed data leads to the definition of the Radon-Nikodym derivative, i.e.,

that is critical for estimating the causal estimand E_d̅K [Y | 𝒢] from observational data.

4.2 Simulation Studies

We performed simulation studies to evaluate the finite sample performance of competing estimators of for the causal parameter β in the marginal model μ(β, 𝒢), including outcome regression (OR), inverse probability weighted (IPW), ordinary double robust and improved double robust estimators. First, we partition covariate history as X¯j=((V0(1))⊤,(V¯j(2))⊤)⊤ , where V0(1) baseline (time-invariant) covariates, and V¯j(2)=(V0(2),V1(2),…,Vj(2)) is the history of time-varying covariate information, respectively. Next, treatment assignment probabilities are generated as Bernoulli random variables with success probability λj(X¯j−1;γ)=H(γj+γX⊤X¯j−1) , j = 1, . . . , K, where H(z) = 1/{1 + exp(−z)} and γX=(γV(1)⊤,γV(2)⊤)⊤ . Intermediate outcomes S_j are similarly simulated as Bernoulli random variables with success probability λjξ(X¯j;α0)=P(Sj=1 | X¯j,A¯j−1=0,S¯j−1=0) , j = 0, . . . , K − 1, where λjξ(X¯j;α)=H(αξ,j+αξX⊤X¯j) and αξX=(αξ,V(1)⊤,αξ,V(2)⊤)⊤ . Finally, the endpoint Y is assumed to follow a linear model, Y = g_K0(Ā_K, L̅_K−1; α) + €, where g_K0(Ā_K, L̅_K−1; α₀) = E(Y | Ā_K, L̅_K−1),