{'fraction_non_alphanumeric': 0.03684541077249753, 'fraction_numerical': 0.008792654843436911, 'mean_word_length': 4.583235932542996, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A coherent likelihood parametrization for doubly robust estimation of a causal effect with missing confounders

Missing data and confounding are two problems researchers face in observational studies for comparative effectiveness. Williamson et al. (2012) recently proposed a unified approach to handle both issues concurrently using a multiply-robust (MR) methodology under the assumption that confounders are missing at random. Their approach considers a union of models in which any submodel has a parametric component while the remaining models are unrestricted. We show that while their estimating function is MR in theory, the possibility for multiply robust inference is complicated by the fact that parametric models for different components of the union model are not variation independent and therefore the MR property is unlikely to hold in practice. To address this, we propose an alternative transparent parametrization of the likelihood function, which makes explicit the model dependencies between various nuisance functions needed to evaluate the MR efficient score. The proposed method is genuinely doubly-robust (DR) in that it is consistent and asymptotic normal if one of two sets of modeling assumptions holds. We evaluate the performance and doubly robust property of the DR method via a simulation study.


Introduction
Confounding bias and missing data are two major analytic challenges in comparative effectiveness research using observational data such as electronic medical records. While each problem has been thoroughly studied separately, consolidated approaches for addressing both issues are lacking. In the absence of missing data, confounding bias must still be adjusted for in order to evaluate causal effects (Hernan et al., 2004). Researchers often use the g-formula for identifying the distribution of counterfactuals from the observed data distribution (Robins, 1986;Snowden et al., 2011). Inverse probability weighting estimators are commonly used and involve modeling the propensity score (Rosenbaum and Rubin, 1983;Robins, 1986;Hernan et al., 2000). Doubly-robust estimators for causal effects have been well established and widely studied (Bang and Robins, 2005;Lunceford and Davidian, 2004;Robins, 2000;Vansteelandt et al., 2007). These estimators are doubly robust in the sense that they are consistent and asymptotic normal if either the treatment mechanism (propensity score) or outcome model is correctly specified, but not necessarily both. These methods are also locally semiparametric efficient because they achieve the semiparametric efficiency bound for the nonparametric model when model misspecification is absent, that is, at the intersection submodel of the union of specified models.
Multiple imputation and inverse probability of censoring weighting are increasingly popular methods for addressing missing data (Rubin, 2004;Li et al., 2013;Seaman and White, 2013). In the context of regression analysis, various weighting schemes to account for missing covariates have been examined previously in the literature (Moore et al., 2009;Lipsitz et al., 1999;Parzen et al., 2002;Tchetgen Tchetgen, 2009). Semiparametric locally efficient methods are also available to address data missing at random i.e., the probability of observing the full data depends on the fully observed data only (Kang and Schafer, 2007;Bang and Robins, 2005). Robins, Rotnitzky and others examined improved augmented inverse weighted estimators within the semiparametric framework (Robins et al., 1994;Robins and Rotnitzky, 1995;Rotnitzky and Robins, 1997;Scharfstein et al., 1999). Additionally, Tsiatis (2007) provides an extensive overview of the state of the art for applying semiparametric theory to missing data.
However, to date, few methods have considered joint inferences about causal effects that are doubly or multiply robust in the presence of missing data and confounding. This setting presents a special challenge in that it involves the nesting of causal inference in the missing data setting, each of which requires, to obtain the parameter of interest, estimating a nuisance parameter while appropriately accounting for the fact that nuisance parameters needed to adjust for selection bias are entangled with nuisance parameters needed to address confounding bias. Entangled in the sense that we now need to account for modeling both the confounder and the missingness of that confounder, which are both typically a nuisance in their own right. Davidian et al. (2005) presented a doubly-robust augmented inverse weighted estimator of the causal effect of exposure when the outcome was missing. In their 2003 textbook, Robins and van der Laan give a unified theory for addressing causal inference in the presence of missing data but do not address specific challenges with identifying an appropriate parametrization for the observed data when addressing both confounding adjustment and incomplete confounder data (van der Laan and Robins, 2003). This paper addresses a special case of that general theory. Williamson et al. (2012) attempt to combine existing methods in order to create a multiply-robust estimator. The authors consider a union of four semiparametric models each of which specifies parametric working models for either the missingness mechanism or the missing covariates to account for missing data, and for either the treatment mechanism or the outcome to account for confounding. Multiplyrobust estimation requires that each submodel of the union model is a semiparametric model in that the correctly specified part is parametric, while the remaining submodels are unrestricted. However, we will show in this paper that the rest of the likelihood is, in fact, restricted in at least one submodel of the union model and therefore the multiply-robust property claimed by Williamson et al. (2012) may not be achievable in reality. An immediate implication of this phenomenon is that in addition to possible lack of compatibility across submodels of the union model, the intersection submodel of the union model may in fact be empty. Therefore, unless one explicitly acknowledges the overlap between components of the union model in the process of model specification, one may in fact rule out the possibility of achieving local efficiency.
In this paper we discuss the difficulty of achieving double robustness in semiparametric missing data when full data nuisance parameters are entangled with nuisance parameters needed to account for data missing at random. We carefully examine the previously suggested multiply-robust method and explain why it may fail to achieve the claimed multiply-robust property. We then propose a solution that carefully identifies the modeling assumptions through an alternative transparent parametrization of the likelihood function, which makes explicit model dependencies between various nuisance functions needed to evaluate the multiply-robust estimating equation for the causal effect of interest. The proposed method is genuinely doubly-robust in that it is consistent and asymptotically normal if one of two sets of modeling assumptions holds. Further, due to the inherent model dependencies, we establish that double-robustness to model misspecification is the best one hope to achieve in this setting. This paper suggests an approach that could easily be adopted in other settings where one may wish to obtain a doubly-robust estimator in the presence of entangled nuisance parameters. While the paper focuses on the effect of treatment on the treated, the proposed approach equally applies to the average causal effect.

Full Data Setting
Let A denote a binary treatment, A ∈ {0, 1}, and let Y be the outcome in view with Y 1 and Y 0 denoting the potential outcomes under treatment and control conditions respectively. Let W denote a set of pretreatment covariates. The parameter of interest is the effect of treatment on the treated on the additive Throughout, we make the following standard causal assumptions in order to identify the effect of treatment on the treated: Assumption 3 Positivity: pr(A=0|W ) pr(A=1|W ) > 0 almost surely.
Assumption 1 states that a person's observed outcome corresponds to her potential outcome for the observed treatment. Assumption 2 states that the treatment assignment is ignorable conditional on covariates W , i.e. W includes all common causes of A, Y 1 and Y 0 . And assumption 3 states that there is no treated subject without an untreated counterpart.
Under assumption 1, θ = E[Y | A = 1] (Angrist and Pischke; Kennedy et al., 2015). Under assumptions 1-3, Ψ is well known to be non-parametrically identified and is The following dual representation of (1) is also of interest where µ is a dominating measure of the distribution of (W, Y ).
Any regular and asymptotically linear estimatorΨ of Ψ satisfies is a zero-mean function, called the i th influence function for Ψ. The influence function characterizes the behavior of the estimator (such as the asymptotic distribution) and under certain condition, may also be used to define an estimating equation to obtain an estimator with the corresponding influence function. For functionals defined on nonparametric models, as will be considered in this paper, there exists a unique influence function that is semiparametric efficient under the nonparametric model. Influence functions were first introduced by Huber (Huber, 1972) in the context of robust statistics and later developed in semiparametric theory, in the sense of Bickel et al. (1993). In either context, influence functions represent the influence of a single observation on the estimator.
The efficient influence function of Ψ in the nonparametric model in which assumptions 1-3 hold, but the form of the observed data likelihood is unrestricted is (Hahn, 1998), In order to use the efficient influence function as an estimating function for Ψ when, as is typically the case in observational studies, W is high dimensional, one must estimate the nuisance functions

Missing Data Setting
Next, suppose that only a subset of covariates, C, of W is fully observed while L is missing for a subset of participants where W = {L, C}. Therefore the observed data can be written as (Y, A, C, RL, R), where R is an indicator function which is equal to 1 when L is observed and is otherwise equal to 0.
Define O = (Y, A, C), the fully observed data. Furthermore, suppose that L is missing at random. The data now presents a non-monotone missingness pattern with respect to missing confounders and missing counterfactual outcomes.
In order to address missing data, we make the following additional assumptions: Assumption 4 π = pr(R = 1 | A, L, C, Y ) > 0 almost surely; Assumption 5 Conditional exchangeability: L ⊥ R | A, C, Y .
Assumption 4 is a positivity assumption and states that there is a positive probability of observing any possible value of (A, C, L, Y ) in the complete cases. Assumption 5 is a missing at random assumption and states that the conditional distribution of L given A, C, Y is the same in incomplete and complete cases.
Under assumptions 1-5, the efficient influence function of Ψ in the nonparametric model in which the observed data distribution is unrestricted is The efficient influence function in equation (3) depends on the following functions: the propensity score, p = pr(A = 1 | L, C), the outcome model, m = m (Y | A, L, C), the missing data mechanism, π = pr(R = 1 | A, Y, C), and the density of L given A, C and Y , t = t(L | A, C, Y ).
The efficient influence function is appealing as a basis for obtaining inferences about Ψ, mainly because of the following multiple robust property: if Ψ is evaluated at the truth, and one of the following statements hold: (i) p and π are evaluated at the truth; (ii) m and π are evaluated at the truth; (iii) p and t are evaluated at the truth; (iv) m and t are evaluated at the truth.
Additionally, at the intersection submodel where all of the models are evaluated at the truth, the variance of ι M iss (Ψ) achieves the semiparametric efficiency bound for the union of models (i)-(iv) at the intersection submodel.
A closely related multiply-robust property of the efficient influence function to account for missing confounders was established for the average causal effect by Williamson et al. (2012).
However, because in practice one must estimate p, m, π, and t under corresponding low dimensional working models, model incompatibility may render the multiply-robust property given above infeasible, as we show next. The approach considers four submodels, p, m, π, and t and four unions of those submodels. These submodels are semiparametric in the sense that in practice, within each submodel of the union model, two models are parametrically specified, but the remaining two are not modeled and left unrestricted. However, this cannot hold as the specified models in at least one of the submodels of the union model places restrictions on components of the likelihood not explicitly modeled in the submodel. When, as is typically the case in practice, each component of the likelihood is eventually modeled separately, conflict may arise in two separate models for the same component of the likelihood, therefore ruling out the possibility for multiple robustness and local efficiency. To explain how this potential conflict in model specification may arise, consider that the joint likelihood of all four models The submodels t(L | A, C, Y ) and m(Y | A, C, L) are not variation independent, as they both encode an association between Y and L given A and C. Similarly, the models for t(L | A, C, Y ) and p = pr(A = 1 | L, C) are not variation independent because both densities encode the association between A and L given C. Because these various functions are not variation independent, a choice of model for one may restrict modeling options for the other.
As a result of the lack of variation independence, multiply-robust cannot be achieved under a coherent parametrization of the likelihood which acknowledges the model dependence revealed above. Furthermore, unless one such parametrization can be established, local efficiency may also not be attainable because the intersection submodel may be empty in presence of conflicting models. We have provided a straightforward illustration of this phenomenon in the supplementary materials. The following section provides a coherent parametrization of the observed data likelihood under which a certain degree of double robustness can be achieved and local efficiency remains a genuine possibility.

Reparametrization
We propose one possible parametrization of the conditional likelihood function, f (L, Y | A, C), which makes explicit the model dependencies between nuisance functions that are needed to evaluate the efficient score given by (4) in order to obtain an estimator, as described later in section 4. The proposed approach is based on a conditional odds ratio symmetric parametrization of a joint conditional distribution.
Following Chen (2007) and Tchetgen  we define the conditional odds ratio function of A and Y given L as where (a 0 , y 0 ) is a reference value. Chen (2007) established that the joint distribution of A and Y given L can be written as This parametrization is attractive because , and f (Y | a 0 , L) are variation independent in that the choice of a parametric model for one component does not restrict available model choices for another and their joint parameter space is the product space of their respective parameter spaces. We repeatedly make use of the variation independent parameterization result from Chen (2007) in the supplementary materials to prove the following result. Result , and χ (A, L | C) are variation independent parameters. Then, f can be written as and (a 0 , l 0 , y 0 ) are reference values.
This theorem gives a variation independent parameterization of the likelihood and makes explicit that to Similarly we can parameterize the propensity score, p as . This reparametrization makes explicit the fact that the propensity score can be expressed in terms of χ(A, L | C) and f (A | l 0 , C), which are variation independent.
Therefore, both f and p require the same specification of χ(A, L | C), so that they are not variation independent. Furthermore both m and t share the odds ratio of L and Y given C and A. This implies that assuming a submodel for m also places a restriction on the submodel for t which cannot remain unrestricted as assumed by Williamson et al. (2012) in their claim to achieve multiple robustness.

Doubly Robust Inference
Let j = j (L | A = 0, C) denote the distribution of L given A = 0 and C and let j (α) = j (L | A = 0, C; α) be a parametric model for j.
and W ≥ 0 so that W is the true conditional odds ratio function for Y and L given A and C (Chen, 2007) and let w (ω) = w (Y, L | A, C; ω) be a parametric model for w. Let χ (β) = χ (A, L | C, β) be a parametric model for χ(A, L | C) and let π (η) = pr(R = 1 | A, C, Y ; η) be a parametric model for π. Let r = r (Y | A, L = 0, C) be a model for the distribution of Y given A, L = 0, and C and let r (θ) = r (Y | A, L = 0, C; θ) be a parametric model for r. Additionally, let h = pr(A = 1 | L = 0, C) and let h (κ) = pr(A = 1 | L = 0, C; κ) be a parametric model for h.
An estimator α,ω,β,θ of the parameters (α, ω, β, θ), can be found by using direct likelihood maximization of the observed data. This entails maximizing the observed data likelihood, An estimator of η,η, can be found by fitting π (η) to the observed datâ where n is the total number of subjects. Finally, κ can be estimated using inverse probability weighting using 1/π (η) as weights in the complete cases.
An alternative approach using a more standard parametrization can sometimes be used, provided that the parametrization can be shown to satisfy the variation dependence described in Theorem 1, which ensures the existence of a joint distribution for (L, A, Y | C). A standard parametrization in this case implies specifying parametric models for p, t, π, and m needed to evaluate the efficient influence function (3).
Let p (λ) be a parametric model for p, t (φ) be a parametric model for t, π (η) be a parametric model for π, and m (ν) be a parametric model for m. An estimator of η,η , can be found by fitting π (η) on the observed data by using, using, for example, a logistic regression of R on A, C, and Y . We can estimate λ by using inverse probability of censoring weighting using 1 π(η) as weights in the complete cases. For example we might fit a weighted logistic regression of A on C, and L. By specifying m and t as normal with constant variance, one may ensure the existence of a corresponding joint distribution of (Y, L) given A and C, provided the mean of Y given L, A, and C is linear in L and likewise the mean model for L given Y , A, and C is linear in Y . Assumption 5, missing at random, allows us to estimate t (φ) using the complete cases by using, for example, a standard linear regression of L on A, C, and Y . Finally, we describe a simple Monte Carlo algorithm to estimate m (ν): 1. Create M duplicates of the data.
2. Where L is missing, "fill in" the missing variable with a random draw from t φ .
3. Stack all M datasets in long format and estimate m (ν). In practice this will involve fitting a standard model for Y on A, L, and C. For example, a standard main effects linear model.
We then have the following result.
Then under standard regularity conditions,Ψ is consistent and asymptotically normal if the implied form of χ (β) is correctly specified and in addition either (i) π (η) and p λ are consistent for π and p or (ii) m (ν) and t φ are consistent for m and t, but not necessarily both. As before, at the intersection submodel where all of the models are evaluated at the truth, the variance ofΨ achieves the semiparametric efficiency bound for the union of models (i) and (ii).
Then, it is straightforward to show that the solution to equation (3) iŝ The asymptotic distribution of the estimator can be found as follows. Let Q R (η) be an individual contribution to the score for η, Q A λ be an individual contribution to the score for λ, Q L φ be an individual contribution to the score for φ, and Q Y (ν) be an individual contribution to the score for ν.
For example, Also let Z Ψ ,λ,ν,η,φ be an individual contribution to the estimating equation for Ψ. Let Ξ = (η, λ, φ, ν) and define Then, under standard regularity conditions, Therefore, a consistent estimator of the asymptotic variance of √ n(Ψ − Ψ) is with all marginal expectations replaced by their empirical counterparts.
Alternatively, we recommend using the nonparametric bootstrap to obtain estimates of the variance.
In their application to the B-Aware trial, Williamson et al. (2012) use models that are compatible with this new parametrization. As a result, the estimating equation under those choices of models is doublyrobust, though not multiply-robust as they claim. However, their simulation models do not satisfy our proposed parametrization and therefore fail to be compatible, thus there is no chance of it being even doubly-robust. Moreover, the favorable simulation results obtained by the authors can be explained by two reasons. The first reason is that the effect of the missing confounder on the exposure was small compared to the fully observed confounders in the model. The second reason for their favorable simulation results is that the model for the missing confounder is only mildly misspecified. In order to misspecify the model for the missing confounder, the authors omit variables with small regression coefficients and therefore little influence in the model while retaining variables with larger coefficients (Williamson et al., 2012). Further, their complete case estimator performs well which is indicative of a favorable data setting.

Simulation Study
We report a simulation study comparing finite sample performance of our doubly-robust estimator to a number of existing methods. We compared our doubly-robust estimator to an estimator that used Monte Carlo direct likelihood maximization, one using inverse probability of censoring weights, as well as a complete case estimator, a naive estimator that drops the missing confounder, and an estimator calculated from the complete dataset where L was observed for all subjects. The last estimator is obviously not feasible in the presence of missing data however provides a benchmark to assess efficiency loss due to missing data.
The observed data were therefore n = 2, 500 realizations of (R, RL, Y, A, C). Many more details concerning the simulation can be found in the supplementary materials.
In simulations, we implemented the following estimators for comparison: standard inverse probability of censoring weights (IPCW) estimation, Monte Carlo direct likelihood maximization (MCDLM) using 100 imputed datasets, complete-case analysis (CC), and a naive estimator (Naive) that drops the missing confounder L completely and evaluates (1) upon substituting an estimate of pr(A = 1 | C) for p.
For the various methods we fitted the following models. For the missingness mechanism π (η), we fitted a logistic regression. Similarly, we fitted a logistic regression for the propensity score, p (λ), using inverse probability weighting with 1 π(η) as weights in the complete cases. For the distribution of the missing variable, t (φ), we fitted a main effects linear model. Finally, for the outcome model, m (ν), we fitted a main effects linear model.
The inverse probability of censoring weights estimator required π (η) as well as p (λ). The Monte Carlo direct likelihood maximization estimator used t (φ) as well as p (λ). The complete case estimator only required p (λ). The Naive estimator required a logistic regression for A with main effects for C alone, p λ = pr A = 1 | C;λ . All these methods were compared to the proposed doubly-robust estimator which required p (λ), π (η), t (φ), and m (ν).
For the complete-case, naive, Monte Carlo direct likelihood maximization, and inverse probability of censoring weights, we calculated the effect of treatment on the treated for each method using equation (1) for Ψ. For the naive estimator the odds, p (λ) / [1 − p (λ)], were replaced withp λ / 1 −p λ and for the inverse probability of censoring weights estimator the odds were estimated with inverse probability weighting. Our proposed estimator was calculated using equation (5).
The misspecified versions of each model were as follows. The missingness mechanism, π, was misspecified by only using C in the regression, π * = pr(R = 1 | C; λ * ). In order to misspecify a model for p or f we simply (incorrectly) set the coefficient on C to 0 in the working model. This form of misspecification was chosen in order to preserve the structure of the odds ratio between A and L given C, χ(A, L | C), wherever it is required as seen in Section 3.
If L and C are strongly correlated, particularly when the coefficient on C in the propensity score model, λ 2 , is small, then not including C in the propensity score will not be far off from the truth as L will likely suffice to account for confounding. However if L and C are weakly correlated, then any imputation of L that sets the coefficient on C to zero will not be far off from the true model that includes C. Therefore we impose a weak correlation for the simulations exploring misspecification of p and a strong correlation for those misspecifying f as described above. For settings where both are misspecified, we impose a weak correlation L and C. We denote the incorrect propensity score as p * and the incorrect joint distribution of Y and L given A and C as f * . π * used in place of π, (e) π * and p * used in place of π and p (f) f * and p * used in place of f and p, (g) f * and π * used in place of f and π, and (h) f * , π * , and p * used in place of f ,π, and p.
Regardless of model misspecification, the naive and complete-case estimators are biased (a-g) as expected. Similarly, the inverse probability of censoring weights estimator is biased when p * or π * are used in place of p or π (c-h) as it requires both and not a model for f . Additionally, the inverse probability of censoring weights estimator tended to have large variance compared to the other estimators, even under correct specification for p and π. The Monte Carlo direct likelihood maximization estimator is biased when p * is used in place of p(c,e, f , h). When f * is used in place of f , but p is correctly used (b, g) the Monte Carlo direct likelihood maximization estimator is biased, but not overly so. This is likely an artifact of the simulation design regarding the correlation between L and C as explained above. Finally, our doubly-robust estimator is only biased under the settings we expected, namely when f is misspecified along with p or π or both misspecified (f -h). Even in settings where the doubly-robust estimator is biased, the bias is less than that of the other biased estimators. In the setting where π * and f * are used in place of π and f , the bias is comparable to the bias in the Monte Carlo direct likelihood maximization estimator. However, this may be an artifact of the simulation design. Overall, despite a few anomalies, the simulations are in line with expectations. The simulation used in Williamson et al. (2012) did not allow for the range of settings we have explored. Furthermore, their model for R did not include Y . As a result, their missing data mechanism assumption was stronger than missing at random and fairly mild such that their complete case estimates had little or no bias.

Discussion
Analysts are commonly faced with missing data when using observational data to estimate a causal effect. This is particularly true in the setting of two stage non-monotone missingness, such as when potential confounding information is missing along with counterfactual outcomes. The difficulty arises when full data nuisance parameters are entangled with nuisance parameters which are needed to account for data missing at random. In such a setting it is unlikely that researchers will know the underlying mechanisms for the missingness and confounding. Therefore, model misspecification is a likely source of bias when using standard statistical analysis methods. In this paper we have explained why the proposed method of Williamson et al. (2012) fails to achieve the claimed multiply-robust property by carefully examining model dependencies. We identified the modeling assumptions through an alternative parametrization of the joint distribution of the outcome and missing confounder in order to understand  where the red line indicates the truth. DR is our proposed doubly-robust estimator, Naive is the naive estimator that drops a missing confounder, CC is the complete-case estimator, IPCW is the inverse probability of censoring weights estimator, and MCDLM is the Monte Carlo direct likelihood maximazation estimator.
the nuisance parameter entanglements. In this paper we propose a coherent likelihood parametrization and an estimator of the effect of treatment on the treated that accounts for both missingness and potential confounding and that is robust to partial model misspecification.
The simulation study supported the conclusion that our proposed estimator is doubly robust and outperformed existing methods but still failed to be multiply robust as we argued on theoretical basis.
Moreover, we only considered a setting in which a single confounder had missing data. It is more common that several variables may be missing possibly in arbitrary patterns across individuals (Little and Rubin;Robins et al., 1994;Sun and Tchetgen Tchetgen, 2017). Therefore it would be important to extend our approach to allow for arbitrary missing data patterns.

Influence Function Derivation in Full Data Setting
Assumption 1 Consistency: Y = Y A almost surely; Assumption 2 No unmeasured confounding: A ⊥ Y 0 | W , where W = (C, L); Consider a function of the observed data, Our goal is to write dΨt dt as E[ι F ull × S(O)] where ι F ull is the influence function and S(O) is the score function where Then, y). We will drop the subscript for ease of notation. Then, Therefore, the influence function for T 1 is pr(A = 1) . Next Therefore the influence function for T 2 is Now we can see that If one were to use the efficient influence function as an estimating equation for Ψ, then

Proof of Double Robustness for Full Data Setting
If one were to use the efficient influence function, as an estimating equation for Ψ, one would need to estimate the nuisance functions p(A | W ) and m (Y | A = 0, L, C). The resulting estimator,Ψ, is double robust for Ψ in that it will be consistent provided we correctly specify a model for the propensity score, Now let p(A = 1 | W ) = p (W ) for ease of notation and suppose p (W ) is correct while letting

Influence Function Derivation for Missing Data Setting
Recall that Then the influence function for the missing data problem is where π = pr(R = 1 | A, Y, C).
We, in theory, need to correctly specify: From the above expressions for ι M iss , let: Case 1 -p and π are correct Additionally, Case 3 -p and t are correct Additionally, Example: Problem with the multiply robust method The multiply robust method established above requires correct specification of p = pr(A = 1 | L, C) or m (Y | A, L, C) andπ = pr(R = 1 | A, Y, C) or t(L | A, C, Y ). The problem with these model specifications are that they inherently assume we can estimate each model independent of these others and that they are not related. However, they are closely related quantities.
It is possible to use the logit link function to model both t and p: Thus we see that both t and p model the association between L and A, but the former is conditional on Y and C, while the later is only conditional on C. There may not be an intersection submodel for the particular choice of the nuisance models. In simulation we can ensure these models are compatible, but in practice we won't realistically be able to make this assumption.
Similarly we can relate m(Y | A, C, L) and t(L | A, C, Y ). The proposed multiply-robust solution assumes we can specify t(L | A, C, Y ) and m(Y | A, C, L) independently of each other. For example suppose we propose that Unless ν 2 = 0, these models are not compatible in that there does not exist a joint distribution for (L, Y ) with the given families as its conditional distributions. Therefore, for this example, we could never have t and b both be correct and Case 4 above could never be true.

Proof of Theorem 1. Detailed Reparametrization of the Likelihood
We must reparameterize the likelihood because the nuisance parameters overlap.
Lemma 1 We can see that the joint distribution of L and Y can be expressed in terms of f (L | a 0 , C) , As Chen (2007) shows, these are all variation independent parameters. This allows us to estimate f (L, Y | A, C) using maximum likelihood.
Similarly we can write: Thus we see that he propensity score can be expressed in terms of χ (A, L | C) and f (A | l 0 , C) and therefore both the propensity score and joint distribution of L and Y given A and C require correct specification of χ (A, L | C).

Closed Form Estimator
Recall: V 1 requires models for π, p, and for m (0, L, C) = µ 0 Y , which are easily estimated as described in the main body of the paper.