Handling missing data when estimating causal effects with targeted maximum likelihood estimation

Abstract Targeted maximum likelihood estimation (TMLE) is increasingly used for doubly robust causal inference, but how missing data should be handled when using TMLE with data-adaptive approaches is unclear. Based on data (1992-1998) from the Victorian Adolescent Health Cohort Study, we conducted a simulation study to evaluate 8 missing-data methods in this context: complete-case analysis, extended TMLE incorporating an outcome-missingness model, the missing covariate missing indicator method, and 5 multiple imputation (MI) approaches using parametric or machine-learning models. We considered 6 scenarios that varied in terms of exposure/outcome generation models (presence of confounder-confounder interactions) and missingness mechanisms (whether outcome influenced missingness in other variables and presence of interaction/nonlinear terms in missingness models). Complete-case analysis and extended TMLE had small biases when outcome did not influence missingness in other variables. Parametric MI without interactions had large bias when exposure/outcome generation models included interactions. Parametric MI including interactions performed best in bias and variance reduction across all settings, except when missingness models included a nonlinear term. When choosing a method for handling missing data in the context of TMLE, researchers must consider the missingness mechanism and, for MI, compatibility with the analysis method. In many settings, a parametric MI approach that incorporates interactions and nonlinearities is expected to perform well.


Introduction
A key component of epidemiological research is causal inference from longitudinal studies, where the objective is often to estimate the average causal effect (ACE) of an exposure on an outcome.(1)(2)(3)(4)(5) For a binary exposure (X=1 exposed; X=0 unexposed), the ACE can be defined as the difference in the average potential outcome if all were exposed versus unexposed.(1)(2)(3)(4)(5) In the absence of missing data, under the assumptions of conditional exchangeability given a vector of measured confounders Z, consistency, and positivity, the ACE is identifiable from observable data by the g-formula E[E(Y|X = 1, Z) − E(Y|X = 0, Z)], where Y is the outcome.(6) Several singly robust approaches, including g-computation and propensity score methods, and doubly robust estimators, including Targeted Maximum Likelihood Estimation (TMLE), are available for ACE estimation in the absence of missing data.Here, we focus on TMLE, which combines models for the outcome and propensity score.(7-9)A detailed description of TMLE is available elsewhere.(7)(8)(9) Briefly, the first step is the same as in g-computation, where a model for the expected outcome conditional on exposure and confounders (E / [Y|X, Z]) is estimated and used to predict outcome for all records under exposure and no exposure.In g-computation these predictions are directly plugged into the g-formula to estimate the ACE.In TMLE, they are updated using information from the propensity score (P / [X = 1|Z]) (the targeting step), before being plugged into the g-formula.(7) The targeting step in TMLE ensures that the estimator is doubly robust, where only one of the two models (outcome or propensity model) needs to be consistently estimated to ensure consistent estimation of the ACE.Unlike with singly robust methods, data-adaptive approaches (e.g. machine learning methods) can be used for the exposure and outcome models in TMLE as long as the Donsker class condition (requiring that outcome and propensity score estimators do not heavily overfit the data) holds.(8)(9)(10) Given these desirable properties, interest in the application of TMLE for ACE estimation is growing.
Missing data are ubiquitous in epidemiological studies and can lead to biased estimates and loss of precision if handled inappropriately.(11) Commonly used approaches for handling missing data in studies using TMLE for ACE estimation include multiple imputation (MI; e.g., (12)), complete-case analysis (CCA; e.g., (13)), extending TMLE to handle missing outcome data (e.g., (14)), and the missing covariate missing indicator (MCMI) approach for handling missing confounder data (e.g., (15)).No study has compared these methods in terms of bias and precision.Also, while a requirement for valid inference with MI is that the imputation model should incorporate all relationships assumed to hold in the analysis method, (11) the optimal implementation of MI when using TMLE with data-adaptive methods for ACE estimation is unknown.Answers to these questions are key to developing guidance for appropriate handling of missing data in the context of the growing use of TMLE in applied epidemiological research.
In the current paper, we sought to address this knowledge gap using a simulation study based on an illustrative example from the Victorian Adolescent Health Cohort Study (VAHCS).Our interest was to compare the performance of readily available missing data methods to inform current practice.We begin by introducing the VAHCS example, then describe methods for handling missing data with TMLE, present the simulation study we conducted to evaluate and compare the performance of these approaches, then illustrate the assessed approaches in the VAHCS example.We conclude with a general discussion.

Illustrative example
Our example was based on a previous investigation using data from VAHCS, a longitudinal cohort study of 1,943 participants (1,000 females), recruited in 1992-1993 at ages14-15 years.(16) Data were collected during participants' adolescence (waves one-six) and in young adulthood (1998, wave seven).Investigators aimed to estimate the ACE of frequent cannabis use in adolescent females on mental health in young adulthood, measured using the Clinical Interview Schedule-Revised (CIS-R).We revisited this analysis using TMLE with data-adaptive approaches, using the same confounders considered in the original investigation: parental divorce, antisocial behaviour, depression and anxiety, alcohol use, parental education, all measured across waves two to six.(17) Table 1 shows descriptive statistics for analysis variables, and age at wave two, which is useful as an auxiliary variable in MI (a predictor of missing values but not included in the analysis method (11)).All variables had some degree of missingness, ranging from 0.1% to 30.8%.

Methods for handling missing data in TMLE
Approaches proposed for handling missing data when estimating the ACE using TMLE in studies like our example are described below.

Non-MI approaches
Complete-case analysis (CCA): only records with complete data for target analysis variables are used.(18) This approach generally leads to loss of precision (19) and, depending on missingness mechanism, may inflict bias.(18) Extended TMLE in a sample with complete exposure and confounders (Ext-TMLE+CEC): records with missing data for Z and X are deleted.From records with complete data on Y, E[Y|X, Z] is estimated.In the targeting step, the predictions of the outcome are updated using information from models for P[X=1|Z] and P[M !=0|X,Z], where M ! is the missingness indicator for the outcome.(20) Updated predictions for the outcome under exposure and no exposure are obtained for all records, regardless of whether they have missing outcome.(20) The model for M !can also be fitted using data-adaptive approaches.With missingness only in the outcome, the extended TMLE method is unbiased under an extended exchangeability assumption (Y " ∐ M !|X, Z and Y " ∐ X|Z for x = 0,1, where Y " is the potential outcome when X = x, and ∐ denotes independence).(21) Extended TMLE plus missing covariate missing indicator (MCMI) approach (Ext-TMLE+MCMI): missing outcome data are handled using the extended TMLE approach, and missing confounder data by including missingness indicators for the incomplete confounders in the confounding adjustment set.Records with missing exposure data are excluded.In settings with complete exposure and outcome data, the MCMI approach can be expected to yield an unbiased ACE estimate under an extended exchangeability assumption (Y " ∐ X|Z, M # forx = 0,1, where M # is the vector of missingness indicators for the incomplete confounders), and the assumption that the exposure or outcome depend on the confounder only when the confounder is observed.(22,23) This assumption is plausible in some settings, such as electronic health record data where the decision to prescribe a medication might be influenced by family history of disease only when the clinician has the relevant information.

MI approaches
The MI fully conditional specification (FCS) framework (24) enables simultaneous handling of missing exposure, confounder, and outcome data.Under this approach, univariate imputation models are specified for each incomplete variable conditional on other variables and imputations are drawn sequentially until convergence.(24) The process is repeated multiple times to generate multiple completed datasets.Analysis is performed within each completed dataset and the results are pooled using Rubin's rules to obtain the final estimate and its standard error (SE).(24) For valid inference with MI, each univariate imputation model should be tailored to be compatible with the analysis method.To achieve this, all analysis variables and complexities such as interaction terms in the target analysis should be included as predictors in each univariate imputation model.(24) There are various possible implementations of MI within the FCS framework: Parametric MI with no interaction (MI-no int): each univariate imputation model is based on a regression model with main-effects terms only -the default in most MI software.In the example and simulations, we considered main-effects logistic regression for the binary variables and predictive mean matching (PMM) for the continuous outcome based on a main-effects linear regression.In PMM, imputed values are drawn using the nearest observed value after fitting the regression (24) which makes it robust to misspecification of the latter, for example in the presence of nonlinear associations.(25) Parametric MI with two-way interactions (MI-2-way int): each univariate imputation model is based on a regression model as above, but including two-way exposure-outcome, exposure-confounder, confounderoutcome, and confounder-confounder interactions.Interaction terms are generated within each cycle of the MI algorithm from current values of relevant variables involved in the interaction term (the so-called "passive" approach in mice in R). (26) Parametric MI with two-, three-, and four-way interactions (MI-higher int): as above, but all univariate imputation models additionally include three-and four-way confounder-confounder interactions as predictors.
MI using classification and regression trees (MI-CART): instead of regression, for each variable with missing data a tree is fitted using a recursive partitioning technique, with all other variables as predictors.
Each record belongs to a donor leaf, from which a randomly selected value for the variable is taken as the imputed value.(27) MI-CART (and MI using random forest (MI-RF), see below) have been proposed to enable imputation that can more flexibly allow for interactions and non-linearities.(27) MI using random forest (MI-RF): for each variable with missing data, multiple bootstrap samples are drawn and for each a separate tree is fitted.Each tree contributes a donor leaf, and a randomly selected value for the variable from all these donors is taken as the imputed value.(27) All MI approaches can be implemented with the mice package in R. (26)

Simulation study
To compare the performance of the described methods for handling missing data, we performed a simulation study based on the VAHCS example (Figure 1).We considered six scenarios.For each, we generated 2,000 datasets of 2,000 records.

Generating the complete data
We used parametric regression models to generate the variables.The values of parameters in the models were determined by fitting similar models to the available data in VAHCS (unless stated otherwise).We considered two data generating scenarios (simple and complex), differing in the confounder-confounder interaction terms used in the data generation models.No exposure-confounder interaction terms were included (no effect modification).Supplementary Table 1 provides the parameter values used for simulating the data and Supplementary Table 2 the descriptive statistics of the variables in the simulated data.
For both scenarios, we generated a continuous auxiliary variable A(age at wave 2) and a set of confounders Z = (Z $ (parental divorce), Z % (antisocial behaviour), Z & (depression and anxiety), Z ' (alcohol use), Z ( (parental education)).The models for generating these variables were (all binary variables coded 0/1 and The scenarios differed in the exposure and outcome generation models. Simple scenario: we used main-effects regression models to generate a binary exposure X(frequent cannabis use) and a continuous outcome Y(log-transformed, standardized CIS-R total score): Complex scenario: we used regression models that included confounder-confounder interactions (excluding interactions with Z % because of the low prevalence (15%)) to generate the exposure and outcome: We inflated the coefficients for the interaction terms, approximately four times larger than values estimated in the VAHCS data.Under both outcome generation models, we set the coefficients for X (θ $ ), i.e., the true value of the ACE, to 0.2.For this effect and 2,000 records, the null hypothesis of no causal effect is formally rejected (p < 0.05) in approximately 80% of the simulated datasets.

Imposing missing data
We considered missingness scenarios depicted by m-DAGs (missingness directed acyclic graphs) A and B (Figure 2).(19, 28).These were selected as they represent scenarios under which our expectations of the performance of methods were quite distinct (see Discussion).
We imposed missingness on Z % , Z & , Z ' , X, Y through generating missingness indicators , M !, coded 0 if the variable was observed; 1 if missing.We considered variables A, Z $ , Z ( , which had <10% missing within VAHCS (Table 1), as fully observed in the simulation study.
For the simple scenario, the models used for generating the missingness indicators according to each m-DAG were: For each missingness indicator, we set the coefficients for all confounders and the exposure to 0.9.We set the coefficient for the outcome to 0 for m-DAG A and to 0.1 for m-DAG B.
For the complex scenario, we considered two sets of models for generating missingness according to each m-DAG.The first set, hereafter complex scenario 1, used the same models as for the simple scenario, detailed above.The second set, hereafter complex scenario 2, used the models: For the complex scenario 2, we set the coefficients for all confounders, the exposure, and the outcome the same values as for the simple scenario, detailed above.We set the coefficients for exposure-confounder interactions (XZ % , XZ & , XZ ' ) to 0.9, and for Y % to 0 for m-DAG A and 0.08 for m-DAG B.
This led to overall 6 scenarios (2 m-DAGs × 3 data generating scenarios).The missingness proportions were the same in all missingness scenarios and approximately the same as in the real VAHCS dataset, except for the outcome, which was increased to 20% (13% in VAHCS).
The overlap between missingness proportions were adjusted so that the proportions of records excluded were 50% for CCA, 40% for Ext-TMLE+CEC, and 30% for Ext-TMLE+MCMI (Supplementary Table 2).

Analysis of the simulated data
The target analysis aimed to estimate the ACE of X on Y using TMLE with data-adaptive methods adjusting for Z $ , Z % , Z & , Z ' , Z ( as confounders.We used the TMLE package in R (20).We fitted the exposure and outcome models using a SuperLearner library that included a range of parametric, semiparametric, and nonparametric methods.(29,30) The SE for TMLE was obtained using the variance of the influence function.(9) The analysis was applied to each simulated incomplete dataset alongside each of the previously described missing data methods.For Ext-TMLE+CEC and Ext-TMLE+MCMI we used the same SuperLearner library for the outcome missingness model.We used the mice package in R to implement MI. (26) Due to computational constraints, for each MI approach, we generated five imputed datasets (see Discussion).( 24) Supplementary Tables 3 and 4 show the variables and interaction terms included in each imputation model for MI-2-way int and MI-higher int.We used the default settings of the mice package for the donor pool for PMM in parametric MI approaches and for the hyperparameters for MI-CART and MI-RF.(26)

Evaluation criteria
We compared the performance of the approaches for handling missing data by calculating the relative bias percent, the empirical standard errors (SE), and the percent error in average model-based SE relative to the empirical SE.For all, Monte-Carlo SEs (MC-SE) were obtained.(31) Analyses were performed in R version 3.6.1.(32)
Biases in complex scenario 1, for both m-DAGs, were similar to those in the simple scenario, except that MI-no int had larger bias than in the simple scenario (33% m-DAG A, 23% m-DAG B).
In the complex scenario 2, for m-DAG A, biases for non-MI methods (<|15%|), MI-CART (-24%) and MI-RF (-58%) were larger than in prior scenarios, while the three parametric MI approaches performed similarly to the complex scenario 1.For m-DAG B, all the non-MI and parametric MI approaches and MI-RF had larger bias than in prior scenarios (-42% to -43% for non-MI approaches, 35% to 43% for parametric MI approaches, -64% for MI-RF), while MI-CART had similar bias as in the simple scenario and complex scenario 1.
The MC-SE for relative bias ranged from 0.81%-1.98%across scenarios and m-DAGs.

Empirical standard error and relative error in model-based standard error
Across scenarios and m-DAGs, the empirical SEs using CCA and Ext-TMLE+CEC were similar (0.13-0.16) and larger than Ext-TMLE+MCMI (0.12-0.14) (Figure 4).The SEs obtained from the parametric MI approaches were similar to each other (0.10-0.14), to Ext-TMLE+MCMI, and MI-CART, except for complex scenario 2 under m-DAG B, where the parametric MI approaches exhibited larger SEs (0.16-0.18).The empirical SEs were similar across scenarios and m-DAGs for MI-CART (0.10-0.13) and MI-RF (0.07-0.09).
MI-RF had the lowest empirical SE in all scenarios and m-DAGs (0.7-0.8).The MC-SE for empirical SEs ranged from 0.001-0.003.
The model SEs were underestimated using non-MI methods and overestimated using MI methods across all scenarios and m-DAGs (Figure 5).The errors were smallest under the simple scenario and largest under the complex scenario 2. Within each scenario, the performance among non-MI approaches was similar.The performance of the MI approaches was similar within each scenario, except MI-RF, which produced model SEs with considerably larger error.The MC-SE for relative % error in model-based SEs ranged from 1.15%-3.02%.

Illustrative example results
We analysed the VAHCS example using the tmle package in R, (20) and applied the eight missing data methods described.Unlike in the simulations, a small proportion of participants had missing data for parental divorce and parental education (Table 1), which were handled in the same way as missing data for the other confounders.Also, the auxiliary variable age had 9.3% missing data, which was multiply imputed in the MI approaches.For the MI approaches, 100 imputations were performed.
The obtained effect sizes were small, with MI-no int yielding somewhat larger effect size and non-MI methods and MI-RF smaller effect sizes (Table 2).The SEs for MI approaches were larger than the non-MI methods, which could be explained by the downward and upward biases in model SEs for non-MI and MI approaches, respectively, observed in our simulation study (Figure 5).For example, using the relative percent error in model SEs averaged over the six scenarios in the simulations, the corrected SEs in the case study would be 0.13 for CCA, 0.13 for MI-no int and 0.11 for MI-RF.
In the VAHCS example, the outcome (mental health in young adulthood) might well have influenced its own missingness, in which case neither of the considered m-DAGs in the simulation study are plausible for our example and we would expect all the considered missing data methods to be biased.(19)

Discussion
We compared methods for handling missing data when estimating the ACE using TMLE with data-adaptive approaches.We considered six scenarios with different exposure and outcome generation models (presence/absence of confounder-confounder interaction terms) and missingness mechanisms (whether the outcome influenced missingness in other variables and presence/absence of interaction/non-linear terms in missingness models).CCA and EXT-TMLE+CEC had small bias under m-DAG A (where the outcome did not influence missingness in other variables), and large bias otherwise.MI-no int -the default in most MI software -had large bias in the complex scenarios 1 and 2 (when the exposure/outcome generation models included interactions) and small bias otherwise (regardless of m-DAG).MI-2-way int and MI-higher int performed best in terms of bias and variance across all settings, except for m-DAG B in complex scenario 2 (where a non-linear outcome term influenced missingness in other variables).MI-RF had consistently large bias across the six scenarios.MI-CART had small bias under m-DAG A in the simple and complex scenario 1, and large bias otherwise.
Based on previous investigations in a setting without effect modification, we determined that for m-DAG A, where the outcome did not influence missingness in any variable, the ACE was identifiable (or "recoverable").(19) Further, because auxiliary variables did not influence missingness in any variable, we expected both CCA and an appropriate implementation of MI to yield low bias (33).Indeed, CCA (and Ext-TMLE+CEC) produced estimates with small bias for m-DAG A across all data generation scenarios.Also, implementations of parametric MI that were approximately compatible with the analysis method (i.e., all parametric MI procedures in the simple scenario, and MI including interaction terms in complex scenarios) returned estimates with little bias, while an inappropriate MI method (e.g., MI-no int in our complex scenarios) was considerably more biased.
Contrary to CCA and Ext-TMLE+CEC, the Ext-TMLE+MCMI approach had higher bias under m-DAG A.
A key assumption under which the MCMI approach has been shown to be unbiased is when the exposure or outcome only depend on the confounder when the confounder is observed.(22,23) We did not consider missingness scenarios where this held, because this assumption is implausible in a prospective cohort study, like VAHCS, where the data are not used for medical decision-making.
For m-DAG B, the ACE was determined to be non-recoverable, but since the outcome did not influence its own missingness, based on a previous simulation study (19), we speculated that an implementation of MI that was tailored to the analysis method may offer some bias reduction compared with CCA.We observed that parametric MI approaches including interactions performed better than all other approaches in terms of bias for the simple scenario and complex scenario 1.For the complex scenario 2, where the missingness models included exposure-confounder interactions and a quadratic term for the outcome, all parametric MI approaches were highly biased.MI-CART outperformed parametric MI approaches in this scenario, but it was still moderately biased.
Within the FCS framework, recursive partitioning techniques, such as CART and RF, have been suggested as alternative approaches that could automatically incorporate interactions and non-linearities in the imputation process.(27) Previous simulation studies have shown that MI using CART performs better than parametric MI without interaction terms.(27,34) However, in these studies, the target analysis was a correctly specified outcome regression model with interactions, and biases in estimates of the main effects were not that different following MI using CART or parametric MI without interaction.These studies imposed missingness in the outcome only (27) or outcome and covariates.(34) In both, missingness depended on fully observed variables.In the present study, the only setting where MI-CART outperformed parametric MI approaches was for m-DAG B in the complex scenario 2. In all scenarios, bias in the ACE approaches, in terms of bias of point estimates as discussed previously, but even more so for bias in variance estimates.A promising alternative approach for obtaining SEs for MI in the presence of incompatibility has been recently proposed using the bootstrap,(38) but we did not explore this because of computational constraints.Incompatibility between the models used in MI and TMLE may also compromise the double robustness property of TMLE.(39) Other approaches, such as EXT-TMLE considered in this paper to handle missing outcome data, and an alternative likelihood parameterization to construct doubly robust estimators in adjusting for confounding and missing data in the presence of missing confounder data,(39) guarantee double robustness under extended assumptions.
Our simulation study was broadly based on VAHCS.We evaluated the performance of missing data methods under various missingness mechanisms.To describe what variables influenced missingness, we used m-DAGs because the standard classification of missing completely at random (MCAR), missing at random (MAR), or missing not at random (MNAR) is difficult to comprehend and substantively assess when there is missingness in multiple variables.Also, although it is possible to estimate key parameters unbiasedly if the MAR assumption holds, MAR is not necessary for unbiased estimation.(19,33) We did not consider missingness mechanisms where outcome influenced its own missingness, under which none of the approaches could be expected to perform well.For each MI approach, due to computational constraints we generated five completed datasets in the simulation study, which is fewer than we would do in practice.(24) We do not expect this to have affected the comparison between MI approaches, but it could have affected comparison of non-MI with MI methods.Our simulated data had a relatively simple structure across the assessed scenarios.Extensions of our study could investigate the performance of these missing data methods for datasets with high-dimensional confounders, binary outcome, and more complex m-DAGs including longitudinal auxiliary variables.

Conclusion
We evaluated the performance of eight available approaches to handle missing data when estimating the ACE using TMLE with data-adaptive approaches under various data generation scenarios and missingness mechanisms.Our results highlight the importance of considering the missingness mechanism and compatibility with the analysis method when choosing a method to handle missing data.In many settings, a parametric MI approach that incorporates interactions and non-linearities is expected to perform well in the context of TMLE with data-adaptive approaches.effect estimated as the difference in the mean potential outcome under exposure and under no exposure Abbreviations -Ext-TML: Eextended targeted maximum likelihood estimation (TMLE) approach; Ext-TMLE+MCMI: extended TMLE plus missing covariate missing indicator (MCMI) approach; MI, no int (linear): parametric multiple imputation (MI) with no interaction -linear regression to impute missing outcome; MI no int: parametric MI with no interaction -predictive mean matching to impute missing outcome; MI, 2-way int: parametric MI with two-way interactions; MI higher int: parametric MI with two-, three-, and four-way interactions; MI, CART: MI using classification and regression trees; MI, RF: MI using random forest  (19) For simplicity of exposition, confounders without missing data (Z ! and Z " ) are presented on a single node and confounders with missing data (Z # , Z $ , Z % ) on another single node.Also, only one missingness indicator has been included for confounders with missing data (M & ), coded as 1 when any of the variables Z # , Z $ , Z % have missing data and as 0 when none has missing data.Auxiliary variable (continuous) Confounders (binary) Exposure (binary) Outcome (continuous) Figure 3 -Relative bias percent (%) in ACE estimation (filled circles with error bars showing ±Monte Carlo standard errors) using different missing data methods for simple and complex scenarios and m-DAGs A and B. For all missing data methods, TMLE was implemented using the SuperLearner including the following methods: mean (the average), glm (generalized linear model), glm.interaction (generalized linear model with 2-way interactions between all pairs of variables), bayesglm (Bayesian generalized linear model), gam (generalized additive model), glmnet (elastic net regression), earth (multivariate adaptive regression splines), rpart (recursive partitioning and regression trees), rpartPrune (recursive partitioning with pruning) and ranger (random forest).Abbreviations -Ext-TML: extended targeted maximum likelihood estimation (TMLE) approach; Ext-TMLE+MCMI: extended TMLE plus missing covariate missing indicator (MCMI) approach; MI-no int (linear): parametric multiple imputation (MI) with no interaction -linear regression to impute missing outcome; MI-no int: parametric MI with no interactionpredictive mean matching to impute missing outcome; MI-2-way int: parametric MI with two-way interactions; MI-higher int: parametric MI with two-, three-, and four-way interactions; MI-CART: MI using classification and regression trees; MI-RF: MI using random forest.Relative % error in model SE

FiguresFigure 1 -Figure 2 -
FiguresFigure 1 -Directed acyclic graph (DAG) used in data generation for the simulation study; abbreviations: CIS-R revised clinical interview schedule

Figure 4 -
Figure4 -Empirical standard error in ACE estimation (filled circles with error bars showing ±Monte Carlo standard errors) using different missing data methods for simple and complex scenarios and m-DAGs A and B. Abbreviations -Ext-TML: Extended targeted maximum likelihood estimation (TMLE) approach; Ext-TMLE+MCMI: extended TMLE plus missing covariate missing indicator (MCMI) approach; MIno int (linear): parametric multiple imputation (MI) with no interaction -linear regression to impute missing outcome; MI-no int: parametric MI with no interaction -predictive mean matching to impute missing outcome; MI-2-way int: parametric MI with two-way interactions; MI-higher int: parametric MI with two-, three-, and four-way interactions; MI-CART: MI using classification and regression trees; MI-RF: MI using random forest.
(36)mates following MI-RF was larger than MI-CART, consistent with Dove et al.'s results.(27)Wespeculatethatthismighthavebeenbecause in the implementation of MI-RF used a small number of randomly preselected predictor variables to split at each node, which might have negatively impacted prediction accuracy.(35)Specifically,therangerpackageusedwithinR'smice package for implementing RF uses, as the default number of preselected variables to split at each node, the square root of the total number of variables in the imputation model, which was 2 in our simulation.(27)CART,ontheotherhand, does not involve variable preselection.(35)Inthisstudy,non-MI approaches underestimated the model SE, which was not surprising.If both the exposure and outcome models are consistently estimated and the Donsker class condition is satisfied, TMLE is an asymptotically linear estimator and its variance can be obtained based on the variance of the influence curve.(9,10)Itis, however, unclear if the Donsker class condition is met when data-adaptive approaches are used for the exposure and outcome models.(36)Thisleads to bias in variance estimation as has been (38)36,37)ere and in other simulation studies.(30,36,37)It is an ongoing area of research to develop approaches to tackle it, such as cross-validated TMLE, which allows asymptotic linearity to be established without the Donsker class condition.(10)Additionally,Rubin'sMIvariance estimator is expected to perform poorly in the presence of incompatibility,(38)which might explain the overestimation of model SEs for the MI approaches.Incompatibility is the key challenge for using MI with TMLE with data-adaptive

Table 1 -
Description of variables used in the case study, their distribution and the proportion with missing data among the VAHCS female participants (n=1,000)

Table 2 -
Estimated average causal effect (ACE) of frequent cannabis use during adolescence on CIS-R score (standardised z-score) using a TMLE approach under different missing data methods within the VAHCS case study Abbreviations -Ext-TML: Extended targeted maximum likelihood estimation (TMLE) approach; Ext-TMLE+MCMI: extended TMLE plus missing covariate missing indicator (MCMI) approach; MIno int (linear): parametric multiple imputation (MI) with no interaction -linear regression to impute missing outcome; MI-no int: parametric MI with no interaction -predictive mean matching to impute missing outcome; MI-2-way int: parametric MI with two-way interactions; MI-higher int: parametric MI with two-, three-, and four-way interactions; MI-CART: MI using classification and regression trees; MI-RF: MI using random forest.Relative % error in model standard error in ACE estimation (filled circles with error bars showing ±Monte Carlo standard errors) using different missing data methods for simple and complex scenarios and m-DAGs A and B. Abbreviations -Ext-TML: Extended targeted maximum likelihood estimation (TMLE) approach; Ext-TMLE+MCMI: extended TMLE plus missing covariate missing indicator (MCMI) approach; MI-no int (linear): parametric multiple imputation (MI) with no interaction -linear regression to impute missing outcome; MI-no int: parametric MI with no interaction -predictive mean matching to impute missing outcome; MI-2-way int: parametric MI with two-way interactions; MI-higher int: parametric MI with two-, three-, and four-way interactions; MI-CART: MI using classification and regression trees; MI-RF: MI using random forest.

Table 6 -
Simulation study results for complete data for simple, intermediate, and complex scenario, using outcome regression, g-computation including all two-way confounder-confounder interactions, and TMLE

Table 7 -
Performance of missing data methods for the 11 assessed missingness mechanisms under the simple scenario

Table 8 -
Performance of missing data methods for the 11 assessed missingness mechanisms under the complex scenario 1

Table 9 -
Performance of missing data methods for the 11 assessed missingness mechanisms under the complex scenario 2