 Research article
 Open Access
 Published:
Performance of the marginal structural cox model for estimating individual and joined effects of treatments given in combination
BMC Medical Research Methodology volume 17, Article number: 160 (2017)
Abstract
Background
The Marginal Structural Cox Model (CoxMSM), an alternative approach to handle timedependent confounder, was introduced for survival analysis and applied to estimate the joint causal effect of two timedependent nonrandomized treatments on survival among HIVpositive subjects. Nevertheless, CoxMSM performance in the case of multiple treatments has not been fully explored under different degree of timedependent confounding for treatments or in case of interaction between treatments. We aimed to evaluate and compare the performance of the marginal structural Cox model (CoxMSM) to the standard Cox model in estimating the treatment effect in the case of multiple treatments under different scenarios of timedependent confounding and when an interaction between treatment effects is present.
Methods
We specified a CoxMSM with two treatments including an interaction term for situations where an adverse event might be caused by two treatments taken simultaneously but not by each treatment taken alone. We simulated longitudinal data with two treatments and a timedependent confounder affected by one or the two treatments. To fit the CoxMSM, we used the inverse probability weighting method. We illustrated the method to evaluate the specific effect of protease inhibitors combined (or not) to other antiretroviral medications on the anal cancer risk in HIVinfected individuals, with CD4 cell count as timedependent confounder.
Results
Overall, CoxMSM performed better than the standard Cox model. Furthermore, we showed that estimates were unbiased when an interaction term was included in the model.
Conclusion
CoxMSM may be used for accurately estimating causal individual and joined treatment effects from a combination therapy in presence of timedependent confounding provided that an interaction term is estimated.
Background
Combining multiple treatments is a common practice in the therapeutic strategy of chronic or infectious diseases in order to strengthen the effect of treatments or to limit the resistance of pathogens to therapies. When an adverse event occurs in a patient taking multiple treatments, the mainstay of therapy would be to discontinue the suspected inducing drug while maintaining others. In this case, precise identification of the causative treatment is essential.
This topic is particularly relevant when performing a safety analysis of treatments in cohort studies. In such studies, the presence of timedependent confounders (i.e. covariates that predict disease progression and treatment initiation) affected by past treatment might lead to biased estimates if conventional regression methods are used [1, 2]. Furthermore, estimating the individual effect of each treatment becomes methodologically challenging when given simultaneously and when treatment changes overtime.
Marginal structural models (MSMs), a class of causal models, have been proposed as a solution to estimate the causal effect of a timedependent treatment in the presence of timedependent confounders [3, 4]. In this approach, the inverse probability of treatment weighted (IPTW) estimation method is used to consistently estimate MSM parameters [5]. In the context of multiple treatments, a seminal work introduced CoxMSM for survival analysis and applied it to estimate the joint causal effect (efficacy) of two timedependent nonrandomized treatments on survival among HIVpositive subjects [6]. IPTW estimation was used to compute stabilized weights related to multiple medication intakes and to balance the treatment groups at each month. The statistically significant beneficial effects observed were consistent with the results of previous randomized clinical trials. More recently, IPTW estimation was used for joint treatment effects of two treatments or the marginal effect of one treatment in a setting where two concurrent treatments are given at a point in time [7]. MSMs have been used to study the direct effect of several exposures on the outcome of interest by controlling interrelation over time between studied exposures and between exposures and timedependent covariates. In this way, Tager et al. [8] simultaneously studied the effects of physical activity and body composition on functional limitation in the elderly by controlling the confusion induced by interrelation of these two variables over time and their relation to other covariates. To more clearly illustrate the differences in methods and their influence on bias, they carried out simulations comparing weighted and unweighted analysis with respect to the true parameters. However, they did not consider interaction in their simulation study. Howe et al. used joint CoxMSM to estimate the joint effects of multiple timevarying exposures (alcohol consumption and injected drug use) on HIV acquisition [9]. LopezGatell et al. used a joint CoxMSM to estimate the effect of incident tuberculosis disease and Highly Active Antiretroviral therapy (HAART) initiation on AIDSrelated mortality [10]. Cole et al. estimated the joint effects of HAART and PCP prophylaxis on time to AIDS or death using marginal structural models [11]. Bodnar et al. estimated the causal effect of 16 different combinations (regimes) of iron treatment throughout pregnancy on the odds of anemia at delivery [12].
Nevertheless, to date, CoxMSM performance in the case of multiple treatments has not been fully explored under different degree of timedependent confounding for treatments or in case of interaction between treatments. While other studies [8, 9, 11] have included interaction between treatments in their analyses, none has specifically focused on the bias generated when interaction is excluded from the estimated model. This latter issue is critical as numerous adverse events are caused by specific drugdrug interactions and would not occur if each drug was taken separately (e.g., interactions with cytochrome P450 3A4 inhibitors and statins).
The goal of this paper is to evaluate the CoxMSM performance for estimating the individual and joined effects of multiple treatments when they are given in combination through simulation studies. For the sake of simplicity we will limit our study to exploring the use of two treatments with a potential interaction between treatments. We will compare results from CoxMSM with estimates obtained using a classic timedependent Cox regression model and provide an application in the context of HIV infection to evaluate the specific effect of protease inhibitors combined (or not) to other antiretroviral medications on the risk of anal cancer in HIVinfected individuals, using CD4 cell count as timedependent confounder.
The paper is structured as follows: Section 2 describes the method used in this work; Section 3 provides the results of simulation studies that estimate the individual effects of two treatments on an adverse event; Section 4 presents an application of the method. We discuss our results in section 5 and finally, we conclude in section 6.
Methods
Notation for the coxMSM
We considered a longitudinal study in which n subjects (labeled i = 1, …, n) entered a study at baseline, given multiple treatments and were followed at regular time intervals from enrollment into the cohort up to M visits or until the event of interest. Visits (labeled m = 0, 1, 2,…, M) were assumed to take place at the beginning of intervals in the form [m, m + 1]. At each interval, the value of the timedependent confounder, the treatments and the event were observed. We used capital letters to represent random variables and lowercase letters to represent possible realizations (values) of random variables. We explored the model with one disease progression marker and considered the case where a subject might be given two treatments. A_{1i} (m) and A_{2i} (m) denote dichotomous variables indicating whether patient I received treatment A_{1} and/or A_{2} at visit m. Accordingly, there are four possible categories for treatment exposure: exposed to both treatments (A_{1i} (m), A_{2i} (m)) = (1, 1), exposed to only one treatment (A_{1i} (m), A_{2i} (m)) = (1, 0) or (A_{1i} (m), A_{2i} (m)) = (0, 1), not exposed (A_{1i} (m), A_{2i} (m)) = (0, 0). We denoted the baseline fixed covariates V = L (0), the timedependent confounder by L_{i} (m), the event (death or side effect) by Y_{i} (m) and the associated failure time variable that may either be exactly observed or interval censored by T_{i}. We used overbar to represent a covariate history up to a visit, i.e. \( {\overline{A}}_i \) (m) = (A_{i} (0), A_{i} (1), … A_{i} (m)) and \( {\overline{L}}_i \) (m) = (L_{i} (0), L_{i} (1), … L_{i} (m)) to indicate treatment and confounder history up to visit m.
The coxMSM with two treatments
We specified the CoxMSM when two treatments are given to a patient:
Where T is the random variable representing a subject’s survival time given the treatment history, \( {\lambda}_{T_{\overline{a}}}\left(mV\right) \) is the hazard of T at visit m among subjects given pretreatment covariates V, λ _{0}(m) is the unspecified baseline hazard at visit m, exp.(β _{1}), exp.(β _{2}) and exp.(β _{3}) are the causal rate ratios for each treatment and their interaction, and exp.(β _{4}) is the rate ratio associated with the vector of baseline covariates.
Simulation method
We performed simulations using a cohort of HIVpositive individuals receiving multiple antiretroviral treatments. We set CD4 count as the timedependent confounding covariate L_{i} and occurrence of an adverse event as outcome Y_{i}. We simulated a data structure where the outcome at visit m depended on current treatment status only. We assumed that the only baseline covariate was the pretreatment value of the confounder L_{i} (0). This section presents data generation and structure.
We assumed that: A_{1i} (m), A_{2i} (m) and L_{i} (m) remained constant during the subsequent interval between visits m and (m + 1); treatment continued once initiated; and there was no loss to followup  thus we considered the case where censoring occurred only at the end of the followup. Figure 1 shows the causal directed acyclic graph corresponding to the structure of simulated data. We considered three different cases of timedependent confounding. In case 1 and 2, the two treatments were predicted by the timedependent confounder and affected the future value of the timedependent confounder but the coefficient of the timedependent confounder in the function of treatment prediction was set to different values for treatment A_{2} in case 1 (= strong confounding) and case 2 (= weak confounding). In case 3, treatment A_{2} was not predicted by the timedependent confounder, which was not affected by that treatment.
Several studies simulated data from CoxMSM under different conditions [13,14,15,16,17,18,19,20]. In our study, we simulated data for two treatments, adapting the data generation processes for one treatment described in Young et al. and Vourli and Touloumi [15, 21]. For each simulated subject we generated: (a) counterfactual survival times, T_{i} ^{0} from an exponential distribution with parameter λ; (b) covariate values at baseline (time 0) as follows: L_{i} (0) = b + c log T_{i} ^{0} + e_{i,0}, where b ∼ N(μ _{b}, σ_{b} ^{2}); e_{i},_{m} ∼ N(0, σ_{e} ^{2}); c is the coefficient that gives the strength of association between the confounder timedependent covariate and the counterfactual survival time; (c) treatments A_{1i} (m) and A_{2i} (m) from a distribution conditional on a function of past variables. For each treatment this function included L_{i} (m), A_{1i} (m1), A_{2i} (m1) and the product of A_{1i} (m1) and A_{2i} (m1). For each m, we generated subsequent values of the covariate L_{i} (m) as a linear function of past variables; (d) finally, we generated the actual survival time T_{i} of each individual.
We considered five different sets of parameters (A, B, C, D, E) for the marginal true effects of treatments on the outcome resulting in a total of 15 subcases (numbered 1A to 3E). Furthermore, for the cases 1 and 2, coefficient of the timedependent confounder in the function of treatment prediction (see Additional file 1) were set to 0.004 (α_{1} = 0.004, ω_{1} = 0.004) and 0.001 (α_{1} = 0.004, ω_{1} = 0.001), implying a strong and a weak confounding, respectively (Table 1). The coefficient c that gives the strength of association between the timedependent confounder and the counterfactual survival time was set to 6 and μ _{b} and σ_{b} were equal to 600 and 200, respectively. σ_{c} was equal to 3.
The event rate varied between 0.1% to 2%. To avoid separation of data, simulations with no event in at least one category of treatment exposure were discarded.
For each set of parameters, we generated 1000 datasets with 1000 patients each. We assumed that each patient in this cohort was followed for 12 months and had monthly clinical followup visits.
Estimation of the cox MSM with two treatments
To fit the CoxMSM and in order to keep the weight variability as low as possible, the stabilized weights (SW) for estimation via IPTW were used in all analyses [4]. Weights related to each treatment and final censoring were computed and multiplied to obtain a final set of weights as follows:
The numerator of the SW is the probability that a subject received observed treatment A_{i} at visit m conditional only on A_{1i}, A_{2i} history and baseline covariate. The denominator is the probability that a subject received observed treatment A_{i} at visit m given A_{1i}, A_{2i} history and timedependent covariate.
Once the weights were computed, we fitted a weighted Cox proportional hazard model to estimate parameters [9]. We used robust variance estimators to estimate standard errors [22]. We implemented this analysis by using the covs option in the timedependent Phreg procedure in SAS [9]. For the situation where interaction was not set to zero, we also examined the model without including the interaction term. To assess the performance of the model, we computed the absolute bias defined as the difference between average simulated estimates and its corresponding true values and the coverage rate defined as the percentage of confidence intervals that included the true value.
Results
Figures 2, 3 and 4 show the bias and the 95% coverage rate of unweighted and weighted treatment effect estimates as the number of events increases for the different cases. Values of mean bias (MB), standard deviations of estimates, root mean squared error (RMSE) and mean coverage rate (MCR) for all cases are presented in the supplementary material.
As shown in Figs. 2, 3 and 4 (see Additional file 2: Table S1 for mean values), weighted analysis yielded the most accurate estimates of the treatment effects. Indeed, weighted analysis yielded unbiased estimates for the treatment effects A_{1}, A_{2} and interaction between treatments in all cases. In contrast, estimates of unweighted analysis were clearly biased in case 1 (for treatment effect A_{1} and A_{2}), case 2 and case 3 (for interaction between treatments). Estimates of unweighted analysis were less biased for interaction between treatment (in case1) and treatment effects A_{1}, A_{2} (in case 2 and case 3). The values of standard deviations were different from RMSE in all cases of unweighted analysis while the weighted analysis produced values of standard deviation identical to that of the RMSE. Furthermore, for estimates of the unweighted analysis, we observed a slight decrease of bias value as the number of events increased in case1 (for treatment effects A_{1} and A_{2} and interaction between treatments), case 2 and 3 (only for interaction between treatments).
For subcases 1B and 1E, weighted estimates obtained when interaction was not included in the model were biased compared to those obtained when the interaction was included in the model (Fig. 5).
Application to exploring the risk of anal cancer associated with exposure to protease inhibitor in HIV1 infected persons from the FHDHANRS CO4 cohort
Recently, Bruyand et al. [23] found a possible association of cumulative PI (protease inhibitor) exposure with a higher risk of anal cancer in HIV1 infected persons. However, these primary analyses did not adjust for CD4 count at treatment initiation and duration of CD4 count <200 cells/μl, known to be associated with the likelihood of receiving PI and the risk of anal cancer [24]. We applied the CoxMSM framework to evaluate the individual and joined effects of PIs given in combination with other antiretroviral treatments (ARVs), on the risk of anal cancer in HIV1infected persons. Data were obtained from the FHDH cohort (French Hospital Database on HIVANRS CO4), a nationwide hospital cohort initiated in 1989 for individuals infected with HIV [25]. We selected all HIV 1infected treatment naïve individuals at enrollment until 2008. Demographic, clinical, laboratory, ARV information, and cancer events were collected at enrollment and at followup visits as reported elsewhere [26]. For illustration purposes, all ARVs other than PIs were grouped in a single category irrespective of drug class. Baseline covariates were age, gender, transmission group, origin (subSaharan vs other), AIDS diagnosis at baseline, CD4 cell count and HIV RNA. Timedependent covariates were AIDS diagnosis, CD4 cell count, HIV RNA. The timedependent confounder was the CD4 cell count. The followup was split into onemonth periods. Treatment and all timedependent covariates were assumed to remain constant within each period. Time zero was the enrollment date in FHDH. Patients were followed until the occurrence of anal cancer, death or the end of followup, whichever occurred first. A total of 72,355 patients (531,823 personyears) were followed. The median age of the study population was 34 years at enrollment in FHDH. Study subjects were 67% male and 79% from SubSaharan origin. Median CD4 cell count and HIV RNA at baseline were 360 cells/μL and 10,095 copies/mL, respectively. The cohort experienced 9972 personyears (PY) of PIs only, 237,323 PY of other ARV and 130,428 PY of PIs and other ARVs, given simultaneously. During the followup, a total of 130 patients (24/100,000 PY) developed anal cancer. The rate of anal cancer was 90/100,000 PY for patients who received PIs only, 27/100,000 PY for those who received other ARVs, 33/100,000 PY for those who received PIs and other ARVs and 9.6/100,000 PY for untreated patients.
To determine whether current CD4 count predicted treatment with PIs and other ARVs, we fitted pooled logistic models for treatment initiation with PIs and other ARVs that included the baseline covariates and the time dependent covariates. CD4 cell count predicted treatment with PIs (Oddsratios (OR) = 2.77 (p < .0001) for low (<200 cells/μL) versus high CD4 cell count (> 500 cells/μL) and OR = 1.39 (p < .0001) for moderate (200–500 cells/μL) versus high CD4 cell count). CD4 cell count also predicted treatment with other ARVs (OR = 5.63 (p < .0001) for low versus high, and 2.00 (p < .0001) for moderate versus high CD4 cell count, respectively). To determine whether the treatments had an impact on the CD4 count, we fitted a linear model for the mean CD4 count (in cells/μl) in the current month given the baseline covariates, PIs and other ARVs in the previous month, and the remaining timedependent covariates in the previous month [6]. As expected, we found that PIs and other ARVs have an impact on the CD4 count, with coefficient estimates by the linear model of 1.44 (p < .0001) and 0.89 (p < .0001), respectively. This exploratory analysis confirmed that CD4 was a potential timedependent confounder affected by past treatment exposure as described in case 1.
Stabilized weights, related to each treatment class (PIs, other ARVs) and censoring, were then constructed using Eqs. (2), (3) and (4). They were estimated using logistic regression models with baseline covariates for the SW numerator and baseline and timedependent covariates for the SW denominator. To reduce the impact of extremely high weights, we truncated the weights at the 1st and the 99th percentiles of their distribution across all personmonths of followup [27]. The SW had a mean of 1.10 and a standard error of 0.37 after truncation at the 1st and the 99th percentiles (Fig. 6).
For the CoxMSM, in addition to treatment variables including interaction, we adjusted for baseline covariates – this weighted model should be considered as the reference model. We also fitted a weighted model without the interaction term. For the standard time dependent Cox model, we adjusted for all baseline covariates and timedependent covariates and interaction was estimated. The product term between PIs and other ARVs would represent interaction in these models only in the absence of bias due to confounding or selection. Conversely, the product term would represent effect measure modification if bias due to confounding was present for only 1 of the 2 treatments [9, 28].
Table 2 presents estimates of hazard ratios for treatment variables. In the reference model, the risk of anal cancer was significantly increased in patients with isolated PI therapy. Based on the weighted model without interaction, none of PIs, other ARVs nor the combination increased the risk of anal cancer. Conversely, all treatment variables appeared to be associated with the risk of anal cancer in the timedependent Cox model – leading to potentially spurious conclusions. Other variables associated with increased risk of anal cancer in the reference model were longer cumulative duration with CD4 count <200 cells/μl and being a MSM (Men who have sex with Men vs Women) – (results not shown – see Additional file 3: Table S2). Our findings suggest that an increased risk of anal cancer, if any, may exist in the specific category of patients taking PI monotherapy.
Discussion
Through simulation study, we explored the performance of the Cox MSM for estimating the individual effects of two treatments given simultaneously. The simulations showed that using a joint CoxMSM in the presence of a time varying confounder yielded unbiased estimates while standard timedependent Cox model yielded biased estimates. Furthermore, we showed the importance of estimating the interaction term when exploring treatment effects from combination therapy.
The strength of our simulation study is twofold: first, we generated data that is suitable for analysis by a CoxMSM and secondly, we applied a data generation process to simulate data for two treatments, while Vourli and Touloumi [15] and Young et al. [15, 21] performed simulations for only one treatment. Furthermore, we generated a data structure where both combined treatments depend on each other by including an interaction term between both treatments in the treatment predictive model. We also considered a realistic situation when a specific adverse event might be caused by two treatments taken simultaneously but not by one treatment taken alone.
Our simulation study has several limitations. First, we considered that the hazard depends only on the current treatments status. However, treatment effects may cumulate over time and depend on the time since exposure [29]. This requires an assessment as to whether the treatment effects cumulate over time when estimating the individual and joined effects of treatments given in combination [18]. Furthermore, with only one timedependent confounder, our simulated setting could be considered unrealistic and too simplistic. Further studies are needed to consider more complex simulated settings with multiple timedependent confounders and complex hazard functions (cumulative treatment). A number of studies have proposed various algorithms of simulating data suitable for fitting CoxMSMs [14, 17, 30] and could be useful in this context. Second, we explored situations where only two treatments or two classes of treatment were administered; however in real life a patient could receive more comedications. Applying this framework to a real situation with more than two treatments could make calculations of stabilized weights more complex as one has to consider multiple and complex interactions between all treatments. Third, our simulations suggested that our results and conclusions are robust with respect to the number of simulated events, and treatment or confounder effects on the hazards. Future simulations should investigate wider ranges of these parameters as well as the potential impact of the sample size, impact of missing values or unmeasured confounder on the results. Fourth, our result confirmed the superiority of the CoxMSM on the standard timedependent Cox model. Other methods could be explored to estimate the individual effect of treatments when given in combination. For example: doubly robust estimation [31, 32], combines inverse probability weighting with regression modeling of the relationship between covariates and the outcome for each treatment in such a way that, as long as either the propensity score or the regression models are correctly specified, the effect of exposure on the outcome will be correctly estimated. Other methods (e.g., targeted maximum likelihood estimation [33, 34], gcomputation, gestimation of structural nested models etc.. [35]) are also potential alternatives. In addition, other choices for estimating weights in multiple treatment settings could be used, e.g. multinomial logistic regression or machine learning methods [36]. Taken together, future studies would be needed to compare our results with these alternative methods. Fifth, comparing the CoxMSM parameter estimate to any conditional treatment effect estimate is not straightforward when noncollapsible measures, such as hazard ratio or odds ratio, are employed [17, 37,38,39]. We did not perform numerical experiments to explore how the marginal and conditional estimates could differ, which is a limitation of our study. However, the difference between the conditional and marginal parameters is expected to be negligible as the event rate in the time intervals under consideration was small [17, 38].
Exchangeability, positivity and correct model specification are three conditions for unbiased estimation of CoxMSM [6, 9, 27, 28]. For the exchangeability, we assumed that the selected covariates are sufficient to adjust for both confounding and selection bias. The limitation is that this is not testable in an observational study. In our study, we did not observe departures from positivity assumption after truncation of weights as the latter is based on a lack of extreme weights [27]. The lack of extreme weights obtained after truncation provides some evidence against model misspecification.
Using the weighted model with interaction, we found a significant association between use of PIs alone and the risk of anal cancer in HIV infected persons. The HR estimates were markedly different from those obtained with the weighted model without interaction – a finding due to a significant interaction between PI and other ARVs (β3 = −1.43 in the weighted model). Compared with the time dependent Cox model, HR estimate for PIs alone was higher in the weighted model with interaction while HR for other treatment variables were lower leading to different conclusions based on statistical testing – however HR from these two models were in the same range of values. This indicates that timedependent confounding might be weak for all treatment variables and that the time dependent Cox dependent estimates are only slightly biased. In previous studies, Bruyand et al. and Chao et al. found an association between PIs use and anal cancer risk [23, 40]. In both cases, multivariable Poisson models were used. The first study did not adjust for CD4 count at initiation and/or cumulative duration of CD4 count <200 cells/μl and the second one adjusted for CD4 count as timedependent covariate but none dealt with complex timedependent confounding. In our application, we used the model (weighted model with interaction) that performed more accurately in our simulations. Nevertheless, the limitation of our application is that we did not take into account the cumulative duration of ARV exposure. This requires further analysis with cumulative duration of ARV as exposure.
In summary, we evaluated the joint CoxMSM for estimating the individual and joined effects of treatments given in combination in observational studies. The CoxMSM performed accurately in a simulation study under all scenarios. Furthermore, the CoxMSM did not perform accurately when an interaction term was not considered in the model. The application of the framework (weighted model with interaction) on real longitudinal data confirmed the results obtained in the simulation study and has shown the utility of the Joint CoxMSM for estimating the individual and joined causal effects of treatments when they are given in combination in observational studies.
Conclusion
CoxMSM may be used for accurately estimating causal individual and joined treatment effects from a combination therapy in presence of timedependent confounding provided that an interaction term is estimated.
Abbreviations
 AIDS:

Acquired immune deficiency syndrome
 ANRS:

Agence Nationale de Recherche sur le Sida et les hépatites virales
 ARVs:

Antiretroviral treatments
 COXMSM:

Marginal structural Cox model
 FHDH:

French Hospital Database on HIV
 HIV:

Human immunodeficiency virus infection
 IPTW:

Inverse Probability of treatment weighted
 MB:

Mean bias
 MCR:

Mean coverage rate
 MSM:

Men who have sex with men
 MSMs:

Marginal structural models
 PI:

Protease inhibitor
 PY:

Personyears
 RMSE:

Root mean squared error
 RNA:

Ribonucleic acid
 SW:

Stabilized weights
References
 1.
Hernan MA, Brumback B, Robins JM. A structural approach to selection bias. Epidemiology. 2004;15:615–25.
 2.
Robins JM. A new approach to causal inference in mortality studies. Mathematical Modelling. 1986;7:1393–512.
 3.
Robins JM. Association, causation, and marginal structural models. Synthese. 1999a;121(1):151–79.
 4.
Robins JM, Hernan MA. Estimation of the causal effects of timevarying exposures. Longitudinal Data analysis. 2008:553–99.
 5.
Robins JM, Hernan MA, Brumback B. Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;11(5):550–60.
 6.
Hernan MA, Brumback B, Robins JM. Marginal structural models to estimate the joint causal effect of nonrandomized treatments. J Am Stat Asso. 2001;96:440–8.
 7.
Ellis AR, Brookhart MA. Approaches to inverseprobabilityoftreatmentweighted estimation with concurrent treatments. J Clin Epidemiol. 2013;66(8 Suppl):S51–6.
 8.
Tager IB, et al. Effects of physical activity and body composition on functional limitation in the elderly: application of the marginal structural model. Epidemiology. 2004;15(4):479–93.
 9.
Howe CJ, et al. Estimating the effects of multiple timevarying exposures using joint marginal structural models: alcohol consumption, injection drug use, and HIV acquisition. Epidemiology. 2012;23(4):574–82.
 10.
LopezGatell H, et al. Effect of tuberculosis on the survival of women infected with human immunodeficiency virus. Am J Epidemiol. 2007;165(10):1134–42.
 11.
Cole SR, et al. Effect of highly active antiretroviral therapy on time to acquired immunodeficiency syndrome or death using marginal structural models. Am J Epidemiol. 2003;158(7):687–94.
 12.
Bodnar LM, et al. Marginal structural models for analyzing causal effects of timedependent treatments: an application in perinatal epidemiology. Am J Epidemiol. 2004;159(10):926–34.
 13.
Havercroft WG, Didelez V. Simulating from marginal structural models with timedependent confounding. Stat Med. 2012;31(30):4190–206.
 14.
Karim ME, et al. On the application of statistical learning approaches to construct inverse probability weights in marginal structural cox models: hedging against weightmodel misspecification. Commun Stat Simul Comput. 2016:1–30.
 15.
Vourli G, Touloumi G. Performance of the marginal structural models under various scenarios of incomplete marker's values: a simulation study. Biom J. 2015;57(2):254–70.
 16.
Westreich D, et al. A simulation study of finitesample properties of marginal structural cox proportional hazards models. Stat Med. 2012;31(19):2098–109.
 17.
Xiao Y, Abrahamowicz M, Moodie EE. Accuracy of conventional and marginal structural cox model estimators: a simulation study. Int J Biostat. 2010;6(2):Article 13.
 18.
Xiao Y, et al. Flexible marginal structural models for estimating the cumulative effect of a timedependent treatment on the hazard: reassessing the cardiovascular risks of Didanosine treatment in the Swiss HIV cohort study. J Am Stat Asso. 2014;109(506):455–64.
 19.
Young JG, et al. Relation between three classes of structural models for the effect of a timevarying exposure on survival. Lifetime Data Anal. 2010;16(1):71–84.
 20.
Young JG, Tchetgen Tchetgen EJ. Simulation from a known cox MSM using standard parametric models for the gformula. Stat Med. 2014;33(6):1001–14.
 21.
Young, J., S. Picciotto, and J.M. Robins, Simulation from structural survival models under complex timevarying data structures. J Am stat asso, 2008.
 22.
Ali RA, Ali MA, Wei Z. On computing standard errors for marginal structural cox models. Lifetime Data Anal. 2014;20(1):106–31.
 23.
Bruyand M, et al. Cancer risk and use of protease inhibitor or nonnucleoside reverse transcriptase inhibitorbased combination antiretroviral therapy: the D: a: D study. J Acquir Immune Defic Syndr. 2015;68(5):568–77.
 24.
Guiguet M, et al. Effect of immunodeficiency, HIV viral load, and antiretroviral therapy on the risk of individual malignancies (FHDHANRS CO4): a prospective cohort study. Lancet Oncol. 2009;10(12):1152–9.
 25.
MaryKrause M, et al. Cohort profile: French hospital database on HIV (FHDHANRS CO4). Int J Epidemiol. 2014;43(5):1425–36.
 26.
Piketty C, et al. Incidence of HIVrelated anal cancer remains increased despite longterm combined antiretroviral treatment: results from the french hospital database on HIV. J Clin Oncol. 2012;30(35):4360–6.
 27.
Cole SR, Hernan MA. Constructing inverse probability weights for marginal structural models. Am J Epidemiol. 2008;168(6):656–64.
 28.
VanderWeele TJ. On the distinction between interaction and effect modification. Epidemiology. 2009;20(6):863–71.
 29.
Csajka C, Verotta D. Pharmacokineticpharmacodynamic modelling: history and perspectives. J Pharmacokinet Pharmacodyn. 2006;33(3):227–79.
 30.
Karim ME, Platt RW. Estimating inverse probability weights using super learner when weightmodel specification is unknown in a marginal structural cox model context. Stat Med. 2017;36(13):2032–47.
 31.
Bang H, Robins JM. Doubly robust estimation in missing data and causal inference models. Biometrics. 2005;61(4):962–73.
 32.
Robins JM, Rotnitzky A, Zhao LP. Estimation of regression coefficients when some Regressors are not always observed. J Am Stat Assoc. 1994;89(427):846–66.
 33.
van Der Laan M. Targeted maximum likelihood based causal inference: part 1. Int J Biostat. 2010;6(2):2.
 34.
van Der Laan M. Targeted maximum likelihood based causal inference: part 2. Int J Biostat. 2010;6(2):3.
 35.
Daniel RM, et al. Methods for dealing with timedependent confounding. Stat Med. 2013;32(9):1584–618.
 36.
McCaffrey DF, et al. A tutorial on propensity score estimation for multiple treatments using generalized boosted models. Stat Med. 2013;32(19):3388–414.
 37.
Austin PC. The performance of different propensity score methods for estimating marginal hazard ratios. Stat Med. 2013;32(16):2837–49.
 38.
Karim ME, et al. Comparison of statistical approaches dealing with timedependent confounding in drug effectiveness studies. Stat Methods Med Res. 2016;
 39.
Pang M, Kaufman JS, Platt RW. Studying noncollapsibility of the odds ratio with marginal structural and logistic regression models. Stat Methods Med Res. 2016;25(5):1925–37.
 40.
Chao C, et al. Exposure to antiretroviral therapy and risk of cancer in HIVinfected persons. AIDS. 2012;26(17):2223–31.
Acknowledgements
Not applicable.
Funding
This work was supported by Agence Nationale de Recherche sur le Sida et les hépatites virales (ANRS). The funding body did not play any role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The data used and analyzed during this study are not publicly available due to confidentiality reasons.
Author information
Affiliations
Contributions
CLN conceived the idea, performed simulations, made analysis and interpretation of results and wrote the draft of the manuscript. HSL made substantial contribution to acquisition of data and was involved in revising critically the manuscript. SG and DC made substantial contribution to analysis and interpretation of results and were involved in revising the manuscript critically. FC conceived the idea and made substantial contribution to analysis and interpretation of results and was involved in writing and revising the manuscript critically. All authors read and gave final approval of the version to be published and agreed to be accountable for all aspects of the work.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1:
Complete data generation. (DOCX 28 kb)
Additional file 2: Table S1.
Mean bias, standard deviation, mean squared error and mean coverage rate of estimates. (DOCX 75 kb)
Additional file 3: Tables S2.
Multivariate parameter estimates for covariate association with the risk of anal cancer in HIVinfected persons: comparison of weighted Cox MSM and standard time dependent Cox models. (DOCX 22 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
LusivikaNzinga, C., SelingerLeneman, H., Grabar, S. et al. Performance of the marginal structural cox model for estimating individual and joined effects of treatments given in combination. BMC Med Res Methodol 17, 160 (2017). https://0doiorg.brum.beds.ac.uk/10.1186/s1287401704341
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s1287401704341
Keywords
 Causal inference
 Timedependent confounding
 Longitudinal data
 Marginal structural models
 Multitherapy