Modeling treatment effect modification in multidrug-resistant tuberculosis in an individual patientdata meta-analysis.

Effect modification occurs while the effect of the treatment is not homogeneous across the different strata 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 direction of the effect at the individual level. However, this task is statistically challenging and typically requires large amounts of data. Investigators may be interested in using the individual patient data from multiple studies to estimate these treatment effect models. Our data arise from a systematic review of observational studies contrasting different treatments for multidrug-resistant tuberculosis, where multiple antimicrobial agents are taken concurrently to cure the infection. We propose a marginal structural model for effect modification by different patient characteristics and co-medications in a meta-analysis of observational individual patient data. We develop, evaluate, and apply a targeted maximum likelihood estimator for the doubly robust estimation of the parameters of the proposed marginal structural model in this context. In particular, we allow for differential availability of treatments across studies, measured confounding within and across studies, and random effects by study.

Introduction

Multidrug-resistant tuberculosis (MDR-TB), a form of tuberculosis (TB) with high mortality, is caused by bacteria resistant to at least the two most effective anti-TB drugs, isoniazid and rifampicin. According to the World Health Organization (WHO) (2020), a global total of 206,030 patients with MDR-TB or rifampicin-resistant (RR)-TB were detected and notified of their infection in 2019, a increase compared to cases in 2018. However, the latest data reported to WHO show a treatment success rate for MDR/RR-TB of globally. Treating MDR-TB is challenging as 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 of drugs (more toxicity and less effect of second-line drugs).[2,3] Patients are typically prescribed a combination of four or more antimicrobial agents depending on the therapeutic phase and drug resistance pattern, if known. In addition, the effect of a treatment regimen may vary by an individual’s characteristics and the specific combination of medications. When the effect and drug resistance pattern may vary from individual to individual, precision medicine can 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 subgroups should be considered. However, this task is statistically challenging and typically requires large amounts of data so that treatment effects may be well-estimated for different combinations of covariate values. One may impose a working model in order to smooth (or summarize) the covariate-specific effects rather than estimate a separate effect for each possible combination of patient covariates.[5-7] When working with observational data one must also adjust for all potential confounders of the treatment-outcome relationship, which can be accomplished via propensity scores and/or outcome regression modeling. One way to model effect modification in a simple binary treatment setting is through a marginal structural model (MSM) for the conditional average treatment effect (CATE).[7,9] This model may be interpreted as the relationship between covariates and the expected treatment effect where the treatment effect is defined through a contrast of counterfactual outcomes. In non-meta-analytical settings, doubly robust estimators have been proposed for the estimation of a parametric MSM for the CATE as well as for nonparametric CATE models.

Due to the large data requirements for estimating effects across patient subgroups, investigators may be interested in using individual participant data (IPD) from multiple studies to fit these treatment effect models. In our study, the data were extracted from 31 observational studies 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 perform an IPD meta-analysis to investigate the impact of different patient and treatment characteristics on the average treatment effect (ATE) of 14 anti-TB medications.

In this project, we propose and evaluate a targeted maximum likelihood estimator (TMLE) to model effect modification for the estimation of the parameters of the CATE MSM, in the IPD meta-analytical context. In Supplemental Materials, we also describe and evaluate a novel augmented inverse probability of treatment weighted estimator (A-IPTW). For estimation of the ATE, TMLE depends on two components: an outcome regression conditional on treatment and covariates; and, weights comprised of the inverse of the propensity score where the propensity score is the probability of treatment conditional on covariates. TMLEs are plug in estimators with asymptotic properties that use a targeting step to optimize the bias-variance tradeoff for the target parameter. In our setting, the estimator we propose allows for differential availability of treatments across studies and random effects by study due to measured and unmeasured characteristics of the study-specific populations.

In the section “Pooled observational studies of MDR-TB,” we describe the MDR-TB data structure and our parameters of interest. The section “Models and algorithms” introduces the TMLE procedure. We also describe a clustered influence function-based variance estimator. In the section “Simulation study,” we present the results of simulation studies to demonstrate the properties of the estimator under different scenarios. Then in the section “MDR-TB data analysis,” we provide the results based on TMLE analysis of effect modification for 14 anti-TB medications using the combined IPD of the 31 observational studies.

Pooled observational studies of MDR-TB

The application data consist of IPD derived from 31 observational studies resulting in a total of 9290 MDR-TB patients. The data are available to any of the data contributors, but not publicly available. The systematic review extended three previous systematic reviews.[15-17] These IPD were collected from cohorts of adults with study years ranging from 1995 to 2009. The information collected for each patient includes demographics (age and sex), past TB history, clinical characteristics (pre-treatment sputum smear results for acid-fast bacilli (AFB) and culture, chest radiography, HIV infection), drug susceptibility test (DST) results, anti-microbial medications given, and outcomes. Our one-stage analytical approach pools all IPD across studies, accounting for clustering and random effects.

Data structure

Outcome

In the pooled dataset, the binary outcome represents the treatment success (the treatment was completed and cured the disease) versus treatment failure (the patient was still culture positive for MDR-TB, experienced a relapse or died). A patient’s outcome realization is defined as lowercase where refers to patient in study and be the set of indices of patients in study where . In the MDR-TB data, there are pooled data from 31 studies, that is, .

Treatments and treatment availabilities

There are 14 antimicrobial agents observed in the data: ethambutol (EMB), ethionamide (ETO), ofloxacin (OFX), pyrazinamide (PZA), kanamycin/amikacin (KM/AM), cycloserine (CS), capreomycin (CAP), para-aminosalicylic acid (PAS), prothionamide (PTO), streptomycin (SM), ciprofloxacin (CIP), later-generation fluoroquinolones (LgFQ), rifabutin (RIF), and group five level drugs (Gp5). LgFQ included levofloxacin, moxifloxacin, gatifloxacin, and sparfloxacin. Gp5 comprised of amoxicillin-clavulanate, macrolides (azithromycin, roxithromycin, and clarithromycin), clofazimine, thiacetazone, imipenem, linezolid, high dose isoniazid, and thioridazine. We use to index the -th antimicrobial agent, where in this case. The binary random variable indicates the exposure to medication with patient realizations .

Not all treatments are observed in each study; we assume that a given treatment was available for a given study’s participants based on whether we have observed that this treatment was taken by any patient in the given study. The binary variable is defined as the treatment availability of treatment in a given study. Specifically, means that the treatment is available to any subject in study , and is true if at least one patient in this study was prescribed the treatment , and otherwise .

Baseline covariates and resistance information

The baseline covariates consist of two study level covariates (the start year of MDR-TB treatment and the income group of the country of the study) and six individual level covariates (age, sex, AFB results, HIV infection, cavitation status on chest radiography, and past TB history).

Resistance information based on DST is defined as the binary variable . In this dataset, drug resistance information is available for eight medications. Thus, we denote if the patient was found to be resistant to the treatment and otherwise which includes the situations where the patient’s infection was susceptible to the treatment or was not known to be resistant to this treatment.

Observed data structure

The observed data can be written as . We will rewrite , similarly for and . Then the data structure is .

Definition of counterfactual notation

In this data, there is differential availability of treatments across studies. In order to define a generalized parameter of interest, we define counterfactual notation under the availability of a given treatment. We define counterfactual exposure to treatment as the patient’s counterfactual usage of treatment had the patient had access. Then we define the counterfactual outcome as the patient’s outcome that would have occurred had treatment been available and taken.

Parameter of interest and assumptions

As in previous work, we aim to estimate the parameters of interest defined with respect to a global population which refers to the union of super-populations specific to each study in this dataset.[20,19] We define a non-parametric structural equation model (NPSEM), which assumes a time ordered data generating structure, in Supplemental Materials Appendix A. Then the counterfactual likelihood is derived under NPSEM and the parameter of interest is defined with respect to this likelihood.

In this project, for a given medication , we are interested in estimating the coefficients of potential effect modifiers of the effect of medication in an MSM for the CATE. To define the corresponding parameters, we recall that patients may take multiple medications concurrently. Given a medication , we conduct the analysis by treating all other medications as confounders. For ease of notation, we define the adjustment set as with realizations for individual patients. Here, the symbol asterisk indicates treatments other than the treatment of interest. The set of potential effect modifiers investigated is a subset of the adjustment set .[7,9] The random covariate vector is . Then we can model the CATE of treatment , denoted , as a linear function of potential effect modifiers such that where the symbol indicates a transpose and . The parameter represents the difference in the expected causal effect related to a single unit change of the effect modifier and is the intercept term.

Identifying the parameter of interest requires some assumptions that allow us to write the parameter in terms of distributions of the observed data.

Consistency: The counterfactual outcome , where treatment was available to and taken by 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 is equal to the observed outcome for those who had not taken the treatment, either due to unavailability or other reasons. In particular, the first consistency assumption may fail if this treatment was only taken temporarily by the patient prior to a more effective antimicrobial being substituted in.

Positivity: For the CATE to be estimable without extrapolation, we also need positivity assumptions. Specifically, the conditional probability of being treated for each patient given the drug’s availability and the conditional probability of not being treated must both be positive. This may fail if contraindications exist in the covariates , rendering treatment with impossible. Furthermore, , the probability of availability of treatment conditional on the study-level covariates must also be positive. We also assume that the probability of treatment availability conditional on all measured covariates is only a function of the study-level covariates. Positivity is violated when, for example, certain studies occurred in a time or country where some drugs were not on the market, but only if time or country is a study-level confounder.

Transportability: The counterfactual outcomes had the patient had access to the treatment and taken it are independent of treatment availability conditional on measured covariates, that is, . This means that we can use the measured covariates to fit models where the treatment is available and use those model fits to estimate overall effects.

In addition, we require unconfoundedness: That the counterfactual outcomes had the patient taken this treatment be independent of the treatment assignment conditional on the measured covariates and treatment availability. that is, . Moreover, the counterfactual outcomes had the patient not taken this treatment, , are independent of the treatment assignment given the measured covariates.

Under these assumptions, the CATE can be written as follows: Both right-hand terms can be estimated from the observed data and thus the CATE is nonparametrically identifiable. A proof is provided in Supplemental Materials Appendix B.

Models and algorithms

Outcome and propensity score models

For ease of notation, we will drop the notation for this section, and consider for a given with the adjustment set for as previously defined. The doubly robust estimators presented in this section require the estimation of two quantities, the conditional outcome expectation and propensity score. We define the former as , the probability of counterfactual treatment success conditional on the baseline covariates . The propensity score is , that is, the conditional probability of treatment given .

In our setting with treatment availability variable and the assumptions listed in the section “Parameter of interest and assumptions,” we note that Thus, we can estimate this quantity by fitting a regression model using the subgroup of subjects who received this treatment, which is necessarily a subset of those who had access to the treatment. On the other hand, includes both patients who did not have access to the given treatment and patients who did but did not receive this treatment.

For the propensity score, the probability of being treated for each patient is The first component may be estimated by fitting a model using patients who had access to the treatment. The second part can be obtained by regressing treatment availability on study level covariates . This second regression takes the study as the unit. Then, we can write the probability of not being treated as .

Efficient influence function

The efficient influence function (EIF) for a specific parameter is the influence function that achieves the efficiency in the given space of semi-parametric models. The EIF defines the linear approximation of any efficient and regular asymptotically linear estimator. The EIF for the coefficients in a working i.i.d. model space is where with normalizing matrix . A proof is provided in Supplemental Materials Appendix C. An equivalent result was also given by Rosenblum and van der Laan.

Targeted maximum likelihood estimator

The general TMLE procedure was proposed by van der Laan and Rubin. Our proposed procedure takes estimates of the conditional expected outcome and updates them using information from the propensity score.

The first step is to produce initial estimates for for and 0, denoted as and , respectively. For each of and 0, we run a weighted logistic regression of Y with offset and covariates corresponding to the set of potential effect modifiers. The weights are for and for , respectively. We then set the updated equal to the predicted values from the above logistic regression. Finally, we fit a linear regression of on the potential effect modifiers in order to obtain the TMLE estimates, , of the parameters of interest.

We note that this TMLE solves the equation . The estimator thus has the properties of double robustness and local efficiency. We also developed a closely related A-IPTW , which is described in Supplemental Materials Appendix D. In finite samples, TMLE has been shown to perform better than a related A-IPTW if the positivity assumption is nearly violated or the true values of the parameters of interest are close to the parameter space boundaries.

Influence function-based variance estimator

Standard errors can be estimated for the TMLE using a large-sample sandwich estimator of the efficient influence function under consistency of and . Let be the true value of and be the TMLE estimate.

Under regularity conditions, we can write the linear approximation of the estimator as where is the influence function at the true values of and for each subject.

In order to estimate the variance of the parameter of interest while taking clustering by studies into account, we only assume independence between studies, and not individuals within the same study. Within study , we denote the dimension variance-covariance matrix of the efficient influence function as . We denote the dimension variance-covariance matrix of the efficient influence function of any two different subjects in study as . Both of these quantities can be estimated from the observed data. Then, for large , the variance-covariance matrix of can be estimated using: where is the size of study . This variance estimator is only valid for a larger number of clusters and consistent estimates of both the outcome and the propensity score. We investigate its finite-sample performance in the simulation study. In past work, we have also found that the nonparametric clustered bootstrap performs well.[26,19]

Simulation study

We conducted simulation studies in R to develop and evaluate TMLE and A-IPTW in the implementation to model effect modification for IPD meta-analysis. We demonstrate the double robustness and finite-sample performance of TMLE and A-IPTW. Results for the A-IPTW estimator are provided in Supplemental Materials Appendix E.

Data generation

For each dataset, we generated two continuous study-level covariates, and , where we treated as unobserved. We also generated three individual-level covariates: continuous and and binary. In this simulation, we included three treatments where each study had access to one, two, or all three. We denote the treatment availability as which was generated conditional on the study-level covariate . Then three treatment indicators were generated based on the values of and . Subjects could take any combination of these treatments. Finally, we generated a continuous outcome in a model with and without random effects, conditional on individual-level covariates, all three treatments, and study-level covariate . Specifically, this model included interactions between the two effect modifiers and and . For random effects by study, we added an additional interaction term between and . Table S1 in Supplemental Material Appendix E displays the full data generating mechanism. The observed data structure for each subject is .

In this simulation study, we only aimed to estimate the parameters representing effect modification of treatment . For both scenarios with the outcomes involving random effects or not, we drew 1000 simulations with studies where each study contained 300 subjects. This resulted in three total sample sizes, .

Analysis

In order to model the effect modification of treatment , we included treatment covariates and as confounders, as explained in section “Parameter of interest and assumptions”. We denote the baseline covariates as . In practice we are not aware of the true set of effect modifiers so we define the potential set as and set .

As we discussed in the section “Parameter of interest and assumptions,” we model the CATE of treatment as a linear function of potential effect modifiers: . Our parameters of interest are thus . In the scenarios without random effects, the true values of these parameters were derived analytically. For the scenarios with random effects, we generated data with large sample sizes ( ) and forced all equal to 1 and 0, respectively, to obtain both counterfactual outcomes. Then we obtained the true values of the parameters by fitting the linear regression of the differences between the two counterfactual outcomes on the potential effect modifiers. Given the data generating mechanism there were two real effect modifiers, and . When outcomes were simulated without random effects, the true corresponding coefficient values were and , respectively. When outcomes were generated with random effects, the true values were and (Tables S2, S3 in Supplemental Materials Appendix E).

We used logistic regressions to estimate each component of , collectively referred to as , and , referred to as . To control the potential sources of sparsity, predicted values for and were truncated at where . To demonstrate the double robustness property of both methods, we ran four different scenarios: (1) correctly specified parametric models for both and ; (2) only the model correctly specified and the model misspecified as a null model; (3) only the model correctly specified and the model misspecified as a null model; (4) both and misspecified as null models.

In this simulation, we applied the TMLE algorithms presented in the section “Targeted maximum likelihood estimator.” Note that the TMLE implementation we used requires that we transform the continuous Y to lie in , and then we reverse-transform at the end of the procedure. The standard errors were estimated by the influence function sandwich estimator, first ignoring clustering and then incorporating clustering as described in the section “Influence function-based variance estimator.” We compared both standard error estimates to the Monte Carlo standard errors. Then based on the clustered standard errors, we constructed Wald-type confidence intervals and computed the corresponding coverage rates which are the percentage of times that the confidence intervals contained the true parameter values.

Results

Figure 1 presents the results of TMLE for the estimation of five potential effect modifiers under the four estimation scenarios and three different sample sizes with outcomes generated without random effects. Figure 2 presents the TMLE results under random effects. The coverage rates are presented in the blue boxes of the two figures. More detailed TMLE results are provided in Supplemental Materials Appendix E (Tables S2 to S4), along with the figures and tables of the results of the A-IPTW estimator (Tables S5 to S7 and Figures S1, S2).

Figure 1.

Error of TMLE estimates under four scenarios and three different sample sizes without random effects. The -axis represents the number of studies. Coverage rates based on the clustered sandwich estimators of the standard error are presented in blue boxes. The four scenarios are as follows: Scenario 1 – both and models are correct; Scenario 2 – model is correct, is null; Scenario 3 – model is null, model is correct; Scenario 4 – both and models are null.

Figure 2.

Error of TMLE estimates under four scenarios and three different sample sizes with random effects. The -axis represents the number of studies for three sample sizes. Coverage rates based on the clustered sandwich estimators of the standard error are presented in blue boxes. The four scenarios are as follows: Scenario 1 – both and models are correct; Scenario 2 – model is correct, is null; Scenario 3 – model is null, model is correct; Scenario 4 – both and models are null.

From Figures 1 and 2, we see that under the first two scenarios, where the model for was correctly specified, the TMLE estimators had no error on average regardless of the presence of random effects. Unsurprisingly, under scenario 4 when all quantities were assigned null models, a small bias was present without random effects and a larger bias was present with random effects. In scenario 3, where the model for was incorrectly specified but the model for was correct, the average error converged to zero as the number of studies increased. The estimates were also more dispersed than for the previous scenarios. This occurred because the model for is estimated using the study as the unit in the analysis, so with few studies, the sample size to fit this model is extremely small. Compared to TMLE, the A-IPTW estimator had greater average error in scenarios 3 and 4, but otherwise performed similarly (Table S8 in Supplemental Materials Appendix E).

The coverage rates, given in the blue boxes of Figures 1 and 2, typically increased with the number of studies. For instance, with random effects with 10 studies, in scenario 1 where all models were correctly specified, the coverage rates 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 in the second and third scenarios.

MDR-TB data analysis

Descriptive statistics of MDR-TB

The combined dataset contained the IPD from 31 observational studies with a total of 9290 patients. After removing 260 ( ) patients without reported outcomes, 9030 patients taking combinations of 14 different antimicrobial agents remained. Study-specific sample sizes ranged from 25 to 2182 patients (Supplemental Materials Appendix F Table S9). Missing values were present in the covariates; the cavity status variable had the most missing values ( ). Appendix F Table S10 and Figure S3 give summaries of the treatment-specific sample sizes and the six individual level covariates. Ofloxacin, pyrazinamide, and cycloserine were the three most prescribed medications. In contrast, fewer than 1000 patients were prescribed later-generation fluoroquinolones 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 male patients was around twice the number of female patients in all treatment groups. A total of 4892 ( ) patients were concentrated in the year-old age group while only 154 ( ) were in the year-old group. For the other individual level covariates, about of patients were diagnosed with TB in the past. Moreover, the majority of 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 displays the number of patients who used a combination of any two medications, with the diagonal indicating the total number of patients taking the corresponding medication. Values ranged between 87 and 4574, indicating at least minimal data support for all pairwise combinations.

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	CAP	CIP	CS	ETO	OFX	PAS	PTO	RIF	SM	PZA	KM/AM	LgFQ	Gp5
EMB	4188	734	565	1617	2472	2995	989	770	1080	634	3518	2706	272	768
CAP	734	1874	482	1719	900	1328	1386	722	325	151	1119	543	205	923
CIP	565	482	968	782	682	236	616	223	350	232	645	481	127	555
CS	1617	1719	782	5629	1941	4195	3573	3004	523	981	3234	2900	745	1833
ETO	2472	900	682	1941	3911	3175	1206	240	397	309	3433	3014	287	670
OFX	2995	1328	236	4195	3175	6464	2750	2566	465	786	4574	4191	192	1262
PAS	989	1386	616	3573	1206	2750	3937	2292	293	732	1962	1816	644	1463
PTO	770	722	223	3004	240	2566	2292	3304	154	749	1564	1532	449	1065
RIF	1080	325	350	523	397	465	293	154	1261	406	1133	332	87	195
SM	634	151	232	981	309	786	732	749	406	1366	870	192	269	339
PZA	3518	1119	645	3234	3433	4574	1962	1564	1133	870	6102	3775	436	930
KM/AM	2706	543	481	2900	3014	4191	1816	1532	332	192	3775	5015	416	1166
LgFQ	272	205	127	745	287	192	644	449	87	269	436	416	866	511
Gp5	768	923	555	1833	670	1262	1463	1065	195	339	930	1166	511	2138

EMB: ethambutol; ETO: ethionamide; OFX: ofloxacin; PZA: pyrazinamide; KM/AM: kanamycin/amikacin; CS: cycloserine; CAP: capreomycin; PAS: para-aminosalicylic acid; PTO: prothionamide; SM: streptomycin; CIP: ciprofloxacin; LgFQ: later-generation fluoroquinolones; RIF: rifabutin; Gp5: group five level drugs.

Analysis and results of MDR-TB

For each medication , the target parameters in this application are the coefficients of the MSMs, where is the set including six individual level covariates and the 13 medications excluding medication . We standardized the continuous variable age. represents the baseline effect for the reference group of female patients with mean age (39 years), negative AFB test, no cavitation on chest radiography, no past TB history, no HIV co-infection and not taking any other of the 13 medications. The associated coefficients represent the difference in the CATE when varying the characteristic by one unit while holding other covariates fixed.

As noted, there are missing values in the individual-level covariates. We used multiple imputation by chained equations with the MICE package in R to produce imputations then used Rubin’s rules to combine the estimates. In each imputed dataset, we followed the procedure described in section “Targeted maximum likelihood estimator” to fit a TMLE. We estimated the component using SuperLearner (SL) which is a methodology that uses cross validation to find an optimal convex combination of the predictions of a library of candidate algorithms defined by the user. We included the following algorithms in the SL library: generalized linear models with penalized maximum likelihood (glmnet function),[32,33] with forward stepwise variable selection (step function), and with a stepwise procedure based on the Akaike Information Criterion (stepAIC function), respectively. Logistic regressions were used for the models and LASSO penalties were added when the logistic regression failed to converge.

In each completed dataset, variance estimates for the coefficients were computed using the sample variance of the influence function following the expression in the section “Influence function-based variance estimator.” Finally, since many comparisons made in this analysis, we performed a multiple testing adjustment of the “significance level” of the p-values to control the false discovery rate via the method of Benjamini and Hochberg with further details given in Supplemental Materials Appendix G.

Figures 3 to 5 show the estimated coefficients, standard errors, and confidence intervals corresponding to the intercept and six individual-level effect modifiers and the other 13 medications. The outputs for estimated coefficients that are statistically significant after adjustment are marked in red color and with an asterisk. Corresponding tables of the numerical results are also provided in Supplemental Materials Appendix F (Tables S11 and S12). From Figure 3, there is no evidence that the six characteristics modify the treatment effect of any drug. In Figure 4 (CS plot), patients prescribed streptomycin might have benefited less from cycloserine than patients not receiving streptomycin. In addition, taking cycloserine was associated with a greater estimated effect of ethionamide while capreomycin, kanamycin/amikacin drugs were associated with lower estimated effects (ETO plot).

Figure 3.

Estimated coefficients and the corresponding confidence interval for 14 medications relative to the intercept and six demographic or clinical covariates. None of the coefficients reached statistical significance. # Larger scale for the -axis of the LgFQ plot.

Figure 5.

Estimated coefficients of potential effect modifiers and the corresponding confidence intervals for PAS, PTO, RIF, SM, PZA, KM/AM, LgFQ and Gp5. None of the coefficients reached statistical significance. # Larger scale for the -axis of the RIF and LgFQ plots.

Figure 4.

Estimated coefficients of potential effect modifiers and the corresponding confidence intervals for EMB, CAP, CIP, CS, ETO and OFX. Significant results are shown in red and indicated with an .

The empirical distributions of the untruncated propensity scores for all drugs are provided in Supplemental Materials Appendix F (Table S13 and Figure S4). For later-generation fluoroquinolones, we noted very large weights which likely yielded the large variability observed in Figure 5, LgFQ plot. In addition, since fewer patients were prescribed later-generation fluoroquinolones among those who were HIV positive (Appendix F Figure S3), the standard errors were inflated for the coefficient of HIV in the fluoroquinolones MSM (Figure 3, LgFQ plot). Also, among those who took rifabutin, only 87/866 subjects were also prescribed later-generation fluoroquinolones, giving rise to the large standard error (Figure 5, RIF plot).

To highlight the differences between meta-analysis and individual study results, we compared the results of Mitnick et al., conducted in 1996–2002[36,37] for the effect modification of ethionamide. Given that all 14 medications were observed in only two studies (Appendix F Table S9), Mitnick et al. (710 subjects) and Tupasi et al. (170 subjects), we chose Mitnick et al. to make the comparison. In Figure 6, we see that a single-study TMLE analysis of the data from the study by Mitnick et al. concluded that prothionamide was associated with greater estimated effects of ethionamide while group 5 drugs and rifabutin were associated with lower estimated effects. Both the meta-analysis and the single study results agreed that cycloserine and capreomycin were effect modifiers (positive and negative, respectively). As expected, the estimation in the meta-analysis had lower standard errors than the individual study (Figure 6). The deviation of the results suggests heterogeneity between different studies in the meta-analysis. It is important to note that analysis of an individual study targets a parameter that is interpreted in the individual study’s super-population while the meta-analysis targets the parameter interpreted in the global population. These parameters may not coincide when there is heterogeneity.

Figure 6.

Estimated coefficients of potential effect modifiers and the corresponding confidence intervals for ethionamide (ETO) in the study by Mitnick et al. (left column) and in the meta-analysis (right column). Significant results are shown in red and indicated with an .

Discussion

In this paper, we developed one-stage doubly-robust methods for the analysis of baseline effect modification in an IPD meta-analysis. The model space that we considered was nonparametric though the parameters of interest were defined through a working linear MSM for the CATE. Our approach allowed us to analyze pooled IPD from multiple studies in order to evaluate how estimated treatment effects may vary depending on the values of patient covariates. 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 success between the patients who used each medication and the overall population. Our current work proposes a novel TMLE and A-IPTW for the estimation of effect modification in the described marginal structural model in an IPD meta-analysis.

Vo et al. illustrated that in IPD meta-analysis, heterogeneity across studies can come from two sources: case-mix heterogeneity, due to effect modification, and beyond case-mix heterogeneity, due to differences in study design and measurement. Our methods address heterogeneity by allowing for differential availability of treatments across studies and random 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 or exposure effects. Doubly robust methods yield consistent estimators even under misspecification of either the treatment or the outcome model. In the simulation study, we demonstrated the double robustness property of both TMLE and A-IPTW. We observed similar performance of TMLE and A-IPTW but we did not investigate near-positivity violations or other scenarios that may differentiate them in finite samples as have others.[25,39] In addition, we demonstrated the double robustness of both methods when there are study-specific random effects for the outcome. Finally, we showed that the proposed confidence intervals, estimated using the influence function sandwich estimator that considers clustering by study, performed well when there were greater than 30 studies in the analysis. Indeed, a limitation of our approach is that it relies on a sufficient number of studies to estimate the generalized parameter interpreted in the global population. In particular, the ability to adjust for confounding by treatment availability depends on fitting a model for treatment availability conditional on study-level covariates, which is limited by the typically small number of studies in a meta-analysis. Indeed, we observe in the simulation study that error may persist when the number of studies is small and the outcome regression model is incorrectly 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 is challenging because of its prolonged duration, toxicity, costs and unsatisfactory outcomes. Second-line TB medicines used for the treatment of drug-resistant TB include injectable drugs (capreomycin, streptomycin, kanamycin, or amikacin), fluoroquinolones (ciprofloxacin, ofloxacin, or later-generation fluoroquinolones), ethionamide, or prothionamide, para-aminosalicyclic acid, and cycloserine. The WHO (2019) revised the guidelines and currently suggest that streptomycin and amikacin are to be considered only if DST results confirm susceptibility and adequate measures to monitor for adverse reactions can be ensured. In addition, capreomycin and kanamycin are not to be included in the treatment of MDR-TB patients. Currently, in general, the selection of antimicrobials for the treatment of an individual patient with TB is based, apart from drug availability, on issues such as the mycobacteria’s drug-resistance pattern and assumptions about the chosen drugs (mechanism of action, potential toxicity, known pharmacological interactions, possible development of drug resistance due to previous use, etc.). However, the actual role played by each drug in the therapeutic outcome is difficult to assess because of the large number of possible combinations of antimicrobials and the potential yet undiscovered interactions affecting their pharmacokinetic and/or pharmacodynamic properties. The statistical approach proposed by the present study takes advantage of the IPD meta-analysis in order to unveil treatment effect modifiers, either regarding medications or other variables at the individual-level, with the aim to better treat and understand not only TB, but also other diseases, including those requiring the concurrent use of several drugs such as arterial hypertension, diabetes, etc.

As discussed in section “Parameter of interest and assumptions”, we require unconfoundedness when estimating the causal effect. However, in practice, this condition is non-testable and we are limited to the covariates collected by the original studies. While most important confounders were collected across all studies, we did not receive DST results from many studies. However, the results were available to clinicians who may have responsively changed or added medications, making DST results an important confounder. Our analyzes adjusted for all available drug sensitivity information.

Our results suggest that cycloserine may enhance the effects of ethionamide but capreomycin and kanamycin or amikacin were associated with reduced effects. In addition, streptomycin may reduce the effect of cycloserine. Our findings support the revised WHO guidelines that injectable drugs are not as beneficial as previously believed. The analyzes informing the 2019 WHO guidelines (Ahmad et al.) showed that kanamycin and capreomycin were associated with worse outcomes. In our results, these drugs were associated with worse outcomes in regimes with ethionamide. For the six individual characteristics, age, sex, acid fast bacilli status, HIV infection, cavitation status on chest radiography, and past TB history, they might be potential for effect modifiers since we did not find any evidence of effect modification. Being that MDR-TB an infectious disease, this finding may simply reflect that antimicrobial drugs play a much larger and definitive role for success than other characteristics at the individual level.

Since our MDR-TB data were identified from studies carried out up to 2009, we have no information about both new and repurposed effective anti-TB drugs.[41,42] Upcoming work will apply our 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 at a time, we cannot directly address how to select combinations of medications that would be expected to optimize the probability of treatment success. Previous work evaluated the causal contrasts between different regimens of concurrent medications in MDR-TB. Future applications should directly address the more challenging question of treatment-treatment interactions on the outcome which would directly allow for the evaluation of optimal medication usage. Other ongoing work in our group involves using LASSO, rather than hypothesis testing, to select the effect modifiers in the linear MSM for the CATE. This may improve upon the current work by utilizing a superior approach to variable selection.

Identifying effect modifiers is an important step for estimating subpopulation causal effects that can help guide treatment decision making for individual patients. However, such analyzes require larger amounts of data than the estimation of average treatment effects. This analytic approach can generate hypotheses for drug combinations which can be tested in randomized controlled trials. We have also contributed by extending existing doubly robust methods to incorporate multiple data sources. Advances in IPD meta-analysis enable researchers to incorporate multiple sources of previously collected observational data in their analyzes in order to increase their power, which is greatly beneficial for the identification of effect modifiers.

Click here for additional data file.

Supplemental material, sj-tex-1-smm-10.1177_09622802211046383 for Modeling treatment effect modification in multidrug-resistant tuberculosis in an individual patientdata meta-analysis by Yan Liu, Mireille E Schnitzer, Guanbo Wang, Edward Kennedy, Piret Viiklepp, Mario H Vargas, Giovanni Sotgiu, Dick Menzies and Andrea Benedetti in Statistical Methods in Medical Research