A Comparison of the Polytomous Logistic Regression and Joint Cox Proportional Hazards Models for Evaluating Multiple Disease Subtypes in Prospective Cohort Studies

There is growing recognition of the importance of distinguishing histologic or molecular subtypes in etiologic studies of cancer. As a result, comparing exposure effects across multiple subtypes of a given cancer is often of primary interest. For example, Schroeder and colleagues (1) examined whether the association of cigarette smoking with bladder cancer differed according to the presence of p53 mutations in the tumor specimen. The conventional statistical approach for analyzing case–control studies with more than one outcome subtype is to fit a polytomous logistic regression model (ref. 2; for simplicity, referred to as polytomous model). In the polytomous model, the associations of genetic or environmental factors with all disease subtypes are estimated simultaneously and inferences comparing associations across subtypes can be conducted (3–7).

Polytomous models have also been used in prospective cohort studies to compare the associations between exposure variables and occurrence of multiple disease subtypes (8–13). In the Iowa Women's Health Cohort Study (8), for example, polytomous models were applied to compare the effects of specific risk factors on the development of different types of breast cancer defined according to the presence or absence of estrogen (ER^+/−) and progesterone receptor (PR^+/−) expression status. However, the performance of the polytomous model, that is, the validity, accuracy, and efficiency of the model, for prospective cohort studies has not been adequately examined in the epidemiologic and statistical literature. A major drawback of applying the polytomous model in prospective cohort studies is that participants are classified only according to whether or not the specific subtype of disease occurred—data on the actual time to development of the different subtypes are not used. An additional limitation of the polytomous method relates to the nature of the control group. For example, to evaluate the association of a specific exposure with ER⁺ breast cancer, the control group includes subjects who did not develop any type of breast cancer during the follow-up period, excluding those who developed ER⁻ breast cancer. Therefore, to approximate the ORs associated with the risk of ER⁺ breast cancer, it is necessary to assume that ER⁻ breast cancer is rare (event rate < 10%) and vice versa. However, when the event is not rare, polytomous models could lead to severely biased estimates of exposure–disease associations.

As an alternative in a prospective cohort study setting, we propose to use a joint Cox proportional hazards model (for simplicity, referred to as joint Cox model) to examine associations of risk factors with multiple disease subtypes. The model accommodates time-to-event data, different baseline hazard functions of each disease subtype, and direct comparison of associations with outcome subtypes. This model was originally proposed for the analysis of competing risk events (14–16). Chatterjee and colleagues (17) extended it to model a disease with a large number of possible subtypes defined by multiple disease characteristics and to allow for the possibility that some cases cannot be categorized into subtypes because of missing characteristics. A 2-stage Cox model and an estimating equation approach were proposed to estimate exposure–disease subtype associations. However, polytomous models remain the most commonly used method for the analysis of multiple disease subtypes because of lack of recognition of potential problems associated with the model and limited familiarity with the formulation and implementation of the 2-stage Cox model. In this article, we formally evaluate the performance of the polytomous model and compare it with the Cox model using both mathematical and simulation-based approaches. For the Cox model, we consider a special case of Chatterjee's model for a limited number of disease subtypes and propose a joint Cox model, which can be directly estimated using existing software.

We then apply the methods to a real data set from the NIH-AARP Diet and Health Study to examine heterogeneity in the association of body size, defined by body mass index (BMI), with breast cancer subtypes characterized by the tumor's ER/PR status (18). In a previous analysis using the polytomous model, BMI was shown to have a significant association only with ER⁺/PR⁺ breast cancer but not with the other subtypes. However, statistical power was limited because of a relatively small number of cases, particularly for subtypes other than ER⁺/PR⁺; in addition, the validity of the polytomous model was not examined. This article presents an extension of the original analysis with 5 more years of follow-up.

Consider a disease with only 2 subtypes, denoted as event 1 and event 2, which are competing risks (e.g., ER⁺ and ER⁻ breast cancer). The hazard functions for time to event 1 and 2, T₁ and T₂, are λ₁(t) and λ₂(t), respectively. Individuals who develop one event are considered censored at the time of the other event. Assume that the hazard for each type of event is affected by a set of common risk factors, X, according to the following Cox proportional hazards model:

where λ_0k(t) is the baseline hazard function and β_k denotes the log HR associated with a unit increase in X (e.g., 1 kg/m² increase in BMI) for the development of the k-th event type (k = 1, 2). By fitting a separate Cox model for each event, one can estimate β₁ and β₂ but not evaluate if β₁ = β₂. The polytomous logistic regression model transforms the time-to-event data to a categorical outcome to examine β₁ = β₂.

Let D denote a polytomous outcome where D = 0 if the subject remains free of disease during follow up; D = k if event k occurs for k = 1, 2. The polytomous model, also known as the multinomial logistic regression model, can be expressed as:

The ratio |$P(D = k|X)/P(D = 0|X)$| in Eq. (2) is referred to as the “odds-like” term for disease type D = k in comparison to the reference category of no disease (D = 0). Note that this odds-like term for event k is not the same as the odds of event k, |$P(D = k|X)/\{1 - P(D = k|X)\}$|⁠, as the odds denominator is the sum of the probability of no event and the probability of all events other than k. It follows from (2) that |$e^{\gamma _k} \equiv {\rm OR}_k^*$|⁠, the ratio of the odds-like terms for event k associated with a unit increase in X, which is again different from the usual OR.

The polytomous model tests whether |$\gamma _1 = \gamma _2$| or equivalently |${\rm OR}_1^* = {\rm OR}_2^*$|⁠. It is commonly assumed that the coefficients obtained using the polytomous models are similar to the corresponding coefficients in the Cox proportional hazards model, that is, |$\gamma _k \approx \beta _k$|⁠. However, as we have shown in Appendix I, this is only true when both events are rare and their baseline hazards functions are proportional. For example, |$\gamma _1 >\beta _1$| when |$\beta _1 >1$| and the event of interest is not rare. Nonetheless, a test of |$\gamma _1 = \gamma _2$| can provide a valid test for β₁ = β₂ if the baseline hazards are proportional. Appendix I shows that as long as the 2 events have proportional baseline hazards functions, |$\gamma _1 - \gamma _2 = \beta _1 - \beta _2$| is true regardless of whether the event rates are rare or not. Furthermore, if the subjects who are lost to follow-up without developing any events are included in the non-event category, the polytomous model continues to provide a valid test for β₁ = β₂ as long as the censoring is uninformative of future events.

If the baseline hazard functions of the different events are not proportional over time, |$\gamma _1 - \gamma _2 = \beta _1 - \beta _2$| does not necessarily hold. It is possible that, for example, the baseline hazards function increases with age for one event but decreases with age for the other event. In this case, an alternative approach for comparing the effects of X on the risks of the different events is needed.

Chatterjee extended the competing risk framework to a 2-stage Cox model to model disease subtypes, in which the subtypes are defined by multiple disease characteristics some of which may be missing for some of the cases (17, 19). However, for the reasons cited above, polytomous models are still the most commonly used methods to compare exposure effects across multiple disease subtypes. Here, we consider a special case of Chatterjee's model in which the disease is categorized into a limited number of subtypes to have sufficient number of cases in each category and no further substructure among the disease subtypes is assumed. We propose a joint Cox proportional hazards modeling approach with an estimating procedure that differs slightly from Chatterjee's model (17) so that standard procedures in commonly used statistical software can be directly applied to evaluate the exposure and disease associations for multiple disease subtypes in a prospective cohort setting.

In a joint Cox regression model, the time-to-event 1 and event 2 are simultaneously modeled. Specifically, each subject contributes 2 potential event times, one for each event. Suppose the follow-up time for a subject is τ. (Appendix II shows the definition of the 2 events.) Marginally, the hazards function for event k is defined the same as the separate Cox model for the event [Eq. (1)]; however, the parameters of interest β₁ and β₂ are estimated simultaneously by stratifying on event type and using a robust variance estimator to account for the correlation between the survival times of the 2 events (20). Specifically, Eq. (1) is reparameterized as follows:

where I_{k = 2} equals to 1 if the event time is for event 2 and 0 otherwise so that the HRs associated with X for the development of events 1 and 2 are |$e^{\beta ^*}$| and |$e^{\beta ^* + \theta}$|⁠, respectively. Note that Eq. (3) is equivalent to Eq. (1) so that β₁ = β* and |$\beta _2 = \beta ^* + \theta$|⁠. Thus, an estimate of θ provides an estimate of β₁ − β₂ and a test of β₁ = β₂ is equivalent to a test of θ = 0, which can be conducted with a Wald test.

We initiated a simulation study to (i) evaluate the polytomous model estimates of the odds-like ratio as an estimate of OR and HR; (ii) evaluate the performance of the polytomous model for comparing the exposure–disease associations between 2 competing events; (iii) evaluate the performance of the joint Cox regression model for estimating and comparing the exposure–disease associations for both events; and (iv) compare the performances of the joint Cox model and the polytomous model when the baseline hazards functions for the 2 events are nonproportional.

The simulation study was designed as follows:

1. Sample size: N = 2,000; for simplicity, we considered X to be a single binary covariate with P(X = 1) = 0.50. For example, X can be a biomarker that was originally continuous but dichotomized at its median level.

2. Time to event: Times to events 1 and 2 were generated from exponential distributions based on Eq. (1) while assuming various magnitudes of baseline event rate and strengths of exposure–disease associations (Table 1). Events 1 and 2 have hazards functions that are proportional to each other in this case. To evaluate the effect of nonproportional baseline hazards, time to the second event was generated from a Weibull distribution with an increasing hazard function.

3. Censoring due to competing events and follow-up length of 10 years was generated; censoring due to loss to follow-up was also considered by generating a random censoring time from an independent exponential distribution (see Appendix II for detail).

Eight scenarios were generated (Table 1): scenarios (a)–(d) considered the cases when times to event 1 and event 2 have proportional hazards functions whereas scenarios(e)–(h) considered the cases when the hazards function for times to event 1 and 2 are not proportional to each other. For scenarios with the proportional hazards functions, (a) and (b) considered rare events whereas (c) and (d) considered non-rare events; (a) and (c) considered X has an effect on event 1 but not on event 2 whereas (b) and (d) considered X has the same effects on both events. The set with nonproportional hazards functions were designed similarly, that is, (e) and (f) for rare events and (g) and (h) for non-rare events and (e) and (g) considered differential effects and (f) and (h) considered same effects for events 1 and 2.

The simulation procedure was repeated 500 times. The performance of each model was summarized according to relative bias, coverage probability of the confidence interval (CI) for estimating OR₂, |$e^{\beta _2}$|⁠, |$e^{\beta _1 - \beta _2}$| and empirical power or empirical type I error for testing equality of exposure–disease associations (for simplicity, we omit results for event 1). A coverage probability close to 95% indicates that the variance estimator is consistent, provided that the parameter estimator is unbiased. The empirical power indicates how well the model is able to detect a true exposure effect or differential exposure effects and the closeness of the empirical type I error to the target level of 5% is a measure of the validity of the statistical test.

We used the polytomous model and the joint Cox model to formally compare the association of BMI with breast cancer by ER and PR status using data from the NIH-AARP Diet and Health Study. The NIH-AARP study was established in 1995 when a baseline screening questionnaire that queried information on medical history and lifestyle characteristics was returned by 566,398 AARP (formerly known as the American Association of Retired Persons) members. Participants were aged 50 to 71 years and resided in 1 of 6 U.S. states (CA, FL, LA, NJ, NC, and PA) or 2 metropolitan areas (Atlanta, GA, and Detroit, MI), selected because they were known to have high-quality cancer registries and a large AARP membership (18).

A total of 54,629 women remained in the analytic cohort after exclusion of premenopausal women, women with previous history of cancer as well as current users of hormone therapy and women residing in FL or PA because these 2 states have substantial proportion of participants with missing ER/PR data. In this cohort, 1,492 cases were diagnosed with breast cancer with known ER/PR status obtained from cancer registry records, in which 246 were classified as ER⁻/PR⁻ breast cancer cases (16.5%), 231 ER⁺/PR⁻ cases (15.5%), 18 ER⁻/PR⁺ cases (0.01%), and 997 ER⁺/PR⁺ cases (66.8%). BMI (kg/m²) was derived from self-reported height and weight on the baseline survey. In statistical models, BMI was expressed as categorical variables (<25.0, 25.0–29.9, and ≥30.0 kg/m²). Covariates included in the multivariable models were consistent with those reported in the prior analysis by Ahn and colleagues (18), namely age, race or ethnic, family history of breast cancer, level of education, age at menarche, age at menopause, age at first birth, parity, smoking status, physical activity, fat intake, alcohol consumption, oophorectomy, and height. The HRs for the association of BMI with the 4 breast cancer subtypes were estimated simultaneously and the difference in HRs were assessed. The analysis was conducted using SAS 9.2 (SAS Institute Inc.).

When hazards functions for events 1 and 2 are proportional and both outcomes are rare [scenarios (a and (b)], the odds-like ratio estimated in the polytomous model is unbiased for the true OR as well as the true HR as indicated by the small relative bias (<2%; Table 2). The variance of the estimator is also unbiased as the coverage probability of CIs is close to 95%. When the events are more common, the odds-like ratio obtained from the polytomous model can be a biased estimate for the true ORs and HRs. The level of the bias increases with incidence rates and the magnitude of the exposure–disease association: when both event rates are above 10%, the relative bias for the odds-like ratio as an estimate for HR₂ is 4.6% when HR₂ = 1 [scenario (c)] and 6.3% when HR₂ = 1.33 [scenario (d)]. Furthermore, the odds-like ratio can be a much poorer estimate of the HR than of the OR: for example, in scenario (d), the relative bias is 2.2% for OR₂ and 6.3% for HR₂; the estimates are worse with nonproportional hazards functions: in scenario (h), the relative bias is 2.6% for OR₂ and 12.3% for HR₂ and the CI coverage probability is 94.4% for OR₂ and 79.4% for HR₂.

In contrast, when the hazards functions for events 1 and 2 are proportional [scenarios (a)–(d)], the estimate of the ratio of the 2 HRs, estimated using the ratio between the 2 odds-like ratios obtained from the polytomous model, is accurate regardless of the disease rate: the relative bias is within 3%, similar to that found in the joint Cox model (see below). This is consistent with the results from the mathematical approach (Appendix I). The variance estimates from the simulation studies are also accurate so that CI coverage probability and type I error rate are adequate. Furthermore, statistical power using this approach is comparable with that of the joint Cox model. Therefore, the use of the polytomous model to compare exposure–disease subtype associations is valid as long as the hazard functions for the disease subtypes are proportional. When the 2 events do not have proportional hazards functions, the polytomous model results in a larger bias [relative bias 4%–9% for scenarios (e)–(h)] in estimating the ratio of HRs and the bias is higher when the disease rate is higher: the relative bias is over 10 times higher than that for the joint Cox model: 8.8% versus 0.6% in scenario (h), although the power and type I error rate for testing β₁ = β₂ are comparable with what was obtained for the joint Cox model.

Table 3 shows that the joint Cox model provides an unbiased point and interval estimates of the true HRs for each event, regardless of the event rates and the proportionality of the hazards functions—the relative bias is small and the coverage probability is always close to 95%. The same is true for the estimation of the ratio of the HRs.

As a sensitivity analysis, we considered cases when one event is rare and the other is not. The same conclusion remains for both models as if both events are common; a less common exposure, that is, P(X = 1) < 0.5 as well as a continuous X were also considered and the results remained similar (not shown).

We first applied the polytomous model and found that compared to leaner women (BMI < 25.0 kg/m²), those who were overweight (25.0–29.9 kg/m²) or obese (≥30.0 kg/m²) had an increased risk of ER⁺/PR⁺ breast tumors (OR, 1.59 and 2.10, respectively) with a statistically significant linear trend (P_trend < 0.0001; Table 4). In contrast, overweight (OR, 1.42) but not obesity was associated with a higher risk of ER⁻/PR⁻ breast tumors without a significant linear trend (P = 0.28). BMI was not associated with risk of ER⁺/PR⁻ breast cancer. As all the breast cancer subtypes are rare, the odds-like ratio provides a reasonable estimate of the ORs and the HRs. We then applied the joint Cox model that also showed that both overweight (HR, 1.55; 95% CI, 1.30–1.85) and obese (HR, 1.76; 95% CI, 1.47–2.10) are significantly associated with ER⁺/PR⁺ cancer and a significant linear trend over BMI groups only among ER⁺/PR⁺ cancers (P < 0.0001) but not for other types (ER⁻/PR⁻, P_trend = 0.69 for ER⁻/PR⁻ and 0.60 ER⁻/PR⁺; Table 4). Among the ER⁻/PR⁻ breast tumors, overweight (HR, 1.43) but not obesity was associated with a higher risk of this subtype as in the polytomous model. Although both the joint Cox model and the polytomous model allow for the other covariates in the model to have the same effects across the different disease subtypes, implementation of such polytomous models is not as straightforward as the joint Cox model using standard procedures in SAS. Therefore, in this analysis, the joint Cox model but not the polytomous model assumed other covariates to have the same effects across the disease subtypes. Consequently, despite the small number of cases, the joint Cox model was able to yield estimates of HRs for BMI associated with ER⁻/PR⁺ breast cancer, albeit with wide CI, whereas the polytomous model failed to converge and produce estimates of ORs. However, one should be cautious about the interpretation of the result for the ER⁻/PR⁺ category because of its limited size.

We next formally evaluated whether the magnitude of the association of BMI with disease differed across subtypes. As the validity of the polytomous model for this analysis depends on the proportional hazards assumption, we examined whether the 4 subtypes have proportional baseline hazards. Graphical approaches as well as statistical tests can be used to assess proportionality. However, the statistical tests are influenced by sample size and cannot provide insights on how the data deviate from proportionality. Therefore, we used a graphical approach based the scaled Schoenfeld residuals (21): nonproportionality of baseline hazards would be suggested by a systematic decreasing or increasing trend in a Lowess smoothed curve over Schoenfeld residuals. Figure 1 indicates a decreasing trend between ER⁻/PR⁻ and ER⁺/PR⁺, and Fig. 2 suggests an increasing trend between ER⁺/PR⁻ and ER⁺/PR⁺. Because the baseline hazards function of ER⁺/PR⁺ appeared to be nonproportional to those of the ER⁻/PR⁻ and ER⁺/PR⁻, the joint Cox model was used to estimate the ratio of the HRs: the HR associated with obesity for ER⁻/PR⁻ is significantly lower than that for ER⁺/PR⁺ (relative HR, 0.59; 95% CI, 0.40–0.88); the relative HR associated with each increase of categories of BMI is 0.78 (0.65–0.94), that is, the linear trend for ER⁺/PR⁺ is 22% stronger than that for ER⁻/PR⁻.

Polytomous models are currently the most commonly used method to compare exposure–disease associations across multiple subtypes of disease; however, the model has several potential limitations when applied to prospective epidemiologic studies. We have provided evidence that supports our recommendation that the joint Cox regression model should be used for time-to-event data. This approach has been under-used for the comparison of exposure effects on different subtypes of disease in epidemiologic studies because lack of awareness of potential problems associated with the polytomous models and lack of familiarity of the data augmentation and implementation of the joint Cox model using standard statistical software. We showed through mathematical and extensive simulation studies that both methods conduct well in comparing exposure and disease subtype associations under the proportionality assumption. However, when the proportionality in hazards functions across disease subtypes is violated, polytomous models can yield biased estimate of differential effects, particularly when at least one event rate exceeds 10%. Joint Cox regression models, in contrast, are a valid approach regardless of disease rates and whether or not proportionality across subtypes satisfies. As the proportional hazards assumption can often be of question, the joint Cox model provides a more robust solution for examining heterogeneity in exposure–cancer subtype associations. Joint Cox regression models are also easy to implement with commonly used statistical software (see Appendix III for codes in SAS and S-plus or R). We note that when adjusting for other covariates or potential confounders in a polytomous model, standard procedures in SAS typically assume that the covariates have differential effects on different subtypes of disease whereas MultinomRob (22) in R and mlogit in STATA allow some or all of the covariates to be constrained to have the same effect across disease subtypes. In contrast, a joint Cox model can be easily implemented using standard procedures in commonly used statistical software, regardless of whether the effects of the covariates are assumed to be the same or different across disease subtypes.

Other alternative approaches to the polytomous models have also been used. For example, Prentice and colleagues used a stratified Cox regression model to compare the effects of a dietary intervention on multiple subtypes of breast cancer in which independence across competing event types was assumed (23, 24). However, a woman with multiple risk factors might be more likely to develop any subtype of breast cancer; therefore, risk of these 2 subtypes could be positively correlated. Consequently, the estimated coefficients for the same covariate could be correlated as well. Failure to take into account these correlations can overestimate the variation of the difference in coefficients, resulting in a loss of power for detecting heterogeneity in disease associations.

When a rare disease such as cancer is considered, a more efficient prospective study design such as a nested case–control study or a case–cohort study is often used. The extension of the joint Cox model to these study designs for the analysis of disease with multiple subtypes would be of great interest for future research.

Although not our main focus, this study also served to further show that BMI has a significant positive association with ER⁺/PR⁺ breast cancer, and the strength of the association for ER⁺/PR⁺ breast tumors is significantly greater than that for other breast cancer subtypes. The current study represents an update of the prior data from the NIH-AARP study and includes 4 times more breast cancer cases. Because the rare disease assumption applied in this study and the departure from proportionality of the baseline hazards was not substantial, results from the polytomous models did not differ materially from that of the joint Cox model.

In conclusion, to our knowledge, this article is the first to formally evaluate and compare the polytomous models and the joint Cox models within a prospective setting. We hope that the result of this study encourage greater use of the joint Cox model for the analysis of cancer subtypes in prospective cohort studies.

Without loss of generality, consider a disease with only 2 subtypes, denoted as event 1 and event 2, which are competing risks for one another. The hazards functions for times to event 1 and 2, T₁ and T₂, are |$\lambda _1 (t)$| and |$\lambda _2 (t)$|⁠, respectively. Assume that the hazard for each type of event is affected by a set of common risk factors, X, according to the following Cox proportional hazards model:

where |$\lambda _{0k} (t)$| is the baseline hazard function for the k-th event type and β_k denotes the log HR associated with a unit increase in X for the development of the k-th event type (k = 1, 2). Let D denote a polytomous outcome where D = 0 if the subject remains free of disease during follow-up; D = k if event k occurs for k = 1, 2. The polytomous logistic regression model, also known as the multinomial logistic regression model, can be expressed as:

Note that

and

and

Assume the baseline hazard functions for event 1 and event 2 are proportional to each other, that is, |$\lambda _{02} (t) = \delta \lambda _{01} (t)$|⁠. Then, we have

and

where |$\Lambda _{01} (\tau)$| and |$\Lambda _{02} (\tau)$| are the cumulative hazards functions by time t for event 1 and 2, respectively, |$\Lambda _{02} (\tau) = \delta \Lambda _{01} (\tau)$| and

Thus,

and

.

Therefore, |$\gamma _1 \ne \beta _1$| and |$\gamma _2 \ne \beta _2$|⁠, however, |$e^{\gamma _1 - \gamma _2} = e^{\beta _1 - \beta _2}$|⁠. That is, a test of |$\gamma _1 = \gamma _2$| is equivalent to a test of |$\beta _1 = \beta _2$|⁠.

When both event rates are rare, |$\Lambda _{01} (\tau)$| and |$\Lambda _{02} (\tau)$| are small, so that |$e^{\gamma _1} \approx e^{\beta _1}$| and |$e^{\gamma _2} \approx e^{\beta _2}$|⁠. Thus only under the rare event assumption we can use γ₁ to estimate β₁ and γ₂ to estimate β₂.

In the presence of independent censoring (e.g., censoring due to early loss to follow-up), we denote the censoring time by C so that

and

Then,

and

so that

Therefore, a polytomous model can still provide a valid estimate of ratio of HRs between event 1 and event 2 in the presence of early loss to follow-up as long as both events have hazards functions, which are proportional to each other.

Without censoring due to early loss of follow-up, the probability of no events by the end of follow-up τ is |$P(D = 0|x) = e^{- \Lambda _1 (\tau) - \Lambda _2 (\tau)}$|⁠, where |$\Lambda _1 (\tau)$| and |$\Lambda _2 (\tau)$| are the cumulative hazards functions up to time τ for events 1 and 2, respectively. With the presence of censoring due to early loss of follow-up, the same probability becomes

where |$\Lambda _c (\tau)$| is the cumulative hazards function for an independent censoring time. We limit the increase in probability of no event with early loss of follow-up to be not more than 5% then the same probability without early loss of follow-up, which implies

For exponential times to event 1 and 2 and time to censoring, the above equation reduces to

where we numerically solve for γ_c with a given γ₁ and γ₂ based on the event rates we assumed for each scenario; for exponential time to event 1 and Weibull for time to event 2, the rate is calculated similarly.

Definition of the joint survival outcome under 3 scenarios

SAS code for the joint Cox model:

• proc phreg covs(aggregate);

• model time*status(0) = exposure expotype;

• strata event2;

• id subject;

• expotype = exposure*event2;

• run;

R or Splus code for the joint Cox model:

• coxph(Surv(time,status)∼exposure*strata(event2)+cluster(subject),robust = T)

No potential conflicts of interest were disclosed.

Conception and design: X. Xue, A.R. Hollenbeck, M.J. Gunter

Development of methodology: X. Xue, M.Y. Kim, M.J. Gunter

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Xue, Y. Park, A.R. Hollenbeck, H.D. Strickler

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): X. Xue, H.D. Strickler, M.J. Gunter

Writing, review, and/or revision of the manuscript: X. Xue, M.Y. Kim, M.M. Gaudet, Y. Park, M. Heo, A.R. Hollenbeck, H.D. Strickler, M.J. Gunter

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Park, A.R. Hollenbeck, H.D. Strickler

The authors thank the referees for their excellent comments, which helped greatly improve the article.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.