The Subtype-Free Average Causal Effect for Heterogeneous Disease Etiology

Studies have shown that the effect an exposure may have on a disease can vary for different subtypes of the same disease. However, existing approaches to estimate and compare these effects largely overlook causality. In this paper, we study the effect smoking may have on having colorectal cancer subtypes defined by a trait known as microsatellite instability (MSI). We use principal stratification to propose an alternative causal estimand, the Subtype-Free Average Causal Effect (SF-ACE). The SF-ACE is the causal effect of the exposure among those who would be free from other disease subtypes under any exposure level. We study non-parametric identification of the SF-ACE, and discuss different monotonicity assumptions, which are more nuanced than in the standard setting. As is often the case with principal stratum effects, the assumptions underlying the identification of the SF-ACE from the data are untestable and can be too strong. Therefore, we also develop sensitivity analysis methods that relax these assumptions. We present three different estimators, including a doubly-robust estimator, for the SF-ACE. We implement our methodology for data from two large cohorts to study the heterogeneity in the causal effect of smoking on colorectal cancer with respect to MSI subtypes.


Introduction
In recent decades, it has become clear that many diseases that share general clinical characteristics evolve through a range of heterogeneous molecular pathologic processes, which may be affected differently by the same exposure. For this reason, classifying a disease into subtypes according to criteria based on molecular characteristics can improve our etiologic understanding of the disease (Ogino et al., 2016).
In this paper, we consider the causal effect smoking may have on having colorectal cancer (CRC). One well-established CRC subtype classification is based on microsatellite instability (MSI). The two mutually-exclusive subtypes are MSI-high (approximately 15% of the cases) and non-MSI-high (approximately 85%). Smoking has been found to be a strong risk factor for the MSI-high subtype, while weaker evidence was found that smoking is a risk factor for the non-MSI-high subtype (Carr et al., 2018;Amitay et al., 2020).
To study etiologic heterogeneity, researchers often use a multinomial regression model in which being healthy or having each of the disease subtypes form the possible values of the outcome. However, the multinomial regression parameters do not correspond to well-defined causal effects, because the multinomial regression parameters are equivalent to the parameters obtained from a series of logistics regressions, each comparing one disease subtype to the healthy controls, resulting in selection bias (Nevo et al., 2021).
This form of bias is not limited to subtype comparisons and is not unique to multinomial regression. More generally, competing events create a challenge in making causal statements. Section 4.2 reviews related approaches and considers their applicability in our setup (Young et al., 2020;Stensrud et al., 2020;Nevo et al., 2021).
In this paper, we use data from two large cohorts to study well-defined and relevant estimands representing the effect of smoking on the two CRC subtypes. We propose an alternative estimand, inspired by the Survivor Average Causal Effect (SACE) (Zhang and Rubin, 2003;Rubin, 2006;Zehavi and Nevo, 2021). The SACE is a causal effect typicallyused to overcome the problem of truncation by death, which occurs when an outcome of interest is undefined for individuals who died before their outcome could be measured. The SACE is defined using principal stratification (Robins, 1986;Frangakis and Rubin, 2002) as the average causal effect of a treatment among individuals that would have survived under both treatment/exposure values. The SACE is not identifiable from the data under standard assumptions. Approaches to identify and estimate SACE include relying on additional strong assumptions, performing sensitivity analyses, and obtaining bounds.
Our proposed approach is based on the same idea as the SACE. First, we define principal strata with respect to the potential outcomes of other disease subtypes. Then, we define the Subtype-Free Average Causal Effect (SF -ACE) to be the average causal effect of the exposure on one disease subtype among the individuals who would have been free of the other MSI disease subtype under either exposure level (smoking status).
One commonly-made assumption in truncation-by-death problems is Monotonicity, namely that treatment/exposure cannot hurt survival (Zhang and Rubin, 2003). In the setting of disease heterogeneity with two subtypes, the definition of Monotonicity is more nuanced for two reasons. First, the assumption needs to be made (or not) for each subtype separately. Second, the disease could occur under both exposure statuses, but the subtype may differ. We term this phenomenon exposure-induced subtype switching and show it plays a key role in the (lack of) identification of the SF -ACE. We further show that unlike the SACE setting, the SF -ACE is non-parametrically identifiable under certain Monotonicity assumptions. Under weaker Monotonicity assumptions, the SF -ACE is no longer identifiable. For this scenario, we develop a sensitivity analysis approach as a function of the unidentifiable exposure-induced subtype switching probabilities.
We consider three SF -ACE estimators. The first two, standardization-based estimator and an inverse probability of treatment weighting (IPTW) estimator, may rely on parametric modeling assumptions. We therefore present a doubly-robust (DR) estimator that will be consistent if at least one of two models, but not necessarily both, is correctly specified.
The rest of the paper is organized as follows. Section 2 describes the epidemiologic problem and the available data. Section 3 presents the notations and the causal framework. Section 4 reviews existing approaches for etiologic disease heterogeneity and discuss available causal estimands. Section 5 presents a detailed study of the proposed causal estimand, the SF -ACE, including its interpretation, assumptions needed for its identification, sensitivity analyses and estimation. Section 6 summarizes the results of our simulation studies. Section 7 presents various analyses of the data to study the causal effect of smoking on the two CRC subtypes. Final conclusions are given in Section 8. The R package TheSFACE is available from CRAN implements our methodology. Reproducibility materials are available from https://github.com/amitSasson/SFACE Reproduce.

Data and problem description
We use data from two large US cohorts. The Nurses Health Study (NHS)  CRC diagnoses were reported by study participants. Formalin-fixed paraffin-embedded tumor tissue specimens were retrieved from hospitals across the US, and a CRC diagnosis was confirmed by the study pathologist (S.O.), who also conducted analyses to determine the MSI status. Further details can be found in Ugai et al. (2022) and references therein.
We set the baseline to be age 60, and define the exposure to be ever versus never smoking by age 60. Such binary exposure definitions are often used to minimize recall bias and measurement error. The outcome was CRC by age 70, subtyped by MSI status. Age 70 was chosen to minimize the impact of death as a competing event.
We removed from the data people who were older than 60 when enrolled in the study, or people who died, had CRC, or were lost to follow up before age 60. As in Ugai et al. (2022), we also removed people with missing baseline smoking data. In total, we were left with 114,947 people at baseline. Then, we removed 5,132 people (4.4%) who have died, without being diagnosed, prior to age 70. We were left with a final sample of 109,815. In this sample, 17,854 people were CRC-free and younger than 70 when the last questionnaire was collected (out of which 94.8% were older than 65). We treated them as CRC-free at age 70, and examined how it may affected our results by considering alternative analyses.
In our sample, the number of exposed (ever smoking) was 58,432 (53.3%) and the number of CRC cases by age 70 was 961 (0.9%). Out of them, 358 were diagnosed with non-MSI-high CRC (37.3%) and 61 were diagnosed with MSI-high CRC (6.3%). The subtype was missing for the remaining 542 (56.4%) diagnosed participants, as is often the case in such studies (Nevo et al., 2018;Liu et al., 2018). To have a clear presentation of the issues at end, we develop our theory and methodology under the assumption that subtypes are known for all disease cases. We address the issue of missing subtypes in our analysis in Section 7 using inverse probability weighting (Liu et al., 2018).
The NHS and HPFS cohorts include a vast number of covariates collected through biennial questionnaires. Table C.7 in the Appendix provides information on key baseline confounders we adjusted for in our analyses: gender (male/female), regular aspirin use (yes/no), first-degree family history of CRC (yes/no), history of lower endoscopy (yes/no), physical activity level (mean metabolic equivalent task score hours per week), body mass index (categorical), alcohol intake (g/day) and total calorie intake (kcal/day).

Preliminaries
Using the potential outcomes framework, let Y i (a) denote the potential disease status had the exposure of individual i = 1, ..., n been set to A i = a for a = 0, 1. In our study, A i = 0 indicates no smoking by age 60 (as captured by the cohort data), while A i = 1 indicates ever smoking by age 60. If individual i would have been diagnosed with disease subtype k = 1, ..., K under exposure value A i = a, then Y i (a) = k, and if they would be disease-free Note that because the subtypes are i (a) ≤ 1. In our study, K = 2, subtype 1 is non-MSI-high CRC, and subtype 2 is MSI-high CRC. So, for example, Y i (0) = 1 if person i would have been diagnosed with non-MSI-high CRC by age 70 had they not been smoking by age 60.
Throughout the paper, we assume there is no interference between individuals and that there are no multiple versions of the exposure leading to different potential outcomes.
Under these assumptions, known collectively as the stable unit treatment value assumption (SUTVA), the observed disease status for person i, hereafter denoted by Y i , with observed exposure value A i , equals to the appropriate potential outcome for that person. That is, . The latter is also known as the consistency assumption.
As in many observational studies, we also assume weak ignorability, namely that rich enough data were collected as part of the NHS and HPFS cohorts such that for a vector of measured covariates X, we have A |= Y (a) | X for a = 0, 1. Another typically needed assumption for identification of causal effects is positivity, meaning that the exposure assignment probability e(x) = Pr(A = a|X = x) is positive for a = 0, 1 and for all x.
In summary, the observed data consists of n i.i.d copies sampled from the distribution of (X, A, Y (1) , ..., Y (K) ), that can be equivalently represented as (X, A, Y ).

Review of existing approaches
A principled approach to causal inference typically starts with a definition of causal effects of interest, followed by the assumptions needed for their identification from observed data, and only then statistical models are discussed. However, as the existing approaches employed in practice for studying etiologic heterogeneity are model-based and are not grounded in the causal inference paradigm, we start with reviewing these models (Section 4.1). Then, we discuss relevant causal effects arising in the presence of competing events (Section 4.2).

Statistical models for disease heterogeneity
Let π k (a, x) = Pr(Y = k|A = a, X = x), k = 0, ..., K, be the probability of subtype k disease (or being subtype-free for k = 0) among those with exposure A = a and confounders X = x. Let also π(a, x) = (π 0 (a, x), ..., π K (a, x)). A commonly-used approach for studying etiologic heterogeneity is built upon the multinomial regression model for k = 1, ..., K and π 0 (a, x) = 1 − K k=1 π k (a, x). Under this model, exp(β k ) is the odds ratio (OR) between the exposure A and the observed subtype-specific outcome Y (k) compared with not having the disease, that is, exp(β k ) = π k (1,x) π 0 (1,x) × π 0 (0,x) π k (0,x) . The parameter exp(β k ) also approximates the risk ratio (RR) when the disease is rare, as in our study.
Model (1) was also extended for studies within which subtypes were not a-priori known and may be defined using multiple, possibly mis-meausred, biomarkers (Chatterjee, 2004;Wang et al., 2015;Nevo et al., 2016). See also Wang et al. (2016); Zabor and Begg (2017) for a survey of and comparison between methods based on this and similar models.
Estimation of Model (1) parameters is typically carried out by standard maximum likelihood approach, with some exceptions when further complications arise (e.g., Chatterjee, 2004;Nevo et al., 2016). In practice, researchers are often interested in whether there is heterogeneity in the effects of A across disease subtypes. When the null hypothesis H 0 : β 1 = ... = β K is rejected, researchers conclude that the exposure A contributes to the risk of at least one subtype differently from its contribution to the risk of other subtypes.
A recently-proposed alternative approach postulates K logistic regressions, each comparing a subtype k disease to the alternative of being free of that subtype, that is, being either disease-free or diagnosed with the disease, but of a different subtype (Sun et al., 2017). Formally, these models can be expressed by for k = 1, ..., K. Under this model, exp( β k ) is the OR between the exposure A and the observed indicator of subtype-specific outcome Because fitting these models separately may result in K k=0 π k = 1, estimation of the parameters in (2) calls for constrained maximization of the models' likelihood. A more computationally attractive approach, developed by Sun et al. (2017), uses constrained Bayesian estimation, and limits the posterior distribution of the parameter vector to have positive mass only for parameter values such that K k=0 π k = 1.
Comparing models (1) and (2), it has been recently argued (Sun et al., 2017) that variation in β k across the subtypes (across k) is not a satisfactory measure of etiologic heterogeneity, because the parameter space of a multinomial regression model requires a risk factor for one subtype to be a risk factor for the other subtypes, and that β k may be preferred, although this issue was also debated (Begg et al., 2018;Sun et al., 2018).

Causal effects for competing risks
To fix ideas, we focus on the effect of A on Y (1) , the first subtype, and start with a naïve approach that does not correspond to a causal effect. Researchers may opt for a comparison of the subtype 1 proportion of cases at each exposure level only among individuals who are free of other disease subtypes. When researchers also adjust for the confounders X, their statistical estimand on the difference scale is and can be analogously defined on the RR or OR scales. Because model (1) is equivalent to K separate logistic regression models, each comparing one disease subtype to the disease-free individuals (Wang et al., 2016), then exp(β k ) is the OR analogue of (3), and approximates the RR as CRC is a rare disease. However, (3) does not correspond to a causal effect, because of selection bias due to common causes of the different subtypes (Nevo et al., 2021). These common causes cannot be classified as classical confounders, as they can be independent of the exposure. If these common causes are also mediators of an exposure-subtype mechanisms that are shared between more than one subtype, even had they been observed, adjusting for them would have eliminated part of the causal effect of the exposure on the subtype (Nevo et al., 2021). Importantly, because MSI-high and non-MSI-high are subtypes of the same disease, CRC, they are likely to have unknown and/or unmeasured common causes. Young et al. (2020) considered possible causal estimands in the presence of competing events. The total effect, T E (1) = E[Y 1 (1) − Y 1 (0)], contrasts the population-level subtype 1 risk had the entire population been exposed, with the population-level subtype 1 risk had the entire population were unexposed. In the presence of measured confounders, this effect is identified from the observed data under SUTVA, weak ignorability and positivity by Unlike (3) is sensitive to the phenomenon of exposure-induced subtype switching. For example, assume some individuals would have been diagnosed with subtype 1 under A = 0 and would have been diagnosed with subtype 2 under A = 1, while the rest of the population will be free of subtype 1 for all A values. Then, the total effect on subtype 1 will be negative, T E (1) < 0, and it would appear as if the exposure protects from subtype 1, even though the individuals "protected" by A would still have the disease, just of another subtype.
The direct effect (Young et al., 2020) is the causal effect comparing subtype 1 rates under two hypothetical joint interventions. The first sets the exposure to A = 1 and eliminates the possibility of disease subtypes other than subtype 1 disease, and the second sets the exposure to A = 0 and also eliminates risks of other disease subtypes. This approach is, however, irrelevant for our study, because an intervention eliminating the risk of a particular CRC subtype is not currently available and is not expected to be in the near future.
Another alternative considers the separable effects (Stensrud et al., 2020), which rely on a conceptual separation of the exposure into the components affecting each disease subtype.
Such separation requires profound understanding of the mechanisms leading to each CRC subtype. Furthermore, the conditions needed for identifiability of these effects are violated if there are unobserved common causes of both disease subtypes (Stensrud et al., 2020).
Such common causes are expected to exist in our study and in many disease subtype studies (Nevo et al., 2021).

The Subtype-Free Average Causal Effect
Our focus in this paper is on the common scenario of two disease subtypes (K = 2), as in our study. We propose a new estimand, inspired by the SACE (Rubin, 2006). We define the SF -ACE on the difference scale by

SF -ACE
(1) to be the SF -ACE on the RR scale. The SF -ACE (1) is the average causal effect among the individuals that would have been free of disease subtype 2 under either exposure level (i.e., regardless of smoking status). Similarly to SF -ACE (1) , let

SF -ACE
(2) be the average causal effect on the difference scale of A on Y (2) among the individuals that would have been free of disease subtype 1 under either exposure level. Let also SF -ACE

Interpretation of the SF -ACE
As explained in Section 4.2, competing events present a challenge in defining causal estimands that shed light on the scientific question of interest. The SF -ACE presents another source of information. A positive SF -ACE (1) D value means that in the sub-population who would have been free of MSI-high CRC regardless of their smoking status, we will expect more non-MSI-high CRC cases had the entire sub-population smoked compared to the scenario the entire sub-population did not smoke before age 60.
Therefore, the SF -ACE offers a way to circumvent non-zero causal effects due to exposure-induced subtype switching. Nevertheless, the above appealing advantage does not come without a price. One disadvantage is that the SF -ACE is defined in a latent subset of the population. This is a general problem with principal stratification approaches.
A second, related, issue is that when contrasting SF -ACE (1) and SF -ACE (2) , we contrast effects defined in two overlapping, but not identical, sub-populations. The SF -ACE (1) is defined within the sub-population {Y (2) (0) = Y (2) (1) = 0}, and the SF -ACE (2) is defined Thus, a non-zero θ D captures changes that might be due to the different population, and as a result, the exact source of heterogeneity in the causal effect is not revealed. Instead, a non-zero θ D presents evidence of heterogeneity and calls for further investigation. Lastly but importantly, in rare disease scenarios, such as CRC, this problem is less concerning because the two sub-populations {Y (2) (0) = Y (2) (1) = 0} and {Y (1) (0) = Y (1) (1) = 0} have almost perfect overlap, as those who would not have CRC of any subtype under any exposure value belong to both groups and comprise the vast majority of the population.

Identifiability
Similar to the SACE in truncation-by-death scenarios, the SF -ACE is not identifiable from the data under the standard assumptions of SUTVA, weak ignorability and positivity. One commonly-made assumption when targeting the SACE is Monotonicity (Zhang and Rubin, 2003), stating that treatment cannot hurt survival. In our study, the definition and plausibility of Monotonicity-like assumptions are more nuanced, and, as it turns out, resulting in new insights regarding principal causal effects and their identification in a competing risks setup. We start with the definition of Subtype Monotonicity (S-Monotonicity). It is clear from the above that the SF -ACE RR is identifiable from the data under SUTVA, weak ignorability, positivity, and S-Monotonicity for both subtypes. The following proposition establishes identification of SF -ACE D , and presents identification formulas under either standardization or IPTW (Hernán and Robins, 2020).

Assumption 1 Subtype Monotonicity (S-Monotonicity
Proposition 1 Under SUTVA, weak ignorability, positivity, and S-Monotonicity for sub- using standardization, or via IPTW by The proofs and analogous expressions for SF -ACE (2) , θ D , and for causal contrasts on the RR scale are all given in Section A of the Appendix. Estimation based on these identification formulas is discussed in Section 5.5.

Relaxing S-Monotonicity
Assuming S-Monotonicity for both subtypes might be too restrictive, even if we believe the exposure cannot protect from the outcome, as S-Monotonicity disallows subtype switching for any individual in the population. For example, consider the scenario the exposure is believed to increase non-MSI-high CRC risk but not MSI-high CRC risk. If individual i would have been diagnosed with MSI-high CRC when unexposed, then under S-Monotonicity for MSI-high CRC, it will be impossible for individual i to be diagnosed with non-MSI-high CRC when exposed. This is even though the exposure is a risk factor only for non-MSI-high CRC.
Therefore, we consider a weaker assumption, termed Disease Monotonicity (D-Monotonicity), that states that an individual who would have been diagnosed with disease subtype k under no exposure, would have also been diagnosed with the disease under exposure, but not necessarily with the same subtype.
D-Monotonicity keeps the premise that the exposure cannot protect from the disease while being less restrictive than S-Monotonicity. As with S-Monotonicity, D-Monotonicity can be modified for scenarios where the exposure is suspected to have a protective effect.
For each disease subtype, we may assume S-Monotonicity, D-Monotonicity, or neither.
With two disease subtypes, we have eight different combinations of possible assumptions.

Under each combination, certain potential outcome profiles of {Y
(1) i (1)} are assumed to be absent from the population. For example, under S-Monotonicity for subtype 1, there is no one in the population with {Y (1) i (1)} = {1, 0, 0, 1}. Note also that because the subtypes are mutually exclusive, it is impossible to i (1) = 1 regardless of the taken monotonicity assumptions. Table 1 presents the nine possible potential outcome profiles, and specifies whether each profile exists under four different combinations of assumptions. This table demonstrates that the exact implications of the different monotonicity assumptions on the potential outcomes in the population are quite explicit. The observed data, however, does not hold all the information required for classifying each individual into their potential outcome profile. Table A.1 in Section A of the Appendix specifies the possible profiles for each participant according to their observed (A, Y (1) , Y (2) ) under each of the four combinations of assumptions.

Potential outcomes profiles
Monotonicity assumptions

A sensitivity analysis approach
By Proposition 1, under S-Monotonicity for both disease subtypes, the SF -ACE (1) , SF -ACE (2) and θ are identifiable from the observed data distribution. If we replace S-Monotonicity with D-Monotonicity or with no monotonicity assumption for at least one subtype, then the causal effects are no longer identifiable, and SF -ACE Therefore, we take a sensitivity analysis approach, as is common when studying causal effects under untestable assumptions. Let λ 1 = Pr[Y (2) (1) = 1|Y (1) (0) = 1] be the proportion of individuals who would have been diagnosed with disease subtype 2 when exposed, out of the individuals who would have been diagnosed with disease subtype 1 when unexposed. That is, λ 1 is the probability of exposure-induced subtype switching from subtype 1 to subtype 2. Note that λ 1 = 0 under S-Monotonicity for subtype 1 and λ 1 ≥ 0 under D-Monotonicity for subtype 1. Let also λ 2 = Pr[Y (1) (1) = 1|Y (2) (0) = 1] be the subtype-switching probability from subtype 2 to subtype 1. While λ 1 and λ 2 are not identifiable from the observed data, they can be bounded. Specifically, because The following two propositions demonstrate that the subtype-switching probabilities play a key role in the identification of the SF -ACE whenever the S-Monotonicity assumptions are not met. Proposition 2 shows that (a) under D-Monotonicity for subtype 1 and S-Monotonicity for subtype 2, both the SF -ACE

Proposition 2 Under SUTVA, weak ignorability, positivity and D-Monotonicity for both subtypes, the SF -ACE
(1) D is identifiable from the data by using standardization, or via IPTW by .
The identification of SF -ACE (2) D is analogous.
The proof is given in Section A of the Appendix. Of note is that if one further assumes S-Monotonicity for one of the subtypes (say subtype 2), then an identification formula of

SF -ACE
(1) D is obtained by setting the appropriate λ k (say λ 2 ) to zero in (8) and (9). A second remark is that, technically, results (8) and (9) require D-Monotonicity for subtype 2 only, and similarly, the identification for subtype 2 requires D-Monotonicity for subtype 1 only. In Section A of the Appendix, we present further identification results for SF -ACE D under D-Monotonicity for one subtype only or without any monotonicity assumptions, and show that in these cases λ 1 and λ 2 do not suffice for identification of SF -ACE D .
Turning to the RR scale, Proposition 3 asserts that SF -ACE RR is identifiable as a function of λ 1 and λ 2 and the observed data without imposing any monotonicity assumptions.

The identification of SF -ACE
(2) RR is analogous.
The proof is given in Section A of the Appendix.
Propositions 2 and 3 serve as basis for our proposed sensitivity analyses. Researchers first decide what is the effect scale of interest, and then for which if any of the two subtypes the S-Monotonicity or D-Monotonicity assumptions are plausible. Then, depending on the taken assumptions, they may set λ 1 or λ 2 to zero, and estimate the SF -ACE (k) and θ as a function of λ 1 and/or λ 2 . Such an analysis can reveal how sensitive the results are to different levels of exposure-induced subtype switching. We implement these analyses in our study of smoking effect on CRC MSI subtypes in Section 7.

Estimation
Recall that for each individual i, the observed data are ( i ).
Estimation of SF -ACE (k) and θ from the data alongside the sensitivity analyses for the S-Monotonicity assumptions require estimation of π k (a, x), k = 1, 2, a = 0, 1 for the standardization-based formulas, or e(x) for the IPTW-based identification formulas. Expectations over X or (X, A, Y (k) ) are estimated by taking the sample analogues.
Non-parametric estimation of π k (a, x) or e(x) when the vector of measured confounders X is of high dimension or includes many continuous variables could be a difficult task, that will yield an estimator with a too large variance to be able to derive conclusions from the data. Therefore, parametric estimation might be a better option, and also aligns with the models used in practice and that were presented in Section 4.1. We first discuss parametric estimation, before implementing a DR estimator that will be consistent even if not all models are correctly-specified.
Starting with standardization, let π(a, x; η 1 ) = (π 0 (a, x; η 1 ), π 1 (a, x; η 1 ), π 2 (a, x; η 1 )) be the parametric model for π(a, X) with parameter vector η 1 . This could be, for example, the regression models (1) or (2). Upon estimatingη 1 , the estimated probabilities π(a, X i ;η 1 ) are calculated in the entire sample, for a = 0, 1, before plugged in the relevant identification formulas. For example, the proposed SF -ACE If there is an unmeasured covariate U affecting both disease subtypes (Y (1) , Y (2) ), model (1) might be misspecified, if one believes the model that is also conditioned on U is a multinomial regression model. However, this is not a main concern in studies involving rare disease, as in our study (Nevo et al., 2021).
Turning to IPTW, let e(x; η 2 ) be a parametric model for e(x), with parameter vector η 2 , for example a logistic regression. Then, after estimatingη 2 , the SF -ACE (1) D estimator is obtained by and similar estimators for SF -ACE (2) D and θ D , or for effects on other scales, can be obtained.
When using the maximum likelihood approach to fit models π(a, x; η 1 ) or e(x, η 2 ), the obtained estimators η 1 or η 2 are consistent and asymptotically normal. Thus, the continuous mapping theorem and the delta method can be used to show that the SF -ACE and θ estimators are also consistent and asymptotically normal. In practice, the bootstrap can be used to estimate the standard errors (SEs) and construct confidence intervals (CIs).
The standardization estimator or the IPTW estimator will be consistent only if their respective models are correctly specified. Therefore, we extend our estimation approach to DR estimation (Robins et al., 1994;Lunceford and Davidian, 2004). That is, we construct an estimator that will be consistent whenever at least one of the two models, π(a, x; η 1 ) or e(x; η 2 ), but not necessarily both, is correctly specified. DR estimators are frequently used in causal inference to reduce dependency on the chosen model. Furthermore, a DR estimator constructed as an augmented-IPTW estimator is also more efficient than IPTW estimators and have desirable theoretical properties (Robins et al., 1994).
Letη 1 andη 2 be the estimators obtained for η 1 and η 2 . As can be seen from the proof of Proposition 1 (Section A of the Appendix), under S-Monotonicity, the SF -ACEs and . Under weak ignorability and positivity, a DR augmented-IPTW (Bang and Robins, 2005;Lunceford and Davidian, 2004). The estimator (12) and analogous and similar equations to obtain estimators for the causal effects that are consistent by the continuous mapping theorem, if at least one of the models, π(a, x; η 1 ) or e(x; η 2 ), is correctly specified. Similar arguments can be used to construct DR estimators when conducting sensitivity analyses based on Proposition 2 and 3.

Simulations
We conducted simulation studies to evaluate the finite-sample performance of our proposed estimators and to compare the SF -ACE to existing estimands. Detailed information about the data generating mechanism and parameters used, as well as more detailed results are given in Section B of the Appendix. For each i = 1, ...n, potential outcomes Y i (a) were simulated using a multinomial regression model with two confounders, X 1i ∼ Ber(0.5), X 2i ∼ N (0, 1), and one unobserved covariate U i ∼ N (0, 1) affecting both subtypes. S-Monotonicity was imposed for both subtypes (without violating the multinomial model assumption, see Section B). The exposure A i was simulated using a logistic regression with To study different aspects of the proposed approach, we conducted three studies. In the first, we fixed the model parameters and varied the sample size between 5,000 and 50,000, to study the methods' performance in different sample sizes. In the second study, we fixed the sample size and varied the effect the unmeasured covariate U i had on Y (2) i to compare between the different causal effects. In the third study, we investigated the methods' performance when the fitted models were misspecified.
In each simulated dataset, we first fitted the multinomial regression (1), and a logistic regression for e(x). Then, we calculated the conditional estimators (3), the TE estimators (4), and the three SF -ACE estimators described in Section 5.5. We used the bootstrap with 200 repetitions to estimate the SEs and to construct Wald-type 95% CIs.
In the first study (Tables B.3 and B.4), the SF -ACE estimators had a small relative bias. The SE of all estimators were similar, and were generally well-estimated. Empirical coverage rates of the CIs were satisfactory. These results were also apparent for the second study ( Figures B.1 and B.2). In the second study, as the effect of U on Y (2) increased, absolute bias of the conditional estimand and the TE on the difference scale, compared to the true SF -ACEs, increased, while the SF -ACE estimators remained unbiased. As expected, because data were simulated under S-Monotonicity for both subtypes, the T E

Data analysis
We conducted several analyses to study the causal effects of ever smoking by age 60 on having CRC MSI subtypes by age 70, using the NHS and HPFS datasets described in Section 2. The main results are presented here, and additional results are given in Section C of the Appendix. In all of our analyses, we used the confounders detailed in Section 2 and Table C.7. Wald-type 95% CIs and p-values were calculated using the bootstrap with 200 repetitions.
As indicated in Section 2, for more than 50% of CRC cases the subtype status was unknown. To minimize the impact of potential selection bias, we used inverse probability weighting for missing subtype (Liu et al., 2018). We first fitted a model for the probability of observing the subtypes (Table C.8). The model included additional variables available for CRC cases (e.g., tumor stage and size). Then, for each CRC case with known subtype, we calculated the weight as the reciprocal of the probability of observing the subtype. The weights were truncated at the 99% percentile and then incorporated into our estimators.
We estimated the TE of smoking on each disease subtype by standardization (4), using a multinomial regression model (1) (Table C.9). The point estimates indicated smoking increases the risk of both CRC MSI subtypes (Tables 2 and C.10), although the results were not significant at the 5% level. The point estimates indicated that the effect on non-MSI-high subtype was larger than the effect on the MSI-high subtype on the difference scale and vice versa on the RR scale. These results are in agreement with the existing literature, because MSI-high is more rare, so while the RR is larger, the number of excess cases due to smoking can be smaller.   The S-Monotonicity assumption might be too restrictive, especially for the non-MSIhigh subtype, for which there is limited information about the relative effect strength (Amitay et al., 2020). From a subject-manner perspective, we could not overlook the possibility of exposure-induced subtype switching for non-MSI-high subtype, as the relative effect of smoking on MSI-high CRC is known to be larger. Therefore, we replaced the S-Monotonicity assumption for non-MSI-high (subtype 1) with D-Monotonicity by fixing λ 2 = 0 and varying the values of λ 1 between 0 and (Figure 1). As λ 1 increased, the estimated SF -ACE for non-MSI-high CRC increased, the estimated SF -ACE for MSI-high CRC decreased. For λ 1 ≥ 0.04, H 0 : θ D = 0 was rejected. This result serves as evidence for heterogeneity in the SF -ACE as long as we believe subtype switching is possible and not extremely rare. Given some a-priori evidence for smoking effects on both subtypes, we also considered an analysis under D-Monotonicity for both subtypes. Figure 2 presents the estimated effects on the difference scale as a function of λ 1 and λ 2 (for which the data did not restrict the range), along with 95% CIs. For nearly the entire range, the effect on non-MSI-high CRC was significant and the effect on MSI-high CRC was not. Sensitivity analysis results were similar for all three estimators (Figures 2, C.4 and C.5). As λ 2 increased, the larger was the minimal value of λ 1 for which heterogeneity can be concluded. For the RR (Figures C.6-C.8), under the combination of low λ 2 and large λ 1 the heterogeneity assumption was rejected.
As described in Section 2, there were 17,854 study participants for which the last available data were before age 70 (and they were not diagnosed with CRC and did not die).
We treated these people as CRC-free at age 70. To assess the impact of this decision, we repeated our analyses while removing these people. In an additional analysis, we again removed these people but used inverse probability of weighting to take them into account.
The results did not change materially, although the estimated effect sizes and the standard errors were generally larger (Tables C.13 and C.14).
To summarize our main findings, we have provided evidence for smoking effect on both CRC MSI subtypes and evidence for heterogeneity under reasonable causal assumptions.
Even though our case numbers were smaller and our follow-up was more limited than previous studies, the estimated SF -ACE lend support of a causal effect on non-MSI-high CRC, for which available evidence was limited.

Discussion
Questions of etiologic heterogeneity are central to our understanding of how diseases evolve.
Viewed through the lenses of competing events, existing approaches largely overlooked formal discussions of causal estimands and identifiability assumptions. On the other hand, the recent progress in causal inference methodology for competing events presented causal estimands that either provide little information about etiology heterogeneity or estimands requiring identifiability assumptions that are implausible in the study of CRC subtypes.
In this paper, to study the effect of smoking on CRC subtypes defined by MSI status, we proposed the SF-ACEs as causal effects that can complement existing methods and provide further insight towards the study of etiologic heterogeneity. One challenge with using multinomial regression with cohort data is the issue of right censoring due to loss to follow-up or administrative censoring. To assess its impact, we considered in this paper a number of analyses. The censoring problem would have been more pronounced if follow-up was longer than 10 years. Therefore, considering time-to-event outcomes is an important topic of future research. Another potential extension will consider continuous exposures instead of binary ones. Specifically, the generalization of the monotonicity assumptions, both in terms of definition and plausibility, is of interest. • Section A includes proofs and additional theory. Section A.1 presents the proofs of Propositions 1-3. Section A.2 illustrates the the connection between observed data and the potential outcomes. Section A.3 discuss identification results under weaker combination of assumptions. Section A.4 provides the complete expression for the DR estimator.
• Section B includes details on the simulations and presents the results. Section B.1 gives additional details about the data generating mechanism of the simulation. Section B.2 presents the details of three different simulation studies. Finally, Section B.3 presents tables and figures summarizing the different simulation results.
• Section C presents details on the data as well as additional analyses and results. Section C.1 includes descriptive information about the data used. Section C.2 gives additional information about the weights used for handling missing subtypes and about the models fitted for estimating effects. Finally, Section C.3 presents additional data analysis results.

A.1 Proof of identification results
Before providing the identification formulas for the SF -ACE, we first show that the SF -ACE equals to the T E when defined on the RR scale under S-Monotonicity for both disease subtypes.
When defined on the RR scale, the SF -ACE (1) equals to (A.14) Equality (A.13) holds by the S-Monotonicity assumption for subtype 2. Equality (A.14) holds because the subtypes are mutually exclusive (numerator) and by the S-Monotonicity assumption for subtype 1 (denominator).

A.1.2 Proof of Proposition 1
Starting from the definition of SF -ACE, we have

SF -ACE
(1) For the standardization-based identification formula, we continue by recalling that under positivity and weak ignorability, Pr[Y (k) (a) = j] is identified by for a, j, k = 0, 1 (Hernán and Robins, 2020). With Pr(Y (k) = j|A = a, x) previously denoted by π k (a, x) for j = 1 and hence equals to 1 − π k (a, x) for j = 0, we substitute these expressions in (A.16) to get

(A.18)
The identification formulas for subtype 2 are analogous. Using standardization,
The difference between the effects is given by substituting the expressions above in θ D =
Finally, define θ RR to be the difference between the effects in the RR scale. So θ RR is given by

A.1.3 Proof of Proposition 2
We start with the effect on subtype 1.

SF -ACE
In (A.19), the numerator follows from the fact the subtypes are mutually exclusive and the denominator is justified by the law of total probability. In (A.20), the numerator holds because of the law of total probability and the denominator follows by noting that the D-Monotonicity assumption for subtype 2 implies Pr[Y (2) (0) = 1, Now, similarly to the proof of Proposition 1, we can continue using standardization, by weak ignorability and positivity. Alternatively, as in the proof of Proposition 1, we can continue from (A.21) using IPTW to obtain

SF -ACE
(1) , again by weak ignorability and positivity, similarly to the standard IPTW-based identification.
The identification formulas for subtype 2 are analogous. Using standardization,
The difference between the effects is given by

SF -ACE
(1) Equality (A.22) holds because the subtypes are mutually exclusive. Equality (A.23) is by the law of total probability. As before, using standardization, we may write
The identification formula of subtype 2 is analogous. Using standardization,

SF -ACE
(2) The difference between the effects is given by θ RR = SF -ACE (1) RR .

An interesting distinction between the results in Propositions 2 and 3 compared to
Proposition 1, is that if one further assumes S-Monotonicity for one of the subtypes only, such that S-Monotonicity holds for one subtype only, say subtype 1 so λ 1 = 0, then the identification formulas for the two subtypes are no longer symmetric. Table 1 of the main text defines the nine possible potential outcome profiles. Table A Observed data Monotonicity assumptions (1, 0, 1) 1, 6 1,6,7 1,6 1,6,7 Table A.3: Connection between observed data and possible potential outcome profiles defined in Table 1 under four combinations of monotonicity assumptions.

A.3.1 S-Monotonicity for one subtype and no monotonicity assumption for the other subtype
We present here identification for the SF -ACE (k) D under S-Monotonicity for subtype 1 and no monotonicity assumption, of any version, for subtype 2. To this end, in addition to λ 2 , we define an additional sensitivity parameter λ 0 This sensitivity parameter is the probability to be CRC-free when exposed among those who would have been diagnosed with subtype 2 had they were unexposed. Under either S-Monotonicity or D-Monotonicity, λ 0 2 = 0.
We begin with the effect on subtype 1,

SF -ACE
(1) Because different assumptions were made for subtype 1 and subtype 2, the effects are not analogous. Turning to the effect on subtype 2,

SF -ACE
(2) Equality (A.27) holds because of the S-Monotonicity assumption for subtype 1. Equality (A.28) holds because the subtypes are mutually exclusive. Equality (A.29) holds because of the law of total probability. As before, standardization, IPTW or DR estimation can be used to identify and estimate the effects.

SF -ACE
(1) Equality (A.30) holds because the subtypes are mutually exclusive. Equality (A.31) holds because of the law of total probability. Equality (A.32) holds because of the law of total probability and because the subtypes are mutually exclusive.

SF -ACE
(2) As before, standardization, IPTW or DR estimation can be used to identify and estimate the effects.

A.4 Details about the doubly-robust estimator
The doubly-robust estimators that are plugged in Equations (A.16) (under S-Monotonicity for both subtypes), (A.21) (under D-Monotonicity for both subtypes), or their RR ana-logues are given by,

B Additional simulation details and results
Section B.1 gives additional details about the data generating mechanism of the simulation. Section B.2 presents the details of three different simulation studies. Finally, Section B.3 presents tables and figures summarizing the different simulation results.

B.1 Data generating mechanism
The data were simulated to resemble an observational study. For each individual i, we first simulated two independent measured confounders, X i = (X i1 , X i2 ), where X 1 ∼ N (0, 1) and X 2 ∼ Ber(0.5) and one normally distributed unmeasured covariate U i ∼ N (0, 1). To simulate the potential outcomes, we used the multinomial regression model for k = 1, 2, and We first simulated the potential outcome under no exposure Y i (0), for each individual i, using model (A.33). To make sure the data followed the S-Monotonicity assumption for both disease subtypes, the potential outcomes under exposure (a = 1) were created as follows: Individuals with Y i (0) > 0, namely that they would have been diagnosed with disease subtype 1 or 2 under a = 0, would automatically be diagnosed with the same disease subtype under a = 1. For individuals who were disease-free under a = 0 the disease subtype under a = 1 was simulated using Equation (A.33) with an adjusted probability using the Law of Total Probability, We simulated the actual exposure status for each individual A i using a logistic regression model, .
The observed outcome Y was obtained for each individual by setting For each individual i, the observed data were (A i , X i , Y i ).

B.2 Additional details about the simulation studies
As indicated in the main text, in each of the studies and scenarios described below, we first estimated π(a, x; η 1 ) using a multinomial regression model and e(x) using a logistic regression model. Then, we computed the estimators for the conditional statistical estimand, the total effect, and the three SF -ACE estimators.
Three studies were investigated. In Study I, we fixed the model parameters and varied the sample size between 5,000 and 50,000. In Study II, we varied δ 2 , the effect of the unmeasured covariate U on Y (2) , while keeping all the other parameters fixed. In Study III, we investigated the estimators' performance when the models fitted were misspecified. To achieve this, a transformation of X 2 was used to create the data, but the analysis assumed X 2 is (logit) linearly connected to the exposure or/and the outcome. Misspecification of the A ∼ X relationship was created by taking log |X 2 | when generating the data instead of X 2 . Misspecification of the Y ∼ A + X model was also created by replacing X 2 with log |X 2 | when generating the data. Table B.4 gives the model parameters and sample sizes used in the three studies. It also gives the obtained true SF -ACE values for each set of parameters yields. In all studies, the prevalence of subtype 1 ranged between 5%-6% and the prevalence of subtype 2 ranged between 2%-4%. The prevalence of the exposure was 49%.

B.3 Simulation Results
Tables B.5 and B.6 present the result for Study I. As can be seen from the Tables, the SF -ACE estimators had a smaller relative bias than the other estimators, even with a small sample size.
The absolute difference between non-SF -ACE estimators and the SF -ACEs on the difference T E and SF -ACE coincide on the RR scale. The relative bias ( Figure B.4) remained generally the same. As all three methods (standardization, IPTW, and DR estimation) were based on correctly-specified models, their finite sample bias was minimal.  Finally, Table B.7 presents results from Study III under four settings: (a) Both π(a, x; η 1 ) and e(x; η 2 ) are correctly specified, (b) The π(a, x; η 1 ) model is correctly-specified and model e(x; η 2 ) is misspecified, (c) The π(a, x; η 1 ) model is misspecified and model e(x; η 2 ) is correctlyspecified, and (d) both models are misspecified. As expected the standardization-based estimator had minimal bias when π(a, x; η 1 ) was correctly specified, and was biased when the model was misspecified. The IPTW-based model had minimal bias whenever model e(x; η 2 ) was correctly specified and produced considerable bias when the model was misspecified. As expected the DR estimator was only biased when both models were misspecified and had minimal bias under scenarios (a), (b), and (c) described above.

C Details on the data analysis and additional results
Section C.1 presents details on the data as well as additional analyses and results. Section C.1 includes descriptive information about the data used. Section C.2 gives additional information about the weights used for handling missing subtypes and about the models fitted for estimating effects. Finally, Section C.3 presents additional data analysis results. Table C.9 gives basic descriptive information about the different covariates X used in our analyses.  (1976-1995, 1996-2000, 2001-2012) and family history of colorectal cancer (yes or no). The estimates and 95% CI of the model coefficients are given in Table C.10. After fitting the model, we removed the individuals who were diagnosed with CRC but their MSI subtype was missing, and used weights based on the inverse of the logistic regression predictions to adjust our sample. We truncated the weights at the 99% percentile to reduce variance. The distribution of the final weights used are given in

C.3 Additional results
Turning to the analysis results. First, we estimated the effects under S-Monotonicity for both disease subtypes. The results for the SF -ACE D are given in Table 2 of the main text. Here, we give the results for the SF -ACE RR in Table C.13, and the results for θ in Table C.14.
We then continued to the analysis under D-Monotonicity for both disease subtypes. The results for the SF -ACE D and θ D using the DR estimator were given in Figures 2 and 3 of the main text. Figures C.6 and C.7 give the corresponding results for standardization and IPTW methods. (1) RR estimator was relatively stable while SF -ACE (2) RR increased as a function of λ 2 and decreased as a function of λ 1 .
Finally, we present the results of an analysis after removing participants for which the last available data were before age 70, with and without inverse probability weighting to take them into account (Tables C.15