# Estimating the Effects of Multiple Time-varying Exposures Using Joint Marginal Structural Models: Alcohol Consumption, Injection Drug Use, and HIV Acquisition

The joint effects of multiple exposures on an outcome are frequently of interest in epidemiologic research. In 2001, Hernán et al (*J Am Stat Assoc.* 2001;96:440–448) presented methods for estimating the joint effects of multiple time-varying exposures subject to time-varying confounding affected by prior exposure using joint marginal structural models. Nonetheless, the use of these joint models is rare in the applied literature. Minimal uptake of these joint models, in contrast to the now widely used standard marginal structural model, is due in part to a lack of examples demonstrating the method. In this paper, we review the assumptions necessary for unbiased estimation of joint effects as well as the distinction between interaction and effect measure modification. We demonstrate the use of marginal structural models for estimating the joint effects of alcohol consumption and injection drug use on HIV acquisition, using data from 1525 injection drug users in the AIDS Link to Intravenous Experience cohort study. In the joint model, the hazard ratio (HR) for heavy drinking in the absence of any drug injections was 1.58 (95% confidence interval = 0.67–3.73). The HR for any drug injections in the absence of heavy drinking was 1.78 (1.10–2.89). The HR for heavy drinking and any drug injections was 2.45 (1.45–4.12). The *P* values for multiplicative and additive interaction were 0.7620 and 0.9200, respectively, indicating a lack of departure from effects that multiply or add. We could not rule out interaction on either scale due to imprecision.

Supplemental Digital Content is available in the text.

From the ^{a} Department of Epidemiology, Center for Population Health and Clinical Epidemiology, Brown University Program in Public Health, Providence, RI; ^{b} Department of Epidemiology, University of North Carolina Gillings School of Global Public Health, Chapel Hill, NC; ^{c} Department of Epidemiology, Johns Hopkins Bloomberg School of Public Health, Baltimore, MD; and ^{d} Division of Infectious Diseases, Department of Medicine, Johns Hopkins School of Medicine, Baltimore, MD.

Submitted 19 June 2011; accepted 13 December 2011; posted 10 April 2012.

Supported by National Institutes of Health grants R01-AA-01759, R01-DA-04334, R01-DA-12568.

The authors reported no other financial interests related to this research.

Supplemental digital content is available through direct URL citations in the HTML and PDF versions of this article (www.epidem.com). This content is not peer-reviewed or copy-edited; it is the sole responsibility of the author.

Correspondence: Chanelle J. Howe, Department of Epidemiology, Center for Population Health and Clinical Epidemiology, Brown University Program in Public Health, 121 South Main St, Providence, RI 02912. E-mail: chanelle_howe@brown.edu.

The joint effects of multiple exposures on an outcome are frequently of interest in epidemiologic research.^{1} Joint effects concern the combined effect of 2 or more exposures on an outcome and whether this combined effect is more or less than expected on a given scale (ie, additive or multiplicative). Departures from expected combined effects are referred to as interaction. In contrast, effect measure modification focuses on whether the effect of a single exposure on an outcome varies with the level of second variable.^{2} For example, in the area of cardiovascular disease epidemiology, if the study objective was to estimate the combined effect of oral contraceptive use and smoking on the risk of myocardial infarction, then the research would attempt to estimate the joint effects of (and interaction between) oral contraceptive use and smoking on the risk of myocardial infarction. However, if the objective was whether the effect of oral contraceptive use on myocardial infarction varied by smoking status, then the aim would be to estimate effect measure modification.

In 2001, Hernán et al^{3} presented methods for estimating the joint effects of multiple time-varying exposures subject to time-varying confounding using joint marginal structural models. We refer to these marginal structural models as “joint” models, given that they examine the potential outcome across the joint distribution of the 2 or more exposures (not to be confused with models for the joint distribution of survival and a repeated measure).^{4} Standard marginal structural models^{5},^{6} avoid removal of indirect effects, and they also avoid the selection bias that can occur with standard adjustment techniques when estimating the effect of a time-varying exposure subject to time-varying confounding affected by prior exposure. Joint marginal structural models have the same methodological advantage when estimating the effect of multiple time-varying exposures, each affected by time-varying confounders affected by prior exposure. Both standard and joint marginal structural models can also be used to account for selection bias due to censoring.

Although estimating joint effects of time-varying exposures is of interest in both chronic and infectious disease epidemiology,^{7} ^{–} ^{10} a decade after the paper by Hernán et al^{3} was published, we found only 3 publications^{11} ^{–} ^{13} that used a joint marginal structural model in the applied biomedical and epidemiologic literature among articles citing that paper. This low uptake of joint marginal structural models, in contrast to the widely used standard marginal structural model,^{6},^{14} ^{–} ^{19} may be due to a lack of examples showing the implementation of the model.

The objective of this paper is to describe and demonstrate joint marginal structural Cox proportional hazards models based on an example from the HIV literature. (Although joint marginal structural models can be fit for other types of regression models,^{2},^{3} these other models are not the focus of this paper.) Previous work^{7},^{8},^{19} has indicated alcohol consumption as well as injection drug use as risk factors for HIV acquisition. To identify groups at particularly high risk for HIV infection, the interaction between risk factors for HIV acquisition has been of interest.^{7},^{8} As depicted in Figure 1, both time-varying alcohol consumption and injection drug use are likely subject to time-varying confounding affected by prior exposure. In particular, sexual risk behaviors and other patterns of drug use are likely common causes of alcohol consumption, injection drug use, and HIV acquisition and may be affected by prior alcohol consumption or injection drug use. Therefore, unbiased estimation of the joint effects of these 2 risk factors on HIV acquisition likely requires use of joint marginal structural models.

We used data from the AIDS Link to Intravenous Experience (ALIVE) cohort study^{20} to demonstrate estimation of the joint effects of alcohol consumption and injection drug use on HIV acquisition. The SAS (SAS Institute, Inc., Cary, NC) code used to estimate these joint effects is provided in the eAppendix (http://links.lww.com/EDE/A567). Assumptions necessary for unbiased estimation of joint effects will be described generally, as well in the context of the ALIVE example. Finally, we will make additional clarifications regarding the distinction between interaction and effect measure modification.

## METHODS

### Study Population

Details regarding recruitment, inclusion criteria, and data collection for this example are provided elsewhere.^{19},^{20} Briefly, the ALIVE cohort enrolled 3627 HIV seronegative adults with a history of injection drug use in Baltimore, MD, between 1988 and 2008. Through extensive community outreach, participants were recruited from venues where injection drug users in Baltimore were known to be frequent, as well as from other local settings, such as sexually transmitted disease clinics, drug abuse treatment programs, homeless shelters, and emergency rooms.

Of 3627 enrollees, 1525 were African American, with complete demographic information who were seronegative at entry and who had at least one seronegative follow-up visit subsequent to baseline, thereby making them eligible for the analysis. Blood samples, as well as information on demographic characteristics, drug and alcohol use history and practices, and sexual behaviors, were collected at study entry and every 6 months thereafter. Blood specimens were tested for HIV antibodies by enzyme-linked immunoassay, and reactive specimens were confirmed by Western blot. The 1525 participants were followed for up to 20 years, yielding 8181 person-years of follow-up, during which time 155 acquired HIV.

### Notation

Let uppercase letters denote random variables, and lowercase letters represent possible realizations of random variables. In the ALIVE cohort of *i* = 1, 2, …, *N* = 1525 participants with *k* _{i} _{(1)} < *k* _{i} _{(2)} < …*k* _{i} _{(} _{M} _{)} ordered measurement times since study entry, *T* _{i} is the time in years from study entry to HIV infection and *C* _{i} is the time in years from study entry to censoring due to dropout, death, or reaching the end of the study. Therefore, *Y* _{i} = min(*T* _{i}, *C* _{i}) and δ_{i} is an indicator of whether HIV infection was observed, where δ_{i} = 1 if *y* _{i} = *t* _{i} and δ_{i} = 0 if *y* _{i} = *c* _{i}. Below, the subscript *i* will be suppressed where possible. Let *X* _{1} _{i}(*k* _{(} _{M} _{)}) and *X* _{2} _{i}(*k* _{(} _{M} _{)}) be the level of alcohol consumption and injection drug use status at time *k* _{(} _{M} _{)}, respectively.

Let *L* _{1} _{i}(*k* _{(} _{m} _{)}) be a vector denoting the levels of potential confounders of the relationship between *X* _{1} and *T* at time *k* _{(} _{m} _{)}, while *Z* _{1} _{i} is a vector subset of *L* _{1} _{i}(*k* _{(} _{m} _{)}) corresponding to time-fixed confounders measured at study entry. Let *L* _{2} _{i}(*k* _{(} _{m} _{)}) be a vector denoting the levels of potential confounders of the relationship between *X* _{2} and *T* at time *k* _{(} _{M} _{)}, while *Z* _{2} _{i} is a vector subset of *L* _{2} _{i}(*k* _{(} _{m} _{)}) corresponding to time-fixed confounders measured at study entry. Let *L* _{3} _{i}(*k* _{(} _{m} _{)}) be a vector denoting the levels of common predictors of the relationship between *T* and *C* at time *k* _{(} _{M} _{)}, whereas *Z* _{3} _{i} is a vector subset of *L* _{3} _{i}(*k* _{(} _{m} _{)}) corresponding to time-fixed common predictors measured at study entry. Covariates *L* _{1}, *L* _{2}, and *L* _{3} were assumed to be the same and included those that were time-fixed (ie, age at entry, sex, and years of education) as well as time-updated pertaining to the prior 6 months (ie, cocaine use; shooting gallery attendance; number of sexual partners; male sex with men; sexually transmitted infections, including genital herpes simplex virus, genital warts, gonorrhea, and syphilis). Selected covariates were variables shown to independently increase the risk of HIV infection among study participants.^{7},^{8}

Allow overbars to represent the history of a given random variable from study entry to *k* _{(} _{M} _{)} (eg, *X[Combining Macron]* _{1} _{i}(*k* _{(} _{M} _{)}) = {*X* _{1} _{i}(*k* _{(1)}), *X* _{1} _{i}(*k* _{(2)}), …, *X* _{1} _{i}(*k* _{(} _{M} _{)})}). Therefore, let

represent the observed time to HIV infection for participant *i* had he or she been assigned alcohol consumption and injection drug use history *x[Combining Macron]* _{1} _{i}(*k* _{(} _{M} _{)}) and *x[Combining Macron]* _{2} _{i}(*k* _{(} _{M} _{)}) rather than the observed alcohol consumption and injection drug use history of

and

. Similarly,

represents the observed time to HIV infection for participant i for observed alcohol consumption and injection drug use history X¯_{1} _{i}(k_{(} _{M} _{)}) and X¯_{2} _{i}(k_{(} _{M} _{)}). Finally, let D_{1} _{i}(k_{(} _{M} _{)}) be an indicator that participant i is censored due to dropout at time k_{(} _{m} _{)},D_{2} _{i}(k_{(} _{m} _{)}) be an indicator that participant i is censored due to death at time k_{(} _{m} _{)},X_{i} = [X_{1} _{i},X_{2} _{i}],L_{i} = [L_{1} _{i},L_{2} _{i},L_{3} _{i}],Z_{i} = [Z_{1} _{i},Z_{2} _{i},Z_{3} _{i}],and D_{i} = [D_{1} _{i},D_{2} _{i}].

### Joint Marginal Structural Cox Proportional Hazards Model

Equation (1) shows the marginal structural Cox proportional hazards model for the joint effects of *X* _{1}(*k* _{(} _{M} _{)}) and *X* _{2}(*k* _{(} _{M} _{)}) on *T* ^{x[Combining Macron]},

where *T* ^{x[Combining Macron]} = *T* ^{x[Combining Macron]1} ^{,x[Combining Macron]2}. Also, λ_{0}(*k* _{(} _{M} _{)}) is the baseline hazard function at time *k* _{(} _{M} _{)} and λ_{Tx[Combining Macron]}(*k* _{(} _{M} _{)}|Z) is the hazard function for *T* ^{x[Combining Macron]} at time *k* _{(} _{M} _{)} conditional on *Z*. The parameter α_{1} is the log hazard ratio for the effect of *X* _{1}(*k* _{(} _{M} _{)}) on *Tx[Combining Macron]* when *X* _{2}(*k* _{(} _{M} _{)}) = 0, α_{2} is the log hazard ratio for the effect of *X* _{2}(*k* _{(} _{M} _{)}) on *T* ^{x[Combining Macron]} when *X* _{1}(*k* _{(} _{M} _{)}) = 0, and α_{3} is the log of the ratio of hazard ratios for the effect of *X* _{1}(*k* _{(} _{M} _{)}) and *X* _{2}(*k* _{(} _{M} _{)}) on *T* ^{x[Combining Macron]}. In addition, α_{4}, α_{5}, and α_{6} are the vectors of log hazard ratios for the time-fixed covariates comprising *Z* _{1}, *Z* _{2}, and *Z* _{3}, respectively while α = [α_{1}, α_{2}, α_{3}, α_{4}, α_{5}, α_{6}].

The joint marginal structural model in equation (1) is fit in a two-part process using inverse probability weights (ie, *W* _{i}) described later. Conditional on the estimated weights, the weighted maximum partial likelihood estimates for [Combining Circumflex Accent]α in equation (1) are obtained by maximizing,

where *I*[*y* _{r} ≥ *y* _{i}] is an indicator of participant *r* being in the risk set at follow-up time *y* _{i}. Given necessary assumptions outlined below, α_{1}, α_{2}, and α_{3} are consistent asymptotically normal estimates of the joint effects of *X* _{1} and *X* _{2} on *T* ^{x1,x2} while correcting for potential time-varying confounding by *L* _{1} and *L* _{2} as well as selection bias due to censoring informed by *L* _{3}. Therefore, α_{3} can be used to consistently assess the presence or absence of multiplicative interaction.

### Inverse Probability Weights

The following stabilized weights can be used to fit the joint marginal structural model in equation (1):

where

The exposure weight *W* ^{X1}(*t*) corrects for potential time-varying confounding by *L* _{1}. The exposure weight *W* ^{X2}(*t*) corrects for potential time-varying confounding by *L* _{2}. The censoring weight *W* ^{D1}(*t*) corrects for potential selection bias due to censoring at dropout informed by *L* _{3}. The censoring weight *W* ^{D2}(*t*) corrects for potential selection bias due to censoring at death informed by *L* _{3}.

For *j* = 1, 2 and logit *p* = ln(*p*/(1 − *p*)), the numerator for the weights in equations (3) and (4) can be estimated from the following pooled logistic regression model:

Note that β_{j0}(*k* _{(} _{M} _{)}) are the time-specific intercepts with time-fixed confounders included in the pooled model. In addition, *Q* _{1}(*k* _{(} _{M} _{)}) = {X¯_{1}(*k* _{(} _{M} _{−1)}),X¯_{2}(*k* _{(} _{M} _{−1)}),*Z* _{1}},*Q* _{2}(*k* _{(} _{M} _{)}) = {X¯_{1}(*k* _{(} _{M} _{)}),X¯_{2}(*k* _{(} _{M} _{−1)}),*Z* _{2}},and β_{j1} is a vector of log hazard ratios for the effect of X¯ and *Z* _{j} on *X* _{j}(*k* _{(} _{M} _{)}). The denominator for the weights in equations (3) and (4) can be estimated from the following pooled logistic regression model:

Note that γ_{j0}(*k* _{(} _{M} _{)}) are the time-specific intercepts with time-fixed and time-varying confounders included in the pooled model. In addition,*U* _{1}(*k* _{(} _{M} _{)}) = {*X[Combining Macron]* _{1}(*k* _{(} _{M−1)}),*X[Combining Macron]* _{2}(*k* _{(} _{M−1)}),*[Combining Macron]L* _{1}(*k* _{(} _{M−1)})},*U* _{2}(*k* _{(} _{M} _{)}) = {*X[Combining Macron]* _{1}(*k* _{(} _{M} _{)}),*X[Combining Macron]* _{2}(*k* _{(} _{M−1)}),*[Combining Macron]L* _{2}(*k* _{(} _{M−1)})},and γ_{j1} is a vector of log hazard ratios for the effect of *X[Combining Macron]* and *[Combining Macron]L* _{j}(*k* _{(} _{M−1)}) on *X* _{j}(*k* _{(} _{M} _{)}).

Similarly, the numerator and denominator for the weights in equations (5) and (6) can be estimated from the following 2 pooled logistic regression models:

where θ_{j0}(*k* _{(} _{M} _{)}) and υ_{j0}(*k* _{(} _{M} _{)}) are time-specific intercepts with solely time-fixed and all common predictors, respectively, included in the pooled model. In addition, θ_{j1} and υ_{j1} are vectors of log hazard ratios for the effect of *X[Combining Macron]*(*k* _{(} _{M} _{)}) on *D* _{j}(*k* _{(} _{M} _{)}) adjusting for solely time-fixed and all common predictors of *T* and *C*, respectively. Finally, θ_{j2} and υ_{j2} are vectors of log hazard ratios for the effect of *Z* _{3} and *[Combining Macron]L* _{3}(*k* _{(} _{M} _{)}) on *D* _{j}(*k* _{(} _{M} _{)}), respectively, adjusting for *X[Combining Macron]*(*k* _{(} _{M} _{)}).

### Data Analysis

In the pooled logistic regression models as well as the joint marginal structural Cox model, *X* _{1} was alcohol consumption as drinks per week in the prior 6 months and *X* _{2} was drug injections per week in the prior 6 months. Drinks per week were categorized as 0–20 and 21–140, with the latter category termed “heavy” alcohol consumption. Drug injections per week were categorized as 0 and 1–39, with the latter category termed “any” drug injections. These categories were selected based on prior research^{19} and to ensure an ample number of HIV cases across the levels of the joint distribution of *X* _{1} and *X* _{2} while minimizing extreme weights and preserving meaningful comparisons.

In all models, predictors were included as either indicator variables for categories or a restricted quadratic spline with knots at the 5th, 35th, 65th, and 95th percentiles.^{21} In the pooled logistic regression models, the time-varying cofounders included in *L* _{1} and *L* _{2} were from the prior measurement time, while the time-varying common predictors included in *L* _{3} were from the most recent measurement time. Weights estimated from the pooled logistic models had a mean of 1.03 (standard deviation = 0.56), with a range of 0.11 to 17.99. These weights are shown over follow-up in Figure 2.

To examine additive interaction, results from the joint marginal structural model with a product term between alcohol consumption and injection drug use were used to estimate the relative excess risk due to interaction (RERI) as well as the corresponding confidence limits and *P* value based on the multivariate delta method. The RERI (= exp(α[Combining Circumflex Accent]_{1} + α[Combining Circumflex Accent]_{2} + α[Combining Circumflex Accent]_{3}) − exp(α[Combining Circumflex Accent]_{1}) − exp(α[Combining Circumflex Accent]_{2}) + 1) represents the increased hazard due to additive interaction where an RERI = 0 corresponds to the null hypothesis of no additive interaction.^{22}

Unadjusted and standard adjusted Cox models were provided for comparison. The standard adjusted model accounted for the same confounders as the weighted model, as well as concurrent values of time-varying confounders, to reflect practices in the literature.^{7},^{8},^{23} The product term between alcohol consumption and injection drug use would represent interaction in these models only in the absence of bias due to confounding, measurement, or selection with respect to both exposures. Conversely, the product term would represent effect measure modification if bias due to confounding, measurement, or selection was present for only 1 of the 2 exposures.^{2} Given that the product term between alcohol consumption and injection drug use in the joint marginal structural model was not statistically significant on the multiplicative scale, a joint marginal structural model as well as unadjusted and standard adjusted Cox models without the corresponding product terms were examined.

We assessed proportional hazards in the weighted Cox model by a joint robust Wald test on the product terms between the indicator for heavy drinks per week and time and between the indicator for any drug injections per week and time in the model without a product term between alcohol consumption and injection drug use. A similar joint Wald test was performed for log time. For the joint Wald test for the product terms with time, *P* = 0.0817; for the joint Wald test for the product terms with log time, *P* = 0.0315. The coefficients corresponding to the *P* value for the joint Wald test in the time model were also closer to the null than the coefficients in the log time model. This suggests nonproportionality in the weighted Cox model with respect to log time and, to a lesser degree, time.

To explore this nonproportionality, the effects of alcohol consumption and injection drug use were examined for follow-up times within the first 3 years of study enrollment and after 3 years, while including a product term between alcohol consumption and injection drug use in the joint weighted model. The duration of 3 years was selected because it was approximately the median time to infection among the HIV cases. To circumvent any occurrences of nonproportionality potentially due to selection bias innate to the hazard ratio,^{24} the cumulative incidence curves (ie, complement of survival functions) were compared by strata of alcohol consumption and injection drug use. All analyses were performed in SAS version 9.2.

## RESULTS

Table 1 describes the characteristics in the prior 6 months of participants at study entry and during follow-up. Among the 1525 participants, 28% were female. At study entry, the median age was 37 (quartiles = 32; 42) years and the median number of years of formal education was 10 (quartiles = 10; 12). During follow-up, 63% of the person-years occurred while participants were consuming any alcohol and 58% occurred while injecting drugs in the prior 6 months. Figure 3 shows alcohol consumption and injection drug use in the prior 6 months by years of follow-up. Alcohol consumption and any drug injecting in the prior 6 months seemed to decrease during follow-up.

Table 2 shows crude, adjusted, and weighted hazard ratios with 95% confidence intervals (CIs) for the joint effects of alcohol consumption and injection drug use in the prior 6 months on HIV acquisition averaged over the 8181 person-years of follow-up. These models each included a product term between alcohol consumption and injection drug use. In the weighted model, the effect of heavy drinking per week in the absence of any drug injections was a hazard ratio of 1.58 (95% CI = 0.67 to 3.73). The effect of any drug injections per week in the absence of heavy drinking was a hazard ratio of 1.78 (1.10 to 2.89). The effect of heavy drinking and any drug injections was a hazard ratio of 2.45 (1.45 to 4.12). The robust *P* value for the Wald test for the product term was 0.7620. The corresponding RERI was 0.08 (95% CI = −1.43 to 1.59), which corresponded to a *P* value of 0.9200. Analogous unadjusted hazard ratio estimates were 1.67 (0.72 to 3.84), 2.07 (1.31 to 3.28), and 3.36 (2.07 to 5.44), respectively, with a *P* value of 0.9474 for the Wald test for the product term. The corresponding standard adjusted hazard ratio estimates were 1.34 (0.57 to 3.13), 1.56 (0.87 to 2.77), and 2.13 (1.15 to 3.95), respectively, with a *P* value of 0.9597 for the Wald test for the product term.

Figure 4 shows the proportion who became HIV-positive by alcohol consumption and injection drug use for the prior 6 months based on the weighted Cox model. Proportions were estimated among men at the median age and years of education of the entire cohort at study entry. Participants who reported both injecting and heavy drinking had the highest proportion of HIV infections during follow-up. Those who reported either injecting or heavy drinking, but not both, had the second highest proportion of HIV infections. Those who did not report injecting or heavy drinking had the lowest proportion of HIV infections.

In the weighted Cox model without a product term between alcohol consumption and injection drug use, the effect of heavy drinking per week was a hazard ratio of 1.41 (0.98 to 2.01). The effect of any drug injections per week was a hazard ratio of 1.72 (1.12 to 2.64). Analogous unadjusted hazard ratio estimates were 1.63 (1.17 to 2.27) and 2.06 (1.38 to 3.07), respectively. The corresponding standard adjusted hazard ratio estimates were 1.36 (0.97 to 1.92) and 1.57 (0.94 to 2.62), respectively.

Table 3 shows the weighted hazard ratio with 95% CI for the joint effects of heavy drinking and any drug injections in the prior 6 months within the first 3 years of follow-up and after. This model included a product term between alcohol consumption and injection drug use. The joint effects were larger in the first 3 years of follow-up compared with between 3 and 20 years. However, the CI for follow-up times within the first 3 years was particularly imprecise.

## DISCUSSION

Among 1525 African American injection drug users in the ALIVE cohort followed for 8181 person-years, the effects of heavy drinking and any drug injections on HIV acquisition seemed not to depart from multiplicative or additive effects. We could not rule out interaction on either scale due to imprecision. If the results of the weighted analysis are corroborated by independent data sources, these results would suggest that alcohol-specific HIV risk reduction interventions should be targeted to both injectors and noninjectors and, similarly, injection-specific interventions should be targeted to both heavy and nonheavy drinkers.

However, based on the weighted model, the hazard ratio for the effect of heavy drinking and any drug injections seemed to be larger within the first 3 years of follow-up. This larger hazard ratio may be due to a waning effect of heavy drinking and injecting on HIV acquisition over time. Evidence against a waning effect was provided in Figure 4, which shows that participants who reported both injecting and heavy drinking consistently had the highest proportion of HIV infections. Alternatively, chance or selection bias due to a greater initial depletion of those susceptible to HIV infection among heavy drinkers and injectors may be the source of the observed attenuation with time.^{24} Selection bias due to a greater initial depletion of susceptible heavy drinkers and injectors would not be accounted for in the censoring weights, given that a measure of innate susceptibility to HIV infection was not available as a covariate when estimating the censoring weights.

The standard adjusted analysis provided hazard ratio estimates that were somewhat attenuated compared with the results from the weighted analysis, both analyses included product terms between alcohol consumption and injection drug use. Specifically, the excess hazard for the effect of heavy drinking and any drug injections from the weighted analysis was 22% (1–1.13/1.45) stronger than results from the standard adjusted analysis. These attenuated hazard ratio estimates may be due to the fact that standard adjustment can remove indirect effects mediated through time-varying confounders and induce selection bias.^{3},^{25} ^{–} ^{27}

For unbiased estimation of effects in the context of standard regression and marginal structural models, the assumptions of exchangeability, consistency, positivity, and correct model specification must be met.^{28} As noted by Hernán et al,^{3} as well as by VanderWeele,^{2} the aforementioned assumptions must also hold for unbiased estimation of joint effects. In the context of the joint marginal structural model, the exchangeability assumption states that there are no unmeasured confounders of the relationship between *X* _{1} and *T* as well as the relationship between *X* _{2} and *T*. In addition, informative censoring by unmeasured factors must be absent. Formally, *f*[*X*(*k* _{(} _{M} _{)})|*X[Combining Macron]*(*k* _{(} _{M} _{−1)}), *[Combining Macron]L*(*k* _{(} _{M} _{−1)})] = *f*[*X*(*k* _{(} _{M} _{)})|*X[Combining Macron]*(*k* _{(} _{M} _{−1)}), *[Combining Macron]L*(*k* _{(} _{M} _{−1)}), *T* ^{x[Combining Macron]}] and *f*[*T*|*[Combining Macron]L*] = *f*[*T*|*[Combining Macron]L*, *C*]. For the ALIVE example, unmeasured factors may exist, given that the indicators used for high-risk behaviors may not fully capture the association between such behaviors and HIV acquisition or censoring among participants.^{19} Measurement error regarding covariates of interest poses an additional threat to exchangeability as does the observed limited number of HIV infections across strata of the joint distribution of alcohol consumption and injection drug use.

Consistency states that *TX[Combining Macron]* = *Tx[Combining Macron],p* if *X[Combining Macron]* = *x[Combining Macron]* for all *p*, where *p* is the mechanism by which the level of *X* is determined. As such, the observed failure time is the same as the failure time that would have occurred if *X* were assigned. It also implies that the mechanism by which the level of *X* is determined does not influence the effect of *X[Combining Macron]* on *TX[Combining Macron]*. In other words, among ALIVE participants, the time to HIV acquisition for the observed level of alcohol consumption and number of drug injections is the time to HIV acquisition that would have occurred if the level of alcohol consumption and drug injections were assigned. In addition, the effect of heavy drinking or injection drug use does not depend on drinking or injection pattern (eg, binge drinking or frequent injecting).

Positivity requires that *f*[*X*(*k* _{(} _{M} _{)})|*X[Combining Macron]*(*k* _{(} _{M} _{−1)}), *[Combining Macron]L*(*k* _{(} _{M} _{−1)}), *[Combining Macron]D*(*k* _{(} _{M} _{)}) = 0, *T ≥ k* _{(M} _{)}] > 0 for all *x* ∈ *X* and *P*(*D*(*k* _{(} _{M} _{)}) = 0|*X[Combining Macron]*(*k* _{(} _{M} _{)}), *[Combining Macron]L*(*k* _{(} _{M} _{)}), *[Combining Macron]D*(*k* _{(} _{M} _{−1)}) = 0, *T* ≥ *k* _{(} _{M} _{)}) > 0 for all *k* _{(} _{M} _{)} ≤ max(*y* _{i}). Therefore, every level of alcohol consumption as well as injection drug use must be possible in each subset of the ALIVE cohort defined by the observed histories of alcohol consumption, injection drug use, and measured confounders. Similarly, there must be participants who are not censored at each observed follow-up time within each subset of the cohort defined by the observed histories of alcohol consumption, injection drug use, and measured common predictors of censoring and HIV acquisition time. In the present setting, we did not expect deterministic departures from positivity, and we observed no evidence of random departures from positivity. The latter is based on a lack of extreme weights.^{28},^{29}

Correct model specification implies that appropriate functional forms are used in the pooled logistic and final weighted models such that *f*[*X*(*k* _{(} _{M} _{)})|*X[Combining Macron]*(*k* _{(} _{M} _{−1)}),*[Combining Macron]L*(*k* _{(} _{M} _{−1)})] = *f*[*X*(*k* _{(m)})|*g* _{0}(*X[Combining Macron]*(*k* _{(} _{M} _{−1)})),*g* _{1}(*[Combining Macron]L*(*k* _{(} _{M} _{−1)}))] and *P*(*D*(*k* _{(} _{M} _{)}) = 0|*X[Combining Macron]*(*k* _{(} _{M} _{)}),*[Combining Macron]L*(*k* _{(} _{M} _{)}),*[Combining Macron]D*(*k* _{(} _{M} _{−1)}) = 0,*T* ≥ *k* _{(} _{M} _{)}) = *P*(*D*(*k* _{(} _{M} _{)}) = 0|*g* _{2}(*X[Combining Macron]*(*k* _{(} _{M} _{)})),*g* _{3}(*[Combining Macron]L*(*k* _{(} _{M} _{)})),*[Combining Macron]D*(*k* _{(} _{M} _{−1)}) = 0,*T* ≥ *k* _{(} _{M} _{)}),where *g* _{0}(*X[Combining Macron]*(*k* _{(} _{M} _{−1)})),*g* _{1}(*[Combining Macron]L*(*k* _{(} _{M} _{−1)})),*g* _{2}(*X[Combining Macron]*(*k* _{(} _{M} _{)})),and *g* _{3}(*[Combining Macron]L*(*k* _{(} _{M} _{)})) are coarsening functions of *X[Combining Macron]*(*k* _{(} _{M} _{−1)}),*[Combining Macron]L*(*k* _{(} _{M} _{−1)}),*X[Combining Macron]*(*k* _{(} _{M} _{)}),and *[Combining Macron]L*(*k* _{(} _{M} _{)}), respectively. The lack of extreme weights in the ALIVE cohort also provides some evidence against model misspecification.^{28},^{29} In addition, positivity and correct model specification were facilitated by careful consideration of all categorizations, particularly for alcohol consumption and injection drug use. Furthermore, indicators for categories and splines in the case of continuous predictors were used in the modeling process.

The aforementioned definitions and assumptions imply that unbiased assessment of interaction requires no unmeasured confounders of the relationship between both exposures and the outcome. In contrast, assessment of effect measure modification requires only no unmeasured confounders of the relationship between the primary exposure and outcome.

The methods presented in this paper pertain to the assessment of interaction using marginal structural models. Such methods are necessary when both exposures are potentially subject to time-varying confounding affected by prior exposure, which can occur in both chronic and infectious disease epidemiology.^{7} ^{–} ^{10} If effect measure modification of a time-varying exposure by a time-fixed secondary covariate is of interest and that exposure is subject to time-varying confounding affected by prior exposure, then a marginal structural model is needed. This model can be a standard marginal structural model with a single set of exposure weights, and a product term between the exposure and the potential time-fixed effect measure modifier can be employed. However, if effect measure modification of a time-varying exposure by a time-varying covariate is of interest and that exposure is subject to time-varying confounding affected by prior exposure, then methods such as G-estimation of structural nested models,^{30} history-adjusted marginal structural models,^{31} or dynamic marginal structural models^{32} are needed.

Important to note is that a standard marginal structural model can also be considered a joint marginal structural model, with censoring playing the role of the second exposure. However, typically, we are interested only in the effect of the exposure while precluding the censoring. Therefore, rather than estimate the entire joint distribution for the exposure and censoring, we estimate the effect of the exposure with censoring set to never.

In conclusion, estimation of joint effects is frequently of interest in epidemiologic research. This paper describes and demonstrates joint marginal structural Cox proportional hazards models in the context of time-varying exposures potentially subject to time-varying confounding affected by prior exposure, with the aim to encourage the use of these valuable techniques.

### ACKNOWLEDGMENTS

We thank Jacquie Astemborski for assistance with the ALIVE data and Petra Sander and Deborah Anderson for editorial assistance, as well as the ALIVE study staff and participants.