Performance of the marginal structural cox model for estimating individual and joined effects of treatments given in combination

Background The Marginal Structural Cox Model (Cox-MSM), an alternative approach to handle time-dependent confounder, was introduced for survival analysis and applied to estimate the joint causal effect of two time-dependent nonrandomized treatments on survival among HIV-positive subjects. Nevertheless, Cox-MSM performance in the case of multiple treatments has not been fully explored under different degree of time-dependent confounding for treatments or in case of interaction between treatments. We aimed to evaluate and compare the performance of the marginal structural Cox model (Cox-MSM) to the standard Cox model in estimating the treatment effect in the case of multiple treatments under different scenarios of time-dependent confounding and when an interaction between treatment effects is present. Methods We specified a Cox-MSM 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 time-dependent confounder affected by one or the two treatments. To fit the Cox-MSM, 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 HIV-infected individuals, with CD4 cell count as time-dependent confounder. Results Overall, Cox-MSM performed better than the standard Cox model. Furthermore, we showed that estimates were unbiased when an interaction term was included in the model. Conclusion Cox-MSM may be used for accurately estimating causal individual and joined treatment effects from a combination therapy in presence of time-dependent confounding provided that an interaction term is estimated. Electronic supplementary material The online version of this article (10.1186/s12874-017-0434-1) contains supplementary material, which is available to authorized users.


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 time-dependent 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 time-dependent treatment in the presence of time-dependent 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 Cox-MSM for survival analysis and applied it to estimate the joint causal effect (efficacy) of two timedependent nonrandomized treatments on survival among HIV-positive 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 time-dependent 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 Cox-MSM to estimate the joint effects of multiple time-varying exposures (alcohol consumption and injected drug use) on HIV acquisition [9]. Lopez-Gatell et al. used a joint Cox-MSM to estimate the effect of incident tuberculosis disease and Highly Active Antiretroviral therapy (HAART) initiation on AIDS-related 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, Cox-MSM performance in the case of multiple treatments has not been fully explored under different degree of time-dependent 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 drug-drug 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 Cox-MSM 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 Cox-MSM with estimates obtained using a classic time-dependent 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 HIV-infected individuals, using CD4 cell count as time-dependent 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.

Notation for the cox-MSM
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 time-dependent confounder, the treatments and the event were observed. We used capital letters to represent random variables and lower-case 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 time-dependent 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. A i (m) = (A i (0), A i (1), … A i (m)) and L i (m) = (L i (0), L i (1), … L i (m)) to indicate treatment and confounder history up to visit m.

The cox-MSM with two treatments
We specified the Cox-MSM when two treatments are given to a patient: Where T is the random variable representing a subject's survival time given the treatment history, λ T a mjV ð Þ 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 time-dependent 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 pre-treatment 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 follow-up -thus we considered the case where censoring occurred only at the end of the follow-up. Figure 1 shows the causal directed acyclic graph corresponding to the structure of simulated data. We considered three different cases of time-dependent confounding. In case 1 and 2, the two treatments were predicted by the time-dependent 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 time-dependent confounder, which was not affected by that treatment.
Several studies simulated data from Cox-MSM 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: ; c is the coefficient that gives the strength of association between the confounder time-dependent 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 (m-1), A 2i (m-1) and the product of A 1i (m-1) and A 2i (m-1). 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 sub-cases (numbered 1A to 3E). Furthermore, for the cases 1 and 2, coefficient of the time-dependent 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 time-dependent 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. Causal directed acyclic graphs corresponding to the structure of simulated data. A 1 and A 2 are the treatments, L is the time-dependent confounder and Y is the outcome. Case 1 and 2 considered all relationship between A 1 , A 2 and L. The time-dependent Confounder was strongly associated to the treatments A 1 and A 2 in the case 1 whereas it was weakly associated to the treatments A 2 in the case 2. Coefficients of the time-dependent confounder in the functions of treatment prediction were set to 0.004 and 0.001, respectively. Case 3: relationship between A 2 and L were not considered. Data were simulated from a marginal structural model as the confounding in the exposures-outcome relationship arises via T 0 as follows: 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 follow-up visits.

Estimation of the cox MSM with two treatments
To fit the Cox-MSM 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 time-dependent 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 time-dependent 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.  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).

Results
For sub-cases 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 HIV-1 infected persons from the FHDH-ANRS 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 HIV-1 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 Cox-MSM 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 HIV1-infected persons. Data were obtained from the FHDH cohort (French Hospital Database on HIV-ANRS CO4), a nationwide hospital cohort initiated in 1989 for individuals infected with HIV [25]. We selected all HIV 1-infected treatment naïve individuals at enrollment until 2008. Demographic, clinical, laboratory, ARV information, and cancer events were collected at enrollment and at follow-up visits as reported elsewhere [26].  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 time-dependent 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 time-dependent 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 time-dependent 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 person-months of follow-up [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 Cox-MSM, 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 time-dependent 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 time-dependent Cox modelleading 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 shownsee 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 Cox-MSM in the presence of a time varying confounder yielded unbiased estimates while standard time-dependent 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 Cox-MSM 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 time-dependent confounders and complex hazard functions (cumulative treatment). A number of studies have proposed various algorithms of simulating data suitable for fitting Cox-MSMs [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 co-medications. 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 Cox-MSM on the standard time-dependent 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], g-computation, g-estimation 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 Cox-MSM parameter estimate to any conditional treatment effect estimate is not straightforward when non-collapsible 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 Cox-MSM [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 interactiona 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 testinghowever 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 time-dependent covariate but none dealt with complex time-dependent 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 Cox-MSM for estimating the individual and joined effects of treatments given in combination in observational studies. The Cox-MSM performed accurately in a simulation study under all scenarios. Furthermore, the Cox-MSM 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 Cox-MSM for estimating the individual and joined causal effects of treatments when they are given in combination in observational studies.