Abstract
Background: Polytomous logistic regression models are commonly used in case–control studies of cancer to directly compare the risks associated with an exposure variable across multiple cancer subtypes. However, the validity, accuracy, and efficiency of this approach for prospective cohort studies have not been formally evaluated.
Methods: We investigated the performance of the polytomous logistic regression model and compared it with an alternative approach based on a joint Cox proportional hazards model using simulation studies. We then applied both methods to a prospective cohort study to assess whether the association of breast cancer with body size differs according to estrogen and progesterone receptor–defined subtypes.
Results: Our simulations showed that the polytomous logistic regression model but not the joint Cox regression model yielded biased results in comparing exposure and disease subtype associations when the baseline hazards for different disease subtypes are nonproportional. For this reason, an analysis of a real data set was based on the joint Cox proportional hazards model and showed that body size has a significantly greater association with estrogen and progesteronepositive breast cancer than with other subtypes.
Conclusions: Because of the limitations of the polytomous logistic regression model for the comparison of exposure–disease associations across disease subtypes, the joint Cox proportional hazards model is recommended over the polytomous logistic regression model in prospective cohort studies.
Impact: The article will promote the use of the joint Cox model in a prospective cohort study. Examples of SAS and Splus programming codes are provided to facilitate use by nonstatisticians. Cancer Epidemiol Biomarkers Prev; 22(2); 275–85. ©2013 AACR.
Introduction
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 followup 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 timetoevent 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 2stage 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 2stage 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 simulationbased 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 NIHAARP 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 followup.
Materials and Methods
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_{1} and T_{2}, are λ_{1}(t) and λ_{2}(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^{2} increase in BMI) for the development of the kth event type (k = 1, 2). By fitting a separate Cox model for each event, one can estimate β_{1} and β_{2} but not evaluate if β_{1} = β_{2}. The polytomous logistic regression model transforms the timetoevent data to a categorical outcome to examine β_{1} = β_{2}.
Polytomous logistic regression model
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 in Eq. (2) is referred to as the “oddslike” term for disease type D = k in comparison to the reference category of no disease (D = 0). Note that this oddslike term for event k is not the same as the odds of event k, , 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 , the ratio of the oddslike terms for event k associated with a unit increase in X, which is again different from the usual OR.
The polytomous model tests whether or equivalently . 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, . 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, when and the event of interest is not rare. Nonetheless, a test of can provide a valid test for β_{1} = β_{2} if the baseline hazards are proportional. Appendix I shows that as long as the 2 events have proportional baseline hazards functions, is true regardless of whether the event rates are rare or not. Furthermore, if the subjects who are lost to followup without developing any events are included in the nonevent category, the polytomous model continues to provide a valid test for β_{1} = β_{2} as long as the censoring is uninformative of future events.
If the baseline hazard functions of the different events are not proportional over time, 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 2stage 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.
Joint Cox proportional hazards model
In a joint Cox regression model, the timetoevent 1 and event 2 are simultaneously modeled. Specifically, each subject contributes 2 potential event times, one for each event. Suppose the followup 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 β_{1} and β_{2} 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 and , respectively. Note that Eq. (3) is equivalent to Eq. (1) so that β_{1} = β* and . Thus, an estimate of θ provides an estimate of β_{1} − β_{2} and a test of β_{1} = β_{2} is equivalent to a test of θ = 0, which can be conducted with a Wald test.
Simulation study
We initiated a simulation study to (i) evaluate the polytomous model estimates of the oddslike 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:
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.
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.
Censoring due to competing events and followup length of 10 years was generated; censoring due to loss to followup 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 nonrare 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 nonrare 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_{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.
Comparison of the polytomous and joint Cox model for estimating the association of BMI with different breast cancer subtypes using data from the NIHAARP Diet and Health Study
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 NIHAARP Diet and Health Study. The NIHAARP 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 highquality 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^{2}) was derived from selfreported 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^{2}). 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.).
Results
Simulation study
When hazards functions for events 1 and 2 are proportional and both outcomes are rare [scenarios (a and (b)], the oddslike 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 oddslike 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 oddslike ratio as an estimate for HR_{2} is 4.6% when HR_{2} = 1 [scenario (c)] and 6.3% when HR_{2} = 1.33 [scenario (d)]. Furthermore, the oddslike 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_{2} and 6.3% for HR_{2}; the estimates are worse with nonproportional hazards functions: in scenario (h), the relative bias is 2.6% for OR_{2} and 12.3% for HR_{2} and the CI coverage probability is 94.4% for OR_{2} and 79.4% for HR_{2}.
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 oddslike 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 β_{1} = β_{2} 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).
Association of BMI with breast cancer subtypes in the NIHAARP Diet and Health Study
We first applied the polytomous model and found that compared to leaner women (BMI < 25.0 kg/m^{2}), those who were overweight (25.0–29.9 kg/m^{2}) or obese (≥30.0 kg/m^{2}) 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 oddslike 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^{−}.
Discussion
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 timetoevent data. This approach has been underused 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 Splus 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 NIHAARP 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.
Appendix I: Proof of the Validity of a Polytomous Model to Compare Exposure–Disease Associations
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_{1} and T_{2}, are and , 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 is the baseline hazard function for the kth event type and β_{k} denotes the log HR associated with a unit increase in X for the development of the kth event type (k = 1, 2). Let D denote a polytomous outcome where D = 0 if the subject remains free of disease during followup; 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, . Then, we have and where and are the cumulative hazards functions by time t for event 1 and 2, respectively, and
Thus, and .
Therefore, and , however, . That is, a test of is equivalent to a test of .
When both event rates are rare, and are small, so that and . Thus only under the rare event assumption we can use γ_{1} to estimate β_{1} and γ_{2} to estimate β_{2}.
In the presence of independent censoring (e.g., censoring due to early loss to followup), 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 followup as long as both events have hazards functions, which are proportional to each other.
Appendix II: A Method to Simulate an Independent Exponential Censoring Time
Without censoring due to early loss of followup, the probability of no events by the end of followup τ is , where and are the cumulative hazards functions up to time τ for events 1 and 2, respectively. With the presence of censoring due to early loss of followup, the same probability becomes
where is the cumulative hazards function for an independent censoring time. We limit the increase in probability of no event with early loss of followup to be not more than 5% then the same probability without early loss of followup, 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 γ_{1} and γ_{2} 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.
Appendix III: Data Structure and Program Code
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)
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
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
Acknowledgments
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.
 Received September 13, 2012.
 Revision received December 7, 2012.
 Accepted December 10, 2012.
 ©2013 American Association for Cancer Research.