Modeling Treatment Effect Modification in Multidrug-Resistant Tuberculosis in an Individual Patient Data Meta-Analysis
Yan Liu, Mireille Schnitzer, Guanbo Wang, Edward Kennedy, Piret Viiklepp, Mario H. Vargas, Giovanni Sotgiu, Dick Menzies, Andrea Benedetti
MM ODELING T REATMENT E FFECT M ODIFICATION IN M ULTIDRUG -R ESISTANT T UBERCULOSIS IN AN I NDIVIDUAL P ATIENT D ATA M ETA -A NALYSIS
Yan Liu , Mireille Schnitzer , , Guanbo Wang , Edward Kennedy ,Piret Viiklepp , Mario H. Vargas , Giovanni Sotgiu , Dick Menzies , andAndrea Benedetti , , Department of Epidemiology, Biostatistics and Occupational Health,McGill University, Montréal, Québec, Canada Faculty of Pharmacy, Université de Montréal, Montréal, Québec, Canada Department of Social and Preventive Medicine,Université de Montréal, Montréal, Québec, Canada Department of Statistics & Data Science, Carnegie Mellon University, Pittsburgh, USA National Institute for Health Development, Estonia Instituto Nacional de Enfermedades Respiratorias, Mexico City, Mexico. Clinical Epidemiology and Medical Statistics Unit, Department of Medical,Surgical and Experimental Sciences, University of Sassari, Sassari, Italy Respiratory Epidemiology and Clinical Research Unit, McGill University Health Centre,Montréal, Québec, Canada Department of Medicine, McGill University, Montréal, Québec, Canada
January 19, 2021 A BSTRACT
Effect modification occurs while the effect of the treatment is not homogeneous across the differentstrata of patient characteristics. When the effect of treatment may vary from individual to individual,precision medicine can be improved by identifying patient covariates to estimate the size and directionof the effect at the individual level. However, this task is statistically challenging and typically requireslarge amounts of data. Investigators may be interested in using the individual patient data (IPD) frommultiple studies to estimate these treatment effect models. Our data arise from a systematic review ofobservational studies contrasting different treatments for multidrug-resistant tuberculosis (MDR-TB),where multiple antimicrobial agents are taken concurrently to cure the infection. We propose amarginal structural model (MSM) for effect modification by different patient characteristics andco-medications in a meta-analysis of observational IPD. We develop, evaluate, and apply a targetedmaximum likelihood estimator (TMLE) for the doubly robust estimation of the parameters of theproposed MSM in this context. In particular, we allow for differential availability of treatments acrossstudies, measured confounding within and across studies, and random effects by study. K eywords Conditional average treatment effect, double robustness, individual patient data, marginal structural model,meta-analysis, multidrug-resistant tuberculosis, targeted maximum likelihood estimation
Multidrug-resistant tuberculosis (MDR-TB), a form of tuberculosis (TB) with high mortality, is caused by bacteriaresistant to at least the two most effective anti-TB drugs, isoniazid and rifampicin. According to the World HealthOrganization (WHO) (2020), a global total of 206 030 patients with MDR-TB or rifampicin-resistant(RR) TB were a r X i v : . [ s t a t . M E ] J a n PREPRINT - J
ANUARY
19, 2021detected and notified of their infection in 2019, a increase compared to cases in 2018. However, the latest datareported to WHO show a treatment success rate for MDR/RR-TB of globally.[1] Treating MDR-TB is challengingas a result of the heterogeneity of patients’ characteristics (age, sex, HIV or other comorbidities), disease characteristics(extent and prior treatment), mycobacteria itself (different patterns of additional resistance) and characteristics ofdrugs (more toxicity and less effect of second-line drugs). [2, 3] Patients are typically prescribed a combination offour or more antimicrobial agents depending on the therapeutic phase and drug resistance pattern, if known.[4] Inaddition, the effect of a treatment regimen may vary by an individual’s characteristics and the specific combination ofmedications. When the effect and drug resistance pattern may vary from individual to individual, precision medicinecan be improved by identifying patient covariates to estimate the size and direction of the effect at the individual level.In other words, identifying effect modifiers and assessing effect modification between different patient subgroupsshould be considered. However, this task is statistically challenging and typically requires large amounts of data so thattreatment effects may be well-estimated for different combinations of covariate values. One may impose a workingmodel in order to smooth (or summarize) the covariate-specific effects rather than estimate a separate effect for eachpossible combination of patient covariates.[5, 6, 7] When working with observational data one must also adjust for allpotential confounders of the treatment-outcome relationship, which can be accomplished via propensity scores and/oroutcome regression modeling.[8] One way to model effect modification in a simple binary treatment setting is througha marginal structural model (MSM) for the conditional average treatment effect (CATE).[7, 9] This model may beinterpreted as the relationship between covariates and the expected treatment effect where the treatment effect is definedthrough a contrast of counterfactual outcomes. In non-meta-analytical settings, doubly robust estimators have beenproposed for the estimation of a parametric MSM for the CATE [7] as well as for nonparametric CATE models. [10]Due to the large data requirements for estimating effects across patient subgroups, investigators may be interested inusing individual participant data from multiple studies to fit these treatment effect models. In our study, the data wereextracted from 31 observational studies [11] which contrasted different treatment regimens for patients with MDR-TB,where multiple antimicrobial agents are taken concurrently by a patient over a long period. Our objective is to performan individual participant data (IPD) meta-analysis [12] to investigate the impact of different patient and treatmentcharacteristics on the average treatment effect (ATE) of 15 anti-TB medications.In this project, we propose a targeted maximum likelihood estimator (TMLE) [13] for the estimation of the parametersof the CATE MSM, in the meta-analytical context. For estimation of ATE, TMLE depends on two components: anoutcome regression conditional on treatment and covariates; and, weights comprised of the inverse of the propensityscore where the propensity score is the probability of treatment conditional on covariates.[14] TMLEs are plug inestimators with asymptotic properties that use a targeting step to optimize the bias-variance trade-off for the targetparameter. In our setting, the estimator we propose allows for differential availability of treatments across studies andrandom effects by study due to measured and unmeasured characteristics of the study-specific populations.In Section 2, we describe the MDR-TB data structure and our parameters of interest. Section 3 introduces the TMLEprocedure. We also describe a clustered influence function-based variance estimator. In Section 4, we present the resultsof simulation studies to demonstrate the properties of the estimator under different scenarios. Then in Section 5, weprovide the results based on TMLE analysis of effect modification for 15 anti-TB medications using the combined IPDof the 31 observational studies.
The application data consist of IPD derived from 31 observational studies resulting in a total of 9290 MDR-TB patients.The updated systematic review [Ahuja] extended three previous systematic reviews.[15, 16, 17] We use a one-stageapproach where we pool all IPD across studies to analyze them in a model accounting for clustering and random effects.This IPD were collected from cohorts of adults and the year of studies ranged from 1995 to 2009. The informationcollected for each patient includes demographics (age and sex), past TB history, clinical characteristics (pre-treatmentsputum smear results for acid-fast bacilli (AFB) and culture, chest radiography, HIV infection, drug susceptibility test(DST) results), anti-microbial medications given, and outcomes.
In the pooled dataset, the binary outcome Y represents the treatment success (the treatment was completed and cured thedisease) versus treatment failure (the patient was still culture positive for MDR-TB, experienced a relapse or died).[18]A patient’s outcome realization is defined as lowercase y ij where ( i, j ) refers to patient i ∈ C j in study j and C j be2 PREPRINT - J
ANUARY
19, 2021the set of indices of patients in study j where j ∈ (1 , , · · · , J ) . In the MDR-TB data, there are pooled data from 31studies, i.e. J = 31 . There are 15 antimicrobial agents observed in the data: ethambutol (EMB), ethionamide (ETO), ofloxacin (OFX),pyrazinamide (PZA), kanamycin (KMA), cycloserine (CS), capreomycin (CAP), para-aminosalicylic acid (PAS),prothionamide (PTO), streptomycin (SM), ciprofloxacin (CIP), amikacin (AMK), later-generation fluoroquinolones(LgFQ), rifabutin (RIF) and group five level drugs (Gp5). LgFQ included levofloxacin, moxifloxacin, gatifloxacinand sparfloxacin.[11] Gp5 comprised of amoxicillin-clavulanate, macrolides (azithromycin, roxithromycin, and clar-ithromycin), clofazimine, thiacetazone, imipenem, linezolid, high dose isoniazid, and thioridazine.[11] We use k toindex the k -th antimicrobial agent, where k = 1 , , . . . , in this case. The binary random variable A ( k ) indicates theexposure to medication k with patient realizations a ( k ) ij .Not all treatments are observed in each study; we assume that a given treatment was available for a given study’sparticipants based on whether we have observed that this treatment was taken by any patient in the given study.[19] Thebinary variable D ( k ) is defined as the treatment availability of treatment k in a given study. Specifically, d ( k ) ij = d ( k ) j = 1 means that the treatment k is available to any subject i in study j , and is true if at least one patient in this study wasprescribed the treatment k , and otherwise d ( k ) ij = 0 . The baseline covariates consist of two study level covariates S (the start year of MDR-TB treatment and the incomegroup of the country of the study) and six individual level covariates W (age, sex, AFB results, HIV infection, cavitationstatus on chest radiography and past TB history).Resistance information based on DST is defined as the binary variable R ( k ) . In this dataset, drug resistance informationis available for eight medications. Thus, we denote r ( k ) = 1 if the patient was found to be resistant to the treatment k and otherwise r ( k ) = 0 which includes the situations where the patient’s infection was susceptible to the treatment orwas not known to be resistant to this treatment. The observed data can be written as O = [ S , W , { A ( k ) , D ( k ) , R ( k ) ; k = 1 , , . . . , } , Y ] . We will rewrite R = { R ( k ) ; k = 1 , · · · , } , similarly for D and A . Then the data structure is O = ( S , D , W , R , A , Y ) . In this data, there is differential availability of treatments across studies. In order to define a generalized parameter ofinterest, we define counterfactual notation under the availability of a given treatment.[19] We define counterfactualexposure to treatment A ( k ) { d ( k ) = 1 } as the patient’s counterfactual usage of treatment k had the patient had access.Then we define the counterfactual outcome Y { d ( k ) = 1 , a ( k ) } = Y { a ( k ) } as the patient’s outcome that would haveoccurred had treatment k been available and taken. As in previous work, we aim to estimate the parameters of interest defined with respect to a global population whichrefers to the union of super-populations specific to each study in this dataset.[20, 19] We define a non-parametricstructural equation model (NPSEM), which assumes a time ordered data generating structure, in the SupplementaryMaterials Appendix A. Then the counterfactual likelihood is derived under NPSEM and the parameter of interest isdefined with respect to this likelihood.In this project, for a given medication k , we are interested in estimating the coefficients of potential effect modifiers V ( k ) of the effect of medication k in an MSM for the CATE.[7] To define the corresponding parameters, we recall thatpatients may take multiple medications concurrently. Given a medication k , we conduct the analysis by treating allother medications as confounders. For ease of notation, we define the adjustment set as X ( k ) = [ S , W , R, { D ( k ∗ ) , A ( k ∗ ) ; f or ∀ k ∗ ∈ (1 , · · · , s.t. k ∗ (cid:54) = k } ] PREPRINT - J
ANUARY
19, 2021with realizations x ( k ) ij for individual patients. Here, the symbol asterisk indicates treatments other than the treatment ofinterest. The set of potential effect modifiers investigated is a subset of the adjustment set X ( k ) . [7, 9] The randomcovariate vector is V ( k ) = { , V ( k )1 , ..., V ( k ) p } . Then we can model the CATE of treatment k , denoted ψ { V ( k ) } , as alinear function of potential effect modifiers such that E [ Y { a ( k ) = 1 } − Y { a ( k ) = 0 }| V ( k ) ] = ψ { V ( k ) ; β V ( k ) } = { V ( k ) } (cid:124) β V ( k ) where the symbol (cid:124) indicates a transpose and β V ( k ) = { β ( k )0 , β ( k )1 , ..., β ( k ) p } ∈ R p +1 . The parameter β ( k ) m , m (cid:54) = 0 represents the difference in the expected causal effect related to a single unit change of the effect modifier V ( k ) m and β ( k )0 is the intercept term.Identifying the parameter of interest requires some assumptions that allow us to write the parameter in terms ofdistributions of the observed data.[19](a) Consistency: The counterfactual outcome Y { d ( k ) = 1 , a ( k ) = 1 } , where treatment k was available to and takenby the patient, is the same as the observed outcome for patients who in fact had access to the treatment and took it.Additionally, the counterfactual outcome had the patient not taken this treatment Y { a ( k ) = 0 } is equal to the observedoutcome for those who had not taken the treatment, either due to unavailability or other reasons. In particular, the firstconsistency assumption may fail if this treatment was only taken temporarily by the patient prior to a more effectiveantimicrobial being substituted in.(b) Positivity: For the CATE to be estimable without extrapolation, we also need positivity assumptions. Specifically, theconditional probability of being treated for each patient given the drug’s availability P r [ A ( k ) { d ( k ) = 1 } = 1 | X ( k ) = x ( k ) , D ( k ) = 1] and the conditional probability of not being treated P r { A ( k ) = 0 | X ( k ) = x ( k ) } must both be positive.This may fail if contraindications exist in the covariates X ( k ) , rendering treatment with A ( k ) impossible. Furthermore, P r { D ( k ) = 1 | S ( k ) = s ( k ) } , the probability of availability of treatment k conditional on the study-level covariates mustalso be positive. We also assume that the probability of treatment availability conditional on all measured covariates isonly a function of the study-level covariates. Positivity is violated when, for example, certain studies occurred in a timeor country where some drugs were not on the market, but only if time or country is a study-level confounder.(c) Transportability: The counterfactual outcomes had the patient had access to the treatment and taken it are independentof treatment availability conditional on measured covariates, i.e. Y ( k ) { d ( k ) = 1 , a ( k ) = 1 } ⊥ D ( k ) | X ( k ) . This meansthat we can use the measured covariates to fit models where the treatment is available and use those model fits toestimate overall effects.(d) In addition, we require unconfoundedness: that the counterfactual outcomes had the patient taken this treatmentbe independent of the treatment assignment conditional on the measured covariates and treatment availability. i.e. Y ( k ) { a ( k ) = 1 } ⊥ A ( k ) | D, X ( k ) . Moreover, the counterfactual outcomes had the patient not taken this treatment, Y ( k ) { a ( k ) = 0 } , are independent of the treatment assignment given the measured covariates.Under these assumptions, the CATE can be written as: ψ { V ( k ) } = E [ E { Y | A ( k ) = 1 , X ( k ) }| V ( k ) ] − E [ E { Y | A ( k ) = 0 , X ( k ) }| V ( k ) ] Both right-hand terms can be estimated from the observed data and thus the CATE is nonparametrically identifiable. Aproof is provided in Supplementary Materials Appendix B.
For ease of notation, we will drop the notation k for this section, and consider A = A ( k ) for a given k with X = X ( k ) the adjustment set for k as previously defined. The doubly robust estimators presented in this section require theestimation of two quantities, the conditional outcome expectation and propensity score. We define the former as Q ( A, X ) = P r { Y ( A = a ) = 1 | X } , the probability of counterfactual treatment success conditional on the baselinecovariates X . The propensity score is g ( A | X ) = P r ( A = a | X ) , i.e. the conditional probability of treatment given X .In our setting with treatment availability variable D and the assumptions listed in Section 2.2, we note that Q (1 , X ) = P r ( Y = 1 | A = 1 , X ) = P r ( Y = 1 | D = 1 , A = 1 , X ) . PREPRINT - J
ANUARY
19, 2021Thus, we can estimate this quantity by fitting a regression model using the subgroup of subjects who received thistreatment, which is necessarily a subset of those who had access to the treatment. On the other hand, Q (0 , X ) = P r ( Y = 1 | A = 0 , X ) includes both patients who did not have access to the given treatment and patients who did butdid not receive this treatment.For the propensity score, the probability of being treated for each patient is: g (1 | X ) = P r ( A = 1 | D = 1 , X ) (cid:124) (cid:123)(cid:122) (cid:125) g P r ( D = 1 | S ) (cid:124) (cid:123)(cid:122) (cid:125) g . The first component g may be estimated by fitting a model using patients who had access to the treatment. The secondpart g can be obtained by regressing treatment availability on study level covariates S . This second regression takesthe study as the unit. Then, we can write the probability of not being treated as g (0 | X ) = 1 − g (1 | X ) = 1 − g · g . The efficient influence function (EIF) for a specific parameter is the influence function that achieves the efficiency inthe given space of semi-parametric models.[21] The EIF defines the linear approximation of any efficient and regularasymptotically linear estimator. The EIF for the coefficients β V in a working i.i.d. model space M is D β V = M − D where D = (cid:20)(cid:26) I A =1 g (1 | X ) − I A =0 g (0 | X ) (cid:27) { Y − Q ( A, X ) } + Q (1 , X ) − Q (0 , X ) − V (cid:124) β V (cid:21) V (1)with normalizing matrix M = − E ( ∂ D ∂ β V ) . A proof is provided in Supplementary Materials Appendix C. An equivalentresult was also given by Rosenblum and van der Laan. [22] The general TMLE procedure was proposed by van der Laan and Rubin.[23] Our proposed procedure takes estimates ofthe conditional expected outcome Q ( a, X ) and updates them using information from the propensity score.The first step is to produce initial estimates for Q ( a, X ) for a = 1 and 0, denoted as Q n (1 , X ) and Q n (0 , X ) ,respectively. For each of a = 1 and 0, we run a weighted logistic regression of Y with offset logit { Q n ( a, X ) } and covariates corresponding to the set of potential effect modifiers. The weights are A/g n (1 | X ) for a = 1 and (1 − A ) /g n (0 | X ) for a = 0 , respectively. We then set the updated Q ∗ n ( a, X ) equal to the predicted values from theabove logistic regression. Finally, we fit a linear regression of Q ∗ n (1 , X ) − Q ∗ n (0 , X ) on the potential effect modifiersin order to obtain the TMLE estimates, ˆ β T MLE V , of the parameters of interest.We note that this TMLE solves the equation (cid:80) Jj =1 (cid:80) i ∈ C j D ij,n (ˆ β T MLE V ) = 0 . The estimator thus has the propertiesof double robustness and local efficiency. We also describe a closely related augmented inverse probability of treatmentweighted estimator (A-IPTW) [24] in Supplementary Materials Appendix D. In finite samples, TMLE has been shownto perform better than A-IPTW if the positivity assumption is nearly violated or the true values of the parameters ofinterest are close to the parameter space boundaries.[25] Standard errors can be estimated for the TMLE using a large-sample sandwich estimator of the efficient influencefunction under consistency of Q n ( a, X ) and g n ( a | X ) .[21] Let β V be the true value of β V and ˆ β V be the TMLEestimate.Under regularity conditions, we can write the linear approximation of the estimator as [19] √ n ( ˆ β V − β V ) ≈ √ n J (cid:88) j =1 (cid:88) i ∈ C j D ij ( β V ) where D ij ( β V ) is the influence function at the true values of Q ( a, X ij ) and g ( a | X ij ) for each subject.5 PREPRINT - J
ANUARY
19, 2021In order to estimate the variance of the parameter of interest while taking clustering by studies into account, weonly assume independence between studies, and not individuals within the same study. Within study j , we denotethe ( p + 1) × ( p + 1) dimension variance-covariance matrix of the efficient influence function as σ j . We denotethe ( p + 1) × ( p + 1) dimension variance-covariance matrix of the efficient influence function of any two differentsubjects in study j as ρ j . Both of these quantities can be estimated from the observed data. Then, for large J , thevariance-covariance matrix of ˆ β V can be estimated using: [26] V ar ( ˆ β V ) ≈ n diag (cid:26) J (cid:88) j =1 (cid:18) (cid:88) i,m ∈ C j ,i (cid:54) = m E [ D ij,n (ˆ β V ) {D mj,n (ˆ β V ) } (cid:124) ]+ (cid:88) i ∈ C j E [ D ij,n (ˆ β V ) {D ij,n (ˆ β V ) } (cid:124) ] (cid:19)(cid:27) = 1 n diag (cid:20) J (cid:88) j =1 (cid:26) n j ( n j − ρ j + n j σ j (cid:27)(cid:21) where n j is the size of study C j . This variance estimator is only valid for larger numbers of clusters and consistentestimates of both the outcome and the propensity score. [19] We investigate its finite-sample performance in thesimulation study. In past work, we have also found that the nonparametric clustered bootstrap performs well.[26, 19] We conducted simulation studies in R to demonstrate the double robustness and finite-sample performance of the TMLE.Results for the A-IPTW estimator are provided in the Supplementary Materials Appendix E.
For each dataset, we generated two continuous study-level covariates, S and S , where we treated S as unobserved.We also generated three individual-level covariates: W continuous and W and W binary. In this simulation, weincluded three treatments where each study had access to one, two or all three. We denote the treatment availabilityas D ( k ) , k = 1 , , which was generated conditional on the study level covariate S . Then three treatment indicators A ( k ) were generated based on the values of S , W , W , W and D ( k ) . Subjects could take any combination of thesetreatments. Finally, we generated a continuous outcome Y in a model with and without random effects, conditionalon individual-level covariates, all three treatments, and study-level covariate S . Specifically, this model includedinteractions between the two effect modifiers W and W and A (1) . For random effects by study we added an additionalinteraction term between S and A (1) . Table S1 in Supplementary Material Appendix E displays the full data generatingmechanism. The observed data structure for each subject is O = ( S , W , W , W , D, A, Y ) .In this simulation study, we only aimed to estimate the parameters representing effect modification of treatment k = 1 .For both scenarios with the outcomes involving random effects or not, we drew 1000 simulations with J ∈ { , , } studies where each study contained 300 subjects. This resulted in three total sample sizes, n = 3000 , , . In order to model the effect modification of treatment k = 1 , we included treatment covariates A (2) and A (3) asconfounders, as explained in Section 2.2. We denote the baseline covariates as X (1) = { S , W , W , W , A (2) , A (3) } .In practice we are not aware of the true set of effect modifiers so we define the potential set as { W , W , W , A (2) , A (3) } and set V (1) = { , W , W , W , A (2) , A (3) } .As we discussed in Section 2.2, we model the CATE of treatment A (1) as a linear function of potential effect modifiers: ψ { V (1) ; β V (1) } = { V (1) } (cid:124) β V (1) = β (1)0 + β (1)1 W + ... + β (1)5 A (3) . Our parameters of interest are thus β V (1) . In thescenarios without random effects, the true values of these parameters were derived analytically. For the scenarios withrandom effects, we generated data with large sample sizes ( ) and forced all A (1) equal to 1 and 0, respectively, toobtain both counterfactual outcomes. Then we obtained the true values of the parameters by fitting the linear regressionof the differences between the two counterfactual outcomes on the potential effect modifiers. Given the data generatingmechanism there were two real effect modifiers, W and W . When outcomes were simulated without random effects,the true corresponding coefficient values were . and . , respectively. When outcomes were generated with randomeffects, the true values were . and . (Tables S2, S3 in the Supplementary Materials Appendix E).6 PREPRINT - J
ANUARY
19, 2021We used logistic regressions to estimate each component of g { A | X (1) } , collectively referred to as g , and ¯ Q { a, X (1) } ,referred to as Q . To control the potential sources of sparsity, predicted values for g and g were truncated at ( α, − α ) where α = 0 . . To demonstrate the double robustness property of both methods, we ran four different scenarios:1) correctly specified parametric models for both Q and g ; 2) only the Q model correctly specified and the g modelmisspecified as a null model; 3) only the g model correctly specified and the Q model misspecified as a null model; 4)both Q and g misspecified as null models.In this simulation, we applied the TMLE algorithms presented in Section 3.3. Note that the TMLE implementationwe used requires that we transform the continuous Y to lie in (0 , , and then we reverse-transform at the end ofthe procedure.[27] The standard errors were estimated by the influence function sandwich estimator, first ignoringclustering and then incorporating clustering as described in Section 3.4. We compared both standard error estimates tothe Monte Carlo standard errors. Then based on the clustered standard errors, we constructed Wald-type confidenceintervals and computed the corresponding coverage rates which are the percentage of times that the confidence intervalscontained the true parameter values.
Figure 1 presents the results of TMLE for the estimation of five potential effect modifiers under the four estimationscenarios and three different sample sizes with outcomes generated without random effects. Figure 2 presents the TMLEresults under random effects. The coverage rates are presented in the blue boxes of the two figures. More detailedTMLE results are provided in the Supplementary Materials Appendix E (Tables S2 - S4), along with the figures andtables of the results of the A-IPTW estimator (Tables S5 - S7 and Figures S1, S2).From Figures 1 and 2, we see that under the first two scenarios, where the model for Q was correctly specified, theTMLE estimators had no error on average regardless of the presence of random effects. Unsurprisingly, under scenario4 when all quantities were assigned null models, a small bias was present without random effects and a larger bias waspresent with random effects. In scenario 3, where the model for Q was incorrectly specified but the model for g wascorrect, the average error converged to zero as the number of studies increased. The estimates were also more dispersedthan for the previous scenarios. This occurred because the model for g is estimated using the study as the unit in theanalysis, so with few studies, the sample size to fit this model is extremely small. Compared to TMLE, the A-IPTWestimator had greater average error in scenarios 3 and 4, but otherwise performed similarly.The coverage rates, given in the blue boxes of Figures 1 and 2, typically increased with the number of studies. Forinstance, with random effects with 10 studies, in scenario 1 where all models were correctly specified, the coveragerates for the five potential effect modifiers were between . and . . With 30 studies, the rates increased to . − . . Then for 50 studies, the coverage rates were . − . . Similar patterns were also obtained inthe second and third scenarios. The combined dataset contained the IPD from 31 observational studies with a total of 9290 patients. After removing260 ( . ) patients without reported outcomes, 9030 patients taking 15 different antimicrobial agents remained.Study-specific sample sizes ranged from 25 to 2182 patients. Missing values were present in the covariates; the cavitystatus variable had the most missing values ( . ). Appendix F Table S8 and Figure S3 give summaries of thetreatment-specific sample sizes and the six individual level covariates. Ofloxacin, pyrazinamide and cycloserine werethe three most prescribed medications. In contrast, fewer than 1000 patients were prescribed amikacin, later-generationfluoroquinolones and ciprofloxacin. Ethambutol and pyrazinamide were widely prescribed in 30 out of 31 studies,while only 14 studies had patients who took ciprofloxacin or later-generation fluoroquinolones. The number of malepatients was around twice the number of female patients in all treatment groups. A total of 4892 ( . ) patientswere concentrated in the − year-old age group while only 154 ( . ) were in the − year-old group. Forthe other individual level covariates, about of patients were diagnosed with TB in the past. Moreover, the majorityof patients had cavity ( . ) and positive AFB ( . ), but there were only coinfected with HIV.In order to assess the data support needed to investigate effect modification by concurrent medication, Table 1 displaysthe number of patients who used a combination of any two medications, with the diagonal indicating the total numberof patients taking the corresponding medication. Values ranged between 58 to 4574, indicating at least minimal datasupport for all pairwise combinations. It should be noted that the patients prescribed amikacin all took kanamycinduring the study period. This indicates a positivity violation for the analysis of the effect of kanamycin if we were totreat amikacin as a potential confounder of kanamycin. It also means that we cannot evaluate the effect modification of7 PREPRINT - J
ANUARY
19, 2021 w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o Figure 1: Error of TMLE estimates under four scenarios and three different sample sizes without random effects. The x -axis represents the number of studies. Coverage rates based on the clustered sandwich estimators of the standarderror are presented in blue boxes. The four scenarios are as follows: Scenario 1 - both Q and g models are correct;Scenario 2 - Q model is correct, g is null; Scenario 3 - Q model is null, g model is correct; Scenario 4 - both Q and g models are null. 8 PREPRINT - J
ANUARY
19, 2021 w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o
75% 85.3% 87.6% 85.5% 91.7% 92.8% 89.9% 92.7% 93.8% 88.3% 94.9% 96.2% 83.8% 89.7% 87.2% w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o
90% 89.1% 88.1% 89.2% 85.9% 74.4% 84.9% 82.9% 80.6% 40.1% 30.5% 18.2% 82.2% 90.9% 91.7% w1 w2 w3 a2 a310 30 50 10 30 50 10 30 50 10 30 50 10 30 50−0.50−0.250.000.250.50 S c ena r i o Figure 2: Error of TMLE estimates under four scenarios and three different sample sizes with random effects. The x -axisrepresents the number of studies for three sample sizes. Coverage rates based on the clustered sandwich estimators ofthe standard error are presented in blue boxes. The four scenarios are as follows: Scenario 1 - both Q and g models arecorrect; Scenario 2 - Q model is correct, g is null; Scenario 3 - Q model is null, g model is correct; Scenario 4 - both Q and g models are null. 9 PREPRINT - J
ANUARY
19, 2021amikacin by kanamycin. Therefore, we pragmatically excluded amikacin both as a confounder and as an effect modifierin the analysis of kanamycin. And in the analysis of amikacin, kanamycin was excluded as a potential effect modifierbut was still included as a confounder.Table 1: Summary of number of patients taking any combinations for any two medications during the treatment period.The diagonal values represent the total number of patients taking each medication.
EMB AMK CAP CIP CS ETO OFX PAS PTO RIF SM PZA KMA LgFQ Gp5EMB 4188
266 734 565 1617 2472 2995 989 770 1080 634 3518 2706 272 768
AMK
181 129 416 167 341 271 187 154 58 340 558 128 339
CAP
734 181
482 1719 900 1328 1386 722 325 151 1119 543 205 923
CIP
565 129 482
782 682 236 616 223 350 232 645 481 127 555 CS ETO
OFX
PAS
989 271 1386 616 3573 1206 2750
PTO
770 187 722 223 3004 240 2566 2292
154 749 1564 1532 449 1065
RIF
406 1133 332 87 195 SM
634 58 151 232 981 309 786 732 749 406
870 192 269 339
PZA
KMA
416 1166
LgFQ
272 128 205 127 745 287 192 644 449 87 269 436 416
Gp5
768 339 923 555 1833 670 1262 1463 1065 195 339 930 1166 511
For each medication k , the target parameters in this application are the coefficients of the MSMs, ψ { V ( k ) ; β V ( k ) } = { V ( k ) } (cid:124) β V ( k ) = β ( k )0 + (cid:88) j =1 V ( k ) j β ( k ) j where V ( k ) = { , V ( k )1 , · · · , V ( k )20 } is the set including six individual level covariates and the 14 medications excludingmedication k . We standardized the continuous variable age. β ( k )0 represents the baseline effect for the reference groupof female patients with mean age (39), negative AFB test, no cavitation on chest radiography, no past TB history, noHIV co-infection and not taking any other of the 14 medications. The associated coefficients β ( k ) j , j = 1 , · · · , represent the difference in the CATE when varying the characteristic V ( k ) j by one unit while holding other covariatesfixed.As noted, there are missing values in the individual-level covariates. We used multiple imputation [28] by chained equa-tions with the MICE package [29] in R to produce imputations then used Rubin’s rules to combine the estimates.[30] Ineach imputed dataset we followed the procedure described in Section 3.3 to fit a TMLE. We estimated the Q componentusing SuperLearner (SL) [31] which is a methodology that uses cross validation to find an optimal convex combinationof the predictions of a library of candidate algorithms defined by the user. We included the following algorithms in theSL library: generalized linear models with penalized maximum likelihood ( glmnet function) [32, 33], with forwardstepwise variable selection ( step function), and with a stepwise procedure based on the Akaike Information Criterion( stepAIC function), respectively. [34] Logistic regressions were used for the g models and LASSO penalties wereadded when the logistic regression failed to converge.In each completed dataset, variance estimates for the coefficients were computed using the sample variance of theinfluence function following the expression in Section 3.4. Finally, since many comparisons made in this analysis, weperformed a multiple testing adjustment of the “significance level” of the p-values to control the false discovery rate viathe method of Benjamini and Hochberg [35] with further details given in the Supplementary Materials Appendix G.Figures 3, 4 and 5 show the estimated coefficients, standard errors and confidence intervals corresponding to theintercept and six individual-level effect modifiers and the other 14 medications. The outputs for estimated coefficientsthat are statistically significant after adjustment are marked in red color and asterisk. Corresponding tables of thenumerical results are also provided in the Supplementary Materials Appendix F (Tables S9 and S10). From Figure 3there is no evidence that the six characteristics modify the treatment effect of any drug. In Figure 4 (AMK plot), we seethat patients prescribed rifabutin had lower estimated effects of amikacin. Patients prescribed streptomycin or amikacinmight have benefited less from cycloserine than patients not taking streptomycin or amikacin (CS plot). In addition,taking cycloserine was associated with a greater estimated effect of ethionamide while capreomycin, kanamycin drugs10 PREPRINT - J
ANUARY
19, 2021were associated with lower estimated effects (ETO plot). Finally, in Figure 5 (PTO plot), we interpret that patientsprescribed ethambutol may have benefited more from prothionamide than patients not taking ethambutol.The empirical distributions of the untruncated propensity scores for all drugs are provided in the SupplementaryMaterials Appendix F (Table S11 and Figure S4). For later-generation fluoroquinolones, we noted very large weightswhich likely yielded the large variability observed in Figure 5, LgFQ plot. In addition, since fewer patients wereprescribed later-generation fluoroquinolones among those who were HIV positive (Appendix F Figure S3), the standarderrors were inflated for the coefficient of HIV in the fluoroquinolones MSM (Figure 3, LgFQ plot). Also, among thosewho took rifabutin, only 87/866 subjects were also prescribed later-generation fluoroquinolones, giving rise to the largestandard error (Figure 5, RIF plot).
In this paper, we developed two related one-stage doubly-robust approaches to the analysis of baseline effect modificationin an IPD meta-analysis. The model space that we considered was nonparametric though the parameters of interestwere defined through a working linear MSM for the CATE. These approaches allowed us to analyze pooled IPD frommultiple studies in order to evaluate how estimated treatment effects may vary depending on the values of patientcovariates. Our past work, which proposed related methods and a TMLE for IPD meta-analysis with multiple treatments,instead estimated a treatment importance metric defined as the difference in adjusted probabilities of treatment successbetween the patients who used each medication and the overall population. [19]Vo et al. illustrated that in IPD meta-analysis, heterogeneity across studies can come from two sources: case-mixheterogeneity, due to effect modification, and beyond case-mix heterogeneity, due to differences in study design andmeasurement.[36] Our methods address this by allowing for differential availability of treatments across studies andrandom effects by study due to measured and unmeasured characteristics of the study-specific populations.In clinical and epidemiological research, model misspecification is always a concern when estimating treatment orexposure effects. Doubly robust methods yield consistent estimators even under misspecification of either the treatmentor the outcome model. In the simulation study, we demonstrated the double robustness property of both TMLE andA-IPTW. We observed similar performance of TMLE and A-IPTW but we did not investigate near-positivity violationsor other scenarios that may differentiate them in finite samples as have others. [25, 37] In addition, we demonstratedthe double robustness of both methods when there are study-specific random effects for the outcome. Finally, weshowed that the proposed confidence intervals, estimated using the influence function sandwich estimator that considersclustering by study, performed well when there were greater than 30 studies in the analysis. Indeed, a limitation of ourapproach is that it relies on a sufficient number of studies to estimate a generalizable parameter. In particular, the abilityto adjust for confounding by treatment availability depends on fitting a model for treatment availability conditional onstudy-level covariates, which is limited by the typically small number of studies in a meta-analysis. Indeed, we observein the simulation study that error may persist when the number of studies is small and the outcome regression model isincorrectly specified, even when the propensity score model components are correctly specified with parametric models.Therefore, our approach should only be undertaken when a larger number of studies are available.The treatment of MDR-TB has been challenging because of its prolonged duration, toxicity, costs and unsatisfactoryoutcomes. [1] Second-line TB medicines used for the treatment of drug-resistant TB include amikacin, ciprofloxacin,cycloserine, capreomycin, ethionamide, later-generation fluoroquinolones, kanamycin, ofloxacin, para-aminosalicyclicacid, prothionamide and streptomycin. We used TMLE to investigate effect modification of each of the 15 observedmedications by other individual-level covariates, including the medications taken concurrently for the treatment ofMDR-TB. Our results suggest that cycloserine may enhance the effects of ethionamide but capreomycin and kanamycinwere associated with reduced effects. In addition, amikacin and streptomycin may reduce the effect of cycloserine andethambutol might boost the effect of prothionamide. For the six individual characteristics, age, sex, acid fast bacillistatus, HIV infection, cavitation status on chest radiography, and past TB history, we did not find any evidence ofeffect modification. However, since our MDR-TB data were identified from studies carried out up to 2009, we have noinformation about both new and repurposed effective anti-TB drugs. [38, 39] Therefore, it would be expected to applyour methods to data from more recently treated patients.Another limitation of our approach is that, because we only considered the effect of intervening on one treatment ata time, we cannot directly address how to select combinations of medications that would be expected to optimizethe probability of treatment success. Previous work evaluated the causal contrasts between different regimens ofconcurrent medications in MDR-TB. [40] Future applications should directly address the more challenging question oftreatment-treatment interactions on the outcome which would directly allow for the evaluation of optimal medicationusage. Other ongoing work in our group involves using LASSO, [33] rather than hypothesis testing, to select the effect11
PREPRINT - J
ANUARY
19, 2021 −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Ethambutol (EMB) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Amikacin (AMK) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Capreomycin (CAP) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Ciprofloxacin (CIP) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Cycloserine (CS) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Ethionamide (ETO) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Ofloxacin (OFX) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Para−aminosalicylic acid (PAS) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Prothionamide (PTO) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Rifabutin (RIF) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Streptomycin (SM) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Pyrazinamide (PZA) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Kanamycin (KMA) −2−101 beta_0 Age Sex AFB Cavity HIV Past_TB
Later−generation fluoroquinolones (LgFQ) −0.250.000.250.50 beta_0 Age Sex AFB Cavity HIV Past_TB
Group 5 drugs (Gp5)
Figure 3: Estimated coefficients and the corresponding confidence interval for 15 medications relative to theintercept and six covariates. None of the coefficients reached statistical significance. y -axis of LgFQ plot. 12 PREPRINT - J
ANUARY
19, 2021 −0.50.00.5 A M K C AP C I P C S E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Ethambutol (EMB) ************** −0.50.00.5 C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Amikacin (AMK) −0.50.00.5 A M K C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Capreomycin (CAP) −0.50.00.5 A M K C AP C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Ciprofloxacin (CIP) ************** ************** −0.50.00.5 A M K C AP C I P E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Cycloserine (CS) ************** ************** ************** −0.50.00.5 A M K C AP C I P C S E M B G p5 K M A Lg F Q O F X PAS P T O P Z A R I F S M Ethionamide (ETO)
Figure 4: Estimated coefficients of potential effect modifications and the corresponding confidence intervals forEMB, AMK, CAP, CIP, CS and ETO. Significant results are shown in red and indicated with ∗ .13 PREPRINT - J
ANUARY
19, 2021 −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q PAS P T O P Z A R I F S M Ofloxacin (OFX) −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X P T O P Z A R I F S M Para−aminosalicylic acid (PAS) ************** −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P Z A R I F S M Prothionamide (PTO) −1.0−0.50.00.51.0 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A S M Rifabutin (RIF) −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O P Z A R I F Streptomycin (SM) −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 K M A Lg F Q O F X PAS P T O R I F S M Pyrazinamide (PZA) −0.50.00.5 A M K C AP C I P C S E M B E T O G p5 Lg F Q O F X PAS P T O P Z A R I F S M Kanamycin (KMA) −2.50.02.5 A M K C AP C I P C S E M B E T O G p5 K M A O F X PAS P T O P Z A R I F S M Later−generation fluoroquinolones (LgFQ) −0.50.00.5 A M K C AP C I P C S E M B E T O K M A Lg F Q O F X PAS P T O P Z A R I F S M Group 5 drugs (Gp5)
Figure 5: Estimated coefficients of potential effect modifications and the corresponding confidence intervals forOFX, PAS, PTO, RIF, SM, PZA, KMA, LgFQ and Gp5. Significant results are shown in red and indicated with ∗ . y -axis of RIF plot and LgFQ plot. 14 PREPRINT - J
ANUARY
19, 2021modifiers in the linear MSM for the CATE. This may improve upon the current work by utilizing a superior approach tovariable selection.Identifying effect modifiers is an important step for estimating subpopulation causal effects that can help guide treatmentdecision making for individual patients. Our findings suggest a way to explore the data and generate hypotheses fordrug combinations which can be tested in randomized controlled trial. However, such analyses require larger amountsof data than the estimation of average treatment effects. We have contributed by extending existing doubly robustmethods to incorporate multiple data sources. Advances in IPD meta-analysis enable researchers to incorporate multiplesources of previously collected observational data in their analyses in order to increase their power, which is greatlybeneficial for the identification of effect modifiers.
References [1] World Health Organization.
Global tuberculosis report 2020 . Geneva: World Health Organization, 2020.[2] Müller B, Borrell S, Rose G, et al. The heterogeneous evolution of multidrug-resistant Mycobacterium tuberculosis.
Trends Genet . 2013; :160-169.[3] Isaakidis P, Casas EC, Das M, et al. Treatment outcomes for HIV and MDR-TB co-infected adults and children:systematic review and meta-analysis. Int J Tuberc Lung Dis : 969-978.[4] World Health Organization. WHO consolidated guidelines on drug-resistant tuberculosis treatment . Geneva:World Health Organization, 2019.[5] Petersen ML, Deeks SG, Martin JN, et al. History-adjusted marginal structural models for estimating time-varyingeffect modification.
Am J Epidemiol : 985-993.[6] Powers S, Qian J, Jung K, et al. Some methods for heterogeneous treatment effect estimation in high dimensions.
Stat Med : 1767-1787.[7] Bahamyirou A, Schnitzer ME, Kennedy EH, et al. Doubly robust adaptive LASSO for effect modifier discovery.Submitted to journal.[8] Hernán MA, Robins JM. Estimating causal effects from epidemiological data. J Epidemiol : 578-586.[9] Zhao Q, Small DS, Ertefaie A. Selective inference for effect modification via the lasso. arXiv preprintarXiv:1705.08020. 2017.[10] Kennedy EH. Optimal doubly robust estimation of heterogeneous causal effects. arXiv preprint arXiv:2004.14497.2020; Apr 29.[11] Ahuja SD, Ashkin D, Avendano M, et al. Multidrug resistant pulmonary tuberculosis treatment regimens andpatient outcomes: an individual patient data meta-analysis of 9,153 patients. PLoS Med (8): e1001300.[12] Mills EJ, Ioannidis JP, Thorlund K, et al. How to use an article reporting a multiple treatment comparisonmeta-analysis. JAMA. :1246-1253.[13] van der Laan MJ, Rose S.
Targeted learning: causal inference for observational and experimental data.
SpringerScience & Business Media; 2011.[14] van der Laan MJ and Luedtke AR. Targeted learning of an optimal dynamic treatment, and statistical inference forits mean outcome.
UC Berkeley Division of Biostatistics Working Paper Series . Working Paper 2014; 317.[15] Johnston JC, Shahidi NC, Sadatsafavi M, et al. Treatment outcomes of multidrug-resistant tuberculosis: asystematic review and meta-analysis.
PloS One (9): e6914.[16] Orenstein EW, Basu S, Shah NS, et al. Treatment outcomes among patients with multidrug-resistant tuberculosis:systematic review and meta-analysis. Lancet Infect Dis : 153–161.[17] Akcakır Y. Correlates of treatment outcomes of multidrug-resistant tuberculosis (MDR-TB): a systematic reviewand meta-analysis . MSc Thesis. McGill University, CA. 2010.[18] Falzon D, Jaramillo E, Schünemann HJ, et al. WHO guidelines for the programmatic management of drug-resistanttuberculosis: 2011.[19] Wang G, Schnitzer ME, Menzies D, et al. Estimating treatment importance in multidrug-resistant tuberculosisusing Targeted Learning: An observational individual patient data network meta-analysis.
Biometrics (3):1007-1016.[20] Schnitzer ME, Steele RJ, Bally M, et al. A causal inference approach to network meta-analysis. J Causal Inference .2016; (2). 15 PREPRINT - J
ANUARY
19, 2021[21] Tsiatis, A.
Semiparametric theory and missing data.
Springer Science & Business Media, 2007.[22] Rosenblum M and van der Laan MJ. Targeted maximum likelihood estimation of the parameter of a marginalstructural model.
Int J Biostat (2): Article 19.[23] van der Laan MJ and Rubin D. Targeted maximum likelihood learning. Int J Biostat (1): Article 11.[24] Scharfstein DO, Rotnitzky A, Robins JM. Adjusting for nonignorable drop-out using semiparametric nonresponsemodels. J Am Stat Assoc :1096-1120.[25] Porter KE, Gruber S, van der Laan MJ, et al. The relative performance of targeted maximum likelihood estimators. Int J Biostat : 1-34.[26] Schnitzer ME, van der Laan MJ, Moodie EE, et al. Effect of breastfeeding on gastrointestinal infection in infants:a targeted maximum likelihood approach for clustered longitudinal data. Ann Appl Stat : 703-725.[27] Gruber S and van der Laan MJ. A targeted maximum likelihood estimator of a causal effect on a boundedcontinuous outcome. Int J Biostat (1): Article 26.[28] Rubin, DB. Multiple imputation for nonresponse in surveys.
New York: Wiley, 1987.[29] Buuren SV, Groothuis-Oudshoorn K. mice: Multivariate imputation by chained equations in R.
J. Stat. Softw
Multiple imputation for nonresponse in surveys.
John Wiley & Sons, 2004.[31] van der Laan MJ, Polley EC, Hubbard AE. Super learner.
Stat Appl Genet Mol Biol
J.Stat. Softw : 1.[33] Tibshirani R. Regression shrinkage and selection via the lasso: a retrospective. J R Stat Soc Series B Stat Methodol : 273-282.[34] Venables WN, Ripley BD. Modern applied statistics with S-PLUS.
Springer Science & Business Media; 2013.[35] Benjamini Y and Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multipletesting.
J R Stat Soc Series B Stat Methodol : 289-300.[36] Vo TT, Porcher R, Chaimani A, et al. A novel approach for identifying and addressing case-mix heterogeneity inindividual participant data meta-analysis. Res Synth Methods : 582–596.[37] Schnitzer ME, Moodie EE, Platt RW. Targeted maximum likelihood estimation for marginal time-dependenttreatment effects under density misspecification. Biostatistics : 1-4.[38] Pym AS, Diacon AH, Tang SJ, et al. Bedaquiline in the treatment of multidrug-and extensively drug-resistanttuberculosis. Eur Respir J : 564-574.[39] Sotgiu G, Centis R, D’Ambrosio L, et al. Efficacy, safety and tolerability of linezolid containing regimens intreating MDR-TB and XDR-TB: systematic review and meta-analysis. Eur Respir J : 1430-1442.[40] Siddique AA, Schnitzer ME, Bahamyirou A, et al. Causal inference with multiple concurrent medications: Acomparison of methods and an application in multidrug-resistant tuberculosis. Stat Methods Med Res28