Eur J Epidemiol (2014) 29:813–820 DOI 10.1007/s10654-014-9953-9
METHODS
Bounds on sufficient-cause interaction Arvid Sjo¨lander • Woojoo Lee • Henrik Ka¨llberg Yudi Pawitan
•
Received: 7 May 2014 / Accepted: 17 September 2014 / Published online: 24 September 2014 Ó Springer Science+Business Media Dordrecht 2014
Abstract A common goal of epidemiologic research is to study how two exposures interact in causing a binary outcome. Sufficient-cause interaction is a special type of mechanistic interaction, which requires that two events (e.g. specific exposure levels from two risk factors) are necessary in order for the outcome to occur. Recently, tests have been derived to establish the presence of sufficientcause interactions, for categorical exposures with at most three levels. In this paper we derive prevalence bounds, i.e. lower and upper bounds on the prevalence of subjects for which sufficient-cause interaction is present. The derived bounds hold for categorical exposures with arbitrary many levels. We apply the bounds to data from a study of gene– gene interaction in the development of Rheumatoid Arthritis. We provide an R-program to estimate the bounds from real data. Keywords Causal inference Counterfactuals Interaction Potential outcomes
Electronic supplementary material The online version of this article (doi:10.1007/s10654-014-9953-9) contains supplementary material, which is available to authorized users. A. Sjo¨lander (&) Y. Pawitan Department of Medical Epidemiology, Karolinska Institute, Stockholm, Sweden e-mail: [email protected] W. Lee Department of Statistics, Inha University, Incheon, Korea H. Ka¨llberg Institute of Environmental Medicine, Karolinska Institute, Stockholm, Sweden
Introduction A common goal of epidemiologic research is to study how two exposures interact in causing a binary outcome. Rothman [1] proposed a conceptual model for causation, often referred to as the ‘sufficient-cause model’ or ‘causal pie model’ [2]. In the sufficient-cause model, the outcome is thought of as brought about by a set of events (e.g. specific exposure levels from two risk factors), referred to as ‘component causes’. A set of component causes is said to constitute a ‘sufficient cause’ of the outcome for a specific subject (or under a certain set of background conditions) if (a) the subject inevitably develops the outcome in the presence of this set of component causes, and (b) the set is minimal, so that the subject does not develop the outcome if the set is reduced (i.e. if not all events occur). If two component causes are parts of the same sufficient cause, so that there is a mechanism through which both components causes are necessary for the outcome to occur, then there is said to be sufficient-cause interaction between the two component causes. For instance, Ka¨llberg et al. [3] studied the gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE in the development of Rheumatoid Arthritis (RA). These genetic alleles are established risk factors for RA, and may thus be considered as ‘component causes’ of RA. If there are subjects who would develop RA when both alleles are present, and would not develop RA if either allele is absent (everything else equal), then these alleles together constitute a sufficient cause of RA, for these subjects. Furthermore, since these alleles are part of the same sufficient cause we would say there is sufficient-cause interaction between the alleles, for these subjects. We return to this example below.
123
A. Sjo¨lander et al.
814
The sufficient-cause model has received significant interest in recent years. VanderWeele and Robins [4, 5] related the sufficient-cause model to the potential outcome model and derived various statistical tests for the presence of sufficient-cause interaction. VanderWeele [6] clarified the relation between statistical interaction (i.e. the presence of product terms in regression models) and sufficient-cause interaction. VanderWeele et al. [7] showed how marginal structural models can be used to test for the presence of sufficient-cause interactions. VanderWeele [8] generalized the definition of sufficient-cause interaction for scenarios with ternary exposures, and derived various statistical tests for the presence of sufficient-cause interaction in such scenarios. Most of the literature on sufficient-cause interactions has focused on hypothesis testing. Once the presence of sufficient-cause interaction has been tested and established, it would often be interesting to estimate its prevalence, i.e. the prevalence of individuals for which there is sufficient-cause interaction. Unfortunately, this prevalence is not generally identifiable [4]. In this paper we derive prevalence bounds on sufficient-cause interactions, that is, lower and upper bounds which are guaranteed to contain the true prevalence. Lower prevalence bounds were presented by VanderWeele et al. [7] for the case where both exposures are binary. We generalize the work by VanderWeele et al. [7] by considering categorical exposures with arbitrary many levels, and by presenting both lower and upper prevalence bounds. The paper is organized as follows. In ‘‘Notation and definitions’’ section we present basic notations and definitions. In ‘‘Bounds’’ section we present the prevalence bounds on sufficient-cause interaction. We demonstrate that the bounds can be made significantly narrower by making an assumption about monotone exposure effects, which is often a reasonable assumption. In ‘‘Relation to previously proposed statistical tests’’ section we show that the bounds can be used to test for the presence of sufficient-cause interaction, and we discuss how such a test relates to the tests previously proposed in the literature. In ‘‘Relative prevalence of sufficient-cause interaction’’ section we discuss the possibility to ‘standardize’ the prevalence of sufficient-cause interaction, to obtain a more meaningful measure when the outcome is rare. In ‘‘Estimation of the bounds’’ section we discuss how the bounds can be estimated using statistical regression models. In ‘‘Application’’ section we show an application of the bounds to data from the Ka¨llberg et al. [3] study of gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE.
123
Notation and definitions We let X and Z denote two categorical exposures of interest, with levels coded as integers 0; . . .; KX and 0; . . .; KZ , respectively. We let Y denote the binary (0/1) outcome of interest, and we let Yxz ðxÞ denote the potential outcome [9, 10] for subject x when exposed to levels X ¼ x and Z ¼ z. We make two definitions, which generalize the definitions given by VanderWeele [8] by allowing for categorical exposures with arbitrary many levels. Definition 1 We say that there is strong sufficient-cause interaction between exposure levels x and z for subject x if Yxz ðxÞ ¼ 1 but Yx0 z ðxÞ ¼ 0 for all x0 6¼ x and Yxz0 ðxÞ ¼ 0 for all z0 6¼ z. Definition 2 We say that there is weak sufficient-cause interaction between exposure levels x and z for subject x if Yxz ðxÞ ¼ 1 but Yx0 z ðxÞ ¼ 0 for some x0 6¼ x and Yxz0 ðxÞ ¼ 0 for some z0 6¼ z. We note that VanderWeele [8] made a distinction between (strong/weak) sufficient-cause interaction and (definite/weak) interdependence. However, he also showed that these two concepts are mathematically equivalent, so we ignore this distinction here. Strong sufficient-cause interaction between x and z implies that any change from level x and z would change the outcome from 1 to 0, whereas weak sufficient-cause interaction between x and some change from level z implies that there is some change from level x and some change from level z that would change the outcome from 1 to 0. For instance, suppose that KX ¼ KZ ¼ 2 and that the potential outcomes for a specific subject is given by Table 1. For this subject there is strong sufficient-cause interaction between x ¼ 2 and z ¼ 2, and there is weak sufficient-cause interaction between x ¼ 2 and z ¼ 2, between x ¼ 1 and z ¼ 0, and between x ¼ 1 and z ¼ 1. Strong and weak sufficient-cause interaction are equivalent when both X and Z are binary, i.e. when KX ¼ KZ ¼ 1. Whether sufficient-cause interaction is present or not depends on how the outcome is coded. To see this, suppose that we reverse the coding of the outcome in Table 1, thus obtaining the potential outcomes in Table 2. Under this reversed coding there no strong sufficient-cause interaction between any levels x and z, and there is weak sufficient-cause interaction between x ¼ 1 and z ¼ 2, between x ¼ 2 and z ¼ 1, and between x ¼ 2 and z ¼ 2. For instance, that two exposures exhibit sufficient-cause interaction with respect to death does not mean that they exhibit sufficient-cause interaction with respect to survival. The coding dependence of sufficient-cause interaction necessitates careful thinking and biological knowledge when coding the outcome and drawing conclusions about sufficient-cause interaction.
Bounds on sufficient-cause interaction Table 1 Potential outcomes for a fictitious subject
Table 2 Potential outcomes for a fictitious subject when the coding of the outcome is reversed
z
815
x 0
1
2
0
0
1
0
1
0
1
0
2
0
0
1
z
x 0
1
2
0
1
0
1
1
1
0
1
2
1
1
0
We let hxz and gxz be the proportion of subjects for which there is strong and weak sufficient-cause interaction, respectively, between x and z, i.e. the prevalence of strong/ weak sufficient-cause interaction between x and z. hxz and gxz are the target parameters. Unfortunately these parameters are non-identifiable, even if both X and Z are randomized. In the next section derive bounds on hxz and gxz . We end this section by noting that sufficient-cause interaction is a special case of causal interaction also referred to as ‘biologic’ interaction [2, 11]. Causal interaction is said to be present if there is at least one level x at which Z has an effect, and there is at least one level z at which X has an effect. Because sufficient-cause interaction is a special case of causal interaction, the prevalence of the former cannot be larger than the prevalence of the latter.
Bounds Assumption-free bounds We let pxz denote the proportion of subjects who would develop the outcome if everybody, contrary to fact, would be exposed to levels x and z; pxz Pr fYxz ðxÞ ¼ 1g. If both X and Z are randomized, then pxz is equal to the proportion of subjects who develop the outcome among those factually exposed to levels x and z; pxz ¼ Pr ðY ¼ 1jX ¼ x; Z ¼ zÞ. In observational studies, this equality does not generally hold due to confounding of the exposure-outcome association. In this section we consider the pxz ’s as known, and we express the bounds on hxz and gxz in terms of these. Below, we discuss how the pxz ’s can be consistently estimated in observational studies, under an assumption of no unmeasured confounding. It can be shown (see online supplementary material) that the following inequalities always hold.
(
)
0
P px0 z z0 6¼z pxz0 9 8 > = < pxz > 1 p x0 z hxz 0 min0 > x 6¼x;z 6¼z> ; : 1 pxz0
max
pxz
P
x0 6¼x
ð1Þ
and max
x0 6¼x;z0 6¼z
0 pxz px0 z pxz0
gxz min
8 > =
pxz
x0 6¼x ð1 px0 z Þ : > ; : P ð1 p 0 Þ > z0 6¼z
xz
ð2Þ We let lðhxz Þ and uðhxz Þ denote the lower and upper bound on hxz in Eq. (1), respectively, and we let lðgxz Þ and uðgxz Þ denote the lower and upper bound on gxz in Eq. (2), respectively. When both exposures are binary (KX ¼ KZ ¼ 1) and x ¼ z ¼ 1, the bounds in Eqs. (1) and (2) simplify to 9 8 p11 > > = < 0 max h11 ¼ g11 min 1 p01 : > > p11 p01 p10 ; : 1 p10 ð3Þ The lower bound in Eq. (3) was reported by VanderWeele et al. [7]. Bounds under monotone exposure effects Sometimes, it may be reasonable to assume that the exposure effects are monotone, in the sense that the exposures are not causative for some subjects but preventative for other. Formally we define monotonicity as
123
A. Sjo¨lander et al.
816
Yxz ðxÞ Yx0 z0 ðxÞ if
x x0 and z z0
for all x:
ð4Þ
Under the monotonicity assumption in Eq. (4) we cannot have strong sufficient-cause interaction for any other x and z than x ¼ KX and z ¼ KZ , since monotonicity rules out that Yx0 z0 ðxÞ ¼ 0 for x0 [ x and z0 [ z. However, weak sufficient-cause interaction is still possible for all x and z. We define /xx0 zz0 ¼ pxz px0 z pxz0 þ px0 z0 . Under the monotonicity assumption it can be shown (see online supplementary material) that the following inequalities hold. (
)
0
h KX KZ /KX ðKX 1ÞKZ ðKZ 1Þ ( ) pKX KZ pðKX 1ÞKZ min pKX KZ pKX ðKZ 1Þ
max
ð5Þ
and 8 >
=
/xx0 zz0 > ; /xx00 zz00 /x0 x00 z0 z00 pxz p0z gxz min : pxz px0 max
x [ x0 [ x00 ;z [ z0 [ z00 > :
ð6Þ
We let lm ðhxz Þ and um ðhxz Þ denote the lower and upper bound on hxz in Eq. (5), respectively, and we let lm ðgxz Þ and um ðgxz Þ denote the lower and upper bound on gxz in Eq. (6), respectively. When both exposures are binary (KX ¼ KZ ¼ 1) and x ¼ z ¼ 1, the bounds in Eqs. (5) and (6) simplify to max
0
p11 p01 p10 þ p00 p11 p01 min : p11 p10
h11 ¼ g11 ð7Þ
The lower bound in Eq. (7) was reported by VanderWeele et al. [7]. The monotonicity assumption can lead to a dramatic improvement of the bounds. As an example, suppose that KX ¼ KZ ¼ 1, p00 ¼ p10 ¼ 0, and p11 ¼ p01 ¼ 0:5. In this example the bounds in Eq. (3) are equal to (0,0.5), whereas the bounds in Eq. (7) are equal to (0,0). Thus, in this example the lower bound coincides with the upper bound under the monotonicity assumption, so that h ¼ 0 becomes identifiable. More generally, we can characterize the situations where the bounds coincide, when KX ¼ KZ ¼ 1. The bounds in Eq. (3) coincide if either p11 ¼ 0, p01 ¼ 1, p10 ¼ 1, p11 ¼ 1 p01 ¼ 1, or p11 ¼ 1 p10 ¼ 1. The bounds in Eq. (7) coincide if either p11 ¼ p01 , p11 ¼ p10 , p00 ¼ p01 , or p00 ¼ p10 .
123
Relation to previously proposed statistical tests The bounds can be used to test for the presence of sufficient-cause interaction; if a particular lower bound is [ 0, then the corresponding sufficient-cause interaction is present. VanderWeele and Robins [5] and VanderWeele [8] have previously proposed various statistical tests for sufficient-cause interactions, for binary exposures, and for combinations of binary and ternary exposures, respectively. All these previous tests are on the form T [ 0, where the test statistic T is some linear combination of probabilities. However, VanderWeele and Robins [5] and VanderWeele [8] made no explicit reference to bounds, and their proposed test statistics have no obvious relation to the prevalence of sufficient-cause interaction. It turns out though, that all proposed test statistics by VanderWeele and Robins [5] and VanderWeele [8] can be identified as left-hand side elements of Eqs. (1), (2), (5), and (6). Thus, all the proposed test statistics can be interpreted as lower bounds on sufficient-cause interaction. For binary exposures, this was proven by VanderWeele et al. [7]. The authors further stated that ‘Similar remarks hold concerning bounds on the prevalence of sufficient-cause interactions when the exposures are categorical or ordinal’, citing VanderWeele [8]. However, VanderWeele [8] did not mention that his proposed test statistics can be interpreted as lower bounds, and VanderWeele et al. [7] did not prove their statement. To the best of our knowledge, it has not been proven elsewhere either. The tests proposed by VanderWeele and Robins [5] and VanderWeele [8] are implied by the bounds in Eqs. (1), (2), (5), and (6). The converse does not hold though; we give here two examples where our bounds provide more information, and are thus more powerful, than the previously proposed tests. The first example is when testing for weak sufficientcause interaction between x ¼ 1 and z ¼ 2, with KX ¼ 1 and KZ ¼ 2, under assumed monotonicity. For this case, VanderWeele [8] proposed the test /1020 [ 0 (Proposition 1 in VanderWeele [8]). However, from Eq. 6 we get the lower bound lm ðg12 Þ ¼ maxð/1020 ; /1021 Þ. For instance, suppose that p12 ¼ 1, p02 ¼ p11 ¼ p10 ¼ 0:75, and p00 ¼ 0. For these figures we obtain /1020 ¼ 1 0:75 0:75 ¼ 0:5, so that the proposed test fails to detect sufficient-cause interaction, but lm ðg12 Þ ¼ /1021 ¼ 1 0:75 0:75 þ 0:75 ¼ 0:25, so that the lower bound is [ 0. The second example is when testing for weak sufficient-cause interaction between x ¼ 2 and z ¼ 2, with KX ¼ KZ ¼ 2, under assumed monotonicity. For this case, VanderWeele [8] proposed the test max ð/2121 ; /2021 ; /2120 ; /2020 Þ [ 0 (Proposition 4 in VanderWeele [8]). However, from Eq. 6 we get the lower bound lm ðg22 Þ ¼ maxð/2121 ; /2021 ; /2120 ; /2020 ; /2020 /1010 Þ. For instance, suppose that p22 ¼ p21 ¼ p12 ¼ 1, p20 ¼ p02 ¼ p11 ¼ p10 ¼ p01 ¼ 0:75, and p00 ¼ 0. For these figures we obtain
Bounds on sufficient-cause interaction
/2121 ¼ 0:25, /2021 ¼ /2120 ¼ 0, and /2020 ¼ 0:5, so that the proposed test fails to detect sufficient-cause interaction, but lm ðh22 Þ ¼ /2020 /1010 ¼ 0:5 ð0:75Þ ¼ 0:25, so that the lower bound is [ 0. We finally note that the tests proposed by VanderWeele and Robins [5] and VanderWeele [8] only apply to exposures with at most three levels, whereas the bounds in Eqs. (1), (2), (5), and (6) can be used to test for sufficient-cause interaction between exposures with arbitrary many levels.
Relative prevalence of sufficient-cause interaction The upper bounds in Eqs. (1) and (2) are close to 0 if pxz is close to 0, which implies that the prevalence of sufficientcause interaction between x and z is close to 0 as well. This is reasonable; if the outcome is rare, then the total exposure effects (on the additive scale) cannot be large, and there is no room for large (additive) interactions either. Thus, if the outcome is rare we would always expect the prevalence of sufficient-cause interaction to be close to 0. The question arises as to whether it is possible to ‘standardize’ the prevalence of sufficient-cause interaction by a meaningful number, to obtain a more informative measure of interaction when the outcome is rare. One option is to standardize by pxz , which by definition is equal to the proportion of subjects for which the outcome is equal to 1 under levels x and z. Thus, hxz =pxz is equal to the proportion of subjects for which any change from levels x and z would change the outcome from 1 to 0, among those subjects for which the outcome is equal to 1 under levels x and z. Similarly, gxz =pxz is equal to the proportion of subjects for which there is some change from levels x and z that would change the outcome from 1 to 0, among those subjects for which the outcome is equal to 1 under levels x and z. We refer to hxz =pxz and gxz =pxz as the ‘relative’ prevalence of strong and weak sufficient-cause interaction, respectively. By dividing the bounds on hxz and gxz with pxz , we obtain bounds on the relative prevalence of sufficient-cause interaction.
Estimation of the bounds To estimate the bounds we first need to estimate the pxz ’s. In observational studies, this requires appropriate control for confounding. We let C denote a set of measured covariates. We say that C is sufficient for confounding control if the potential outcome Yxz ðxÞ is conditionally independent of ðX; ZÞ, given C [4]. Under this condition, we have that Pr fYxz ðxÞ ¼ 1jCg ¼ Pr ðY ¼ 1jX ¼ x; Z ¼ z; CÞ, from which it follows that pxz ¼ EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. If C is categorical with few levels, we
817
may estimate Pr ðY ¼ 1jX ¼ x; Z ¼ z; CÞ non-parametrically for each level of C, and then average over the sample distribution of C to obtain an estimate of EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. In real studies, C may often contain many covariates, of which some may be continuous or categorical with many levels. In such cases, we may use a parametric model for PrðY ¼ 1jX; Z; CÞ, e.g. a logistic regression model. We can fit this model through maximum likelihood, and we can use the fitted model to estimate PrðY ¼ 1jX ¼ x; Z ¼ z; CÞ for each observed value of C (i.e. for each subject in the data set). Finally, we can averaged these estimates to yield an estimate of EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. Alternatively, we can use inverse probability weighting, which requires a regression model for the conditional distribution of X and Z, given C [7].
Application In this section we show an application of the bounds on sufficient-cause interaction to real data. Ka¨llberg et al. [3] studied the gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE in the development of RA, using pooled data from three large case–control studies. The authors coded PTPN22 R620W as binary (none/any) and HLADRB1 SE as either binary (none/any allele) or ternary (none/ single/double allele). They defined the outcome as anti-CCP positive RA (a genetically distinct sub phenotype of RA which is characterized by presence of antibodies toward citrullinated peptides), and found a statistically significant (p value\0:001) additive interaction, for both the binary and the ternary coding of HLA-DRB1 SE. The authors concluded that ‘the observed interaction thus indicates the existence of a disease mechanism for anti-CCP positive RA that requires the simultaneous presence of the HLA-DRB1 SE alleles and the PTPN22 R620W allele’. We re-analyzed data from one of the case–control studies, the Swedish Epidemiological Investigation of Rheumatoid Arthritis (EIRA) study. Originally, these data comprised 1,183 cases and 793 controls, but were later extended to 1,426 cases and 868 controls; we here analyzed the extended data. The allele distribution across cases and controls is given by Table 3. In their statistical analysis, Ka¨llbert et al. [3] controlled for sex and age. Arguably, age is not a confounder in this setting, since genetic alleles are not associated with age. However, for illustrational purposes, and to facilitate a comparison with the results in Ka¨llberg et al. [3], we controlled for sex and age as well. We used a logistic regression model for RA, including main effects of PTPN22 R620W, HLA-DRB1 SE, sex, and age, as well as all possible interactions terms between PTPN22 R620W and HLA-DRB1 SE. We only
123
A. Sjo¨lander et al.
818 Table 3 EIRA case–control data on PTPN22 R620W and HLADRB1 SE PTPN22
HLA-DRB1
No. of cases
No. of controls
None
None
271
327
None
Single
509
276
Table 4 Bounds for strong sufficient-cause interaction between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, not assuming monotonicity x
z
hxz =pxz
hxz Bounds
Confidence interval
Bounds
Confidence interval
None
Double
239
67
Any
None
107
117
0
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
Any
Single
207
66
0
1
0.00,0.01
0.00,0.01
0.00,1.00
0.00,1.00
Any
Double
93
15
0
2
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
Sum = 868
1
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
1
1
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
1
2
0.00,0.02
0.00,0.04
0.00,1.00
0.00,1.00
Sum = 1,426
The case refers to anti-CCP positive Rheumatoid Arthritis
considered a ternary coding of HLA-DRB1 SE. We used weighted logistic regression, in which each case was given the weight p=p and each control was given the weight ð1 pÞ=ð1 p Þ, where p and p is the population and sample prevalence of RA, respectively. These weights account for the oversampling of cases due to the case–control sampling scheme. From the figures in Table 3 we have that p ¼ 1;426=ð1;426 þ 868Þ ¼ 0:62. The population prevalence of RA varies from population to population; we used the population prevalence of RA in United States, which is well documented and approximately equal to 0.006 [12]. In case–control settings where the disease prevalence is less well documented we would recommend to carry out a sensitivity analysis, by varying the prevalence over a range of plausible values and estimate the bounds for each value separately. From the weighted logistic regression model we obtained estimates of the pxz ’s, which were subsequently plugged into Eqs. (1), (2), (5), and (6) to obtain estimated bounds on the prevalence of strong and weak sufficient-cause interaction, with and without the monotonicity assumption, for each combination of x (level of PTPN22 R620W) and z (level of HLA-DRB1 SE). We also estimated bounds on the relative prevalence of strong and weak sufficient-cause interaction, as described above. To address sampling variability we constructed confidence intervals for the bounds as follows. For each estimated lower bound and for each estimated upper bound we calculated a 95 % bootstrap confidence interval, based on 1000 bootstrap replications, using the percentile method [13]. We constructed a confidence interval for the bounds by joining the lower 95 % confidence limit for the lower bound with the upper 95 % confidence limit for the upper bound. It is easy to show that a confidence interval constructed in this way has at least 95 % probability of covering the interval defined by the true bounds. Tables 4 and 5 display the estimated bounds for strong and weak sufficient-cause interaction, respectively, not assuming monotonicity. The results for strong and weak sufficient-cause interaction are very similar. The lower
123
Table 5 Bounds for weak sufficient-cause interaction, between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, not assuming monotonicity x
z
gxz =pxz
gxz Bounds
Confidence interval
Bounds
Confidence interval
0
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
0
1
0.00,0.01
0.00,0.01
0.00,1.00
0.00,1.00
0
2
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
1
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
1 1
1 2
0.00,0.01 0.01,0.02
0.00,0.02 0.00,0.04
0.15,1.00 0.26,1.00
0.00,1.00 0.00,1.00
bounds on hxz and gxz are equal to 0.00 for all allele combinations, except for the lower bound on gxz at x ¼ 1; z ¼ 2 which is equal to 0.01. However, at x ¼ 1; z ¼ 2 the lower confidence limit for gxz is equal to 0.00 as well. Thus, the presence of sufficient-cause interaction cannot be ruled in for any allele combination. The upper bound is close to 0 for all allele combinations, which means that if sufficient-cause interaction is present, then its prevalence is most likely small. This is not surprising, given that RA is a rare disease. The bounds on the hxz =pxz and gxz =pxz are quite uninformative, ranging from 0.00 to 1.00 at most levels x and z. Under monotonicity, the only allele combination at which there can be strong sufficient-cause interaction is x ¼ 1; z ¼ 2. The estimated bounds on h12 are equal to 0.00,0.01, with confidence interval equal to 0.00,0.03. The estimated bounds on h12 =p12 are equal to 0.19,0.42 with confidence interval equal to 0.00,0.69. Table 6 displays the estimated bounds for weak sufficient-cause interaction, assuming monotonicity. The upper bounds on gxz are all close to 0. However, at x ¼ 1 the lower bounds on gxz =pxz are all far from 0, as well as the corresponding lower confidence limits.
Bounds on sufficient-cause interaction
819
Table 6 Bounds for weak sufficient-cause interaction between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, assuming monotonicity x
z
gxz =pxz
gxz Bounds
Confidence interval
Bounds
Confidence interval
0
0
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
0
1
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
0
2
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
1
0
0.01,0.01
0.00,0.03
0.40,0.42
0.19,0.72
1
1
0.01,0.01
0.00,0.03
0.40,0.42
0.19,0.69
1
2
0.01,0.01
0.00,0.03
0.40,0.42
0.20,0.70
In conclusion, under monotonicity assumption, we see evidence of sufficient-cause interaction between the risk allele of PTPN22 gene with all the alleles of HLA-DRB1. This interaction is specific, in the sense that it does not occur with the non-risk alleles of PTPN22 gene, and accounts for a substantial proportion of the observed risk of RA due to the combination of the two genes. Sjo¨lander et al. [11] estimated bounds on the prevalence of causal interaction between PTPN22 R620W and HLADRB1 SE in the development of RA, using the same data as we have used here. These authors found a statistically significant presence of causal interaction, regardless of whether they assumed monotonicity or not. Given that the prevalence of sufficient-cause interaction is generally smaller than the prevalence of causal interaction, the power to detect sufficient-cause interaction will generally be smaller than the power to detect causal interaction. It is therefore not surprising that the presence of causal interaction did turn out statistically significant, whereas the presence of sufficientcause interaction did not in most cases.
Discussion In this paper we have generalized the definitions of strong and weak sufficient-cause interaction to cover scenarios with multi-level categorical exposures. We have derived bounds on the prevalence of sufficient-cause interaction, and we have shown that the bounds can be improved upon by making an assumption about monotone exposure effects. We have discussed how the bounds relate to previously proposed statistical tests, and we have shown how the bounds can be estimated from data using regression models. Estimation of the bounds requires some programing. We have implemented an R-function, bounds.suff, which estimates the bounds using a specified logistic regression model for the outcome, fitted through maximum likelihood, to control for measured covariates. The R-function is
available at the web page http://www.meb.ki.se/personal/ arvsjo/. We describe its usage in the online supplementary material. In line with the literature [4–8] we have focused on the proportion of subjects for which there is strong or weak sufficient-cause interaction between two particular exposure levels x and z. Other possible target parameters are the proportion of subjects for which there is strong and weak sufficient-cause interaction, respectively, between any x 2 ð0; . . .; KX Þ and z 2 ð0; . . .; KZ Þ, say h and g. We can show that the bounds on hxz and gxz can be used to bound h and g as well. Without assuming monotonicity we have that maxx;z flðhxz Þg P P h xz uðhxz Þ and maxx;z flðgxz Þg g xz uðgxz Þ. Under monotonicity we have that maxx;z flm ðhxz Þg h P P xz um ðhxz Þ and maxx;z flm ðgxz Þg g xz um ðgxz Þ. It is possible that these simple bounds can be made sharper, and we recognize this as an interesting topic for future research. Although the prevalence of sufficient-cause interaction is generally not identifiable, there may be special circumstances under which it can be identified. Liao and Lee [14] and Lee [15] demonstrated that the ‘causal pie weight’ and the ‘completions potential of sufficient component causes’, which are strongly related to the prevalence of sufficientcause interaction, can both be identified under certain parametric assumptions which may sometimes be reasonable. One important direction for future research is to explore whether there are plausible assumptions under which the prevalence of sufficient-cause interaction can be identified. Acknowledgments Arvid Sjo¨lander acknowledges financial support from The Swedish Research Council (340-2012-6007). Woojoo Lee acknowledges financial support from Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF2013R1A1A1061332). Conflict of interest of interest.
The authors declare that they have no conflict
References 1. Rothman K. Causes. Am J Epidemiol. 1976;104(6):587–92. 2. Rothman K, Greenland S, Lash T. Modern epidemiology. 3rd ed. Philadelphia: Lippincott Williams & Wilkins; 2008. 3. Ka¨llberg H, Padyukov L, Plenge R, Ro¨nnelid J, Gregersen P, van der Helm-van A, et al. Gene–gene and gene–environment interactions involving HLA-DRB1, PTPN22, and smoking in two subsets of rheumatoid arthritis. Am J Hum Genet. 2007;80(5):867. 4. VanderWeele T, Robins J. The identification of synergism in the sufficient-component-cause framework. Epidemiology. 2007; 18(3):329–39. 5. VanderWeele T, Robins J. Empirical and counterfactual conditions for sufficient cause interactions. Biometrika. 2008;95(1):49–61.
123
820 6. VanderWeele T. Sufficient cause interactions and statistical interactions. Epidemiology. 2009;20(1):6–13. 7. VanderWeele T, Vansteelandt S, Robins J. Marginal structural models for sufficient cause interactions. Am J Epidemiol. 2010; 171(4):506–14. 8. Vanderweele T. Sufficient cause interactions for categorical and ordinal exposures with three levels. Biometrika. 2010;97(3): 647–59. 9. Rubin D. Estimating causal effects of treatments in randomized and nonrandomized studies. J Educ Psychol. 1974;66(5): 688–701. 10. Pearl J. Causality: models, reasoning, and inference. 2nd ed. New York: Cambridge University Press; 2009.
123
A. Sjo¨lander et al. 11. Sjo¨lander A, Lee W, Ka¨llberg H, Pawitan Y. Bounds on causal interactions for binary outcomes. Biometrics. 2014. doi:10.1111/ biom.12166. 12. Helmick C, Felson D, Lawrence R, Gabriel S, Hirsch R, Kwoh C, Liang M, Kremers H, Mayes M, Merkel P, et al. Estimates of the prevalence of arthritis and other rheumatic conditions in the united states: part I. Arthritis Rheum. 2008;58(1):15–25. 13. Efron B, Tibshirani R. An introduction to the bootstrap. New York: Chapman & Hall; 1993. 14. Liao S, Lee W. Weighing the causal pies in case–control studies. Ann Epidemiol. 2010;20(7):568–73. 15. Lee W. Completion potentials of sufficient component causes. Epidemiology. 2012;23(3):446–53.
METHODS
Bounds on sufficient-cause interaction Arvid Sjo¨lander • Woojoo Lee • Henrik Ka¨llberg Yudi Pawitan
•
Received: 7 May 2014 / Accepted: 17 September 2014 / Published online: 24 September 2014 Ó Springer Science+Business Media Dordrecht 2014
Abstract A common goal of epidemiologic research is to study how two exposures interact in causing a binary outcome. Sufficient-cause interaction is a special type of mechanistic interaction, which requires that two events (e.g. specific exposure levels from two risk factors) are necessary in order for the outcome to occur. Recently, tests have been derived to establish the presence of sufficientcause interactions, for categorical exposures with at most three levels. In this paper we derive prevalence bounds, i.e. lower and upper bounds on the prevalence of subjects for which sufficient-cause interaction is present. The derived bounds hold for categorical exposures with arbitrary many levels. We apply the bounds to data from a study of gene– gene interaction in the development of Rheumatoid Arthritis. We provide an R-program to estimate the bounds from real data. Keywords Causal inference Counterfactuals Interaction Potential outcomes
Electronic supplementary material The online version of this article (doi:10.1007/s10654-014-9953-9) contains supplementary material, which is available to authorized users. A. Sjo¨lander (&) Y. Pawitan Department of Medical Epidemiology, Karolinska Institute, Stockholm, Sweden e-mail: [email protected] W. Lee Department of Statistics, Inha University, Incheon, Korea H. Ka¨llberg Institute of Environmental Medicine, Karolinska Institute, Stockholm, Sweden
Introduction A common goal of epidemiologic research is to study how two exposures interact in causing a binary outcome. Rothman [1] proposed a conceptual model for causation, often referred to as the ‘sufficient-cause model’ or ‘causal pie model’ [2]. In the sufficient-cause model, the outcome is thought of as brought about by a set of events (e.g. specific exposure levels from two risk factors), referred to as ‘component causes’. A set of component causes is said to constitute a ‘sufficient cause’ of the outcome for a specific subject (or under a certain set of background conditions) if (a) the subject inevitably develops the outcome in the presence of this set of component causes, and (b) the set is minimal, so that the subject does not develop the outcome if the set is reduced (i.e. if not all events occur). If two component causes are parts of the same sufficient cause, so that there is a mechanism through which both components causes are necessary for the outcome to occur, then there is said to be sufficient-cause interaction between the two component causes. For instance, Ka¨llberg et al. [3] studied the gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE in the development of Rheumatoid Arthritis (RA). These genetic alleles are established risk factors for RA, and may thus be considered as ‘component causes’ of RA. If there are subjects who would develop RA when both alleles are present, and would not develop RA if either allele is absent (everything else equal), then these alleles together constitute a sufficient cause of RA, for these subjects. Furthermore, since these alleles are part of the same sufficient cause we would say there is sufficient-cause interaction between the alleles, for these subjects. We return to this example below.
123
A. Sjo¨lander et al.
814
The sufficient-cause model has received significant interest in recent years. VanderWeele and Robins [4, 5] related the sufficient-cause model to the potential outcome model and derived various statistical tests for the presence of sufficient-cause interaction. VanderWeele [6] clarified the relation between statistical interaction (i.e. the presence of product terms in regression models) and sufficient-cause interaction. VanderWeele et al. [7] showed how marginal structural models can be used to test for the presence of sufficient-cause interactions. VanderWeele [8] generalized the definition of sufficient-cause interaction for scenarios with ternary exposures, and derived various statistical tests for the presence of sufficient-cause interaction in such scenarios. Most of the literature on sufficient-cause interactions has focused on hypothesis testing. Once the presence of sufficient-cause interaction has been tested and established, it would often be interesting to estimate its prevalence, i.e. the prevalence of individuals for which there is sufficient-cause interaction. Unfortunately, this prevalence is not generally identifiable [4]. In this paper we derive prevalence bounds on sufficient-cause interactions, that is, lower and upper bounds which are guaranteed to contain the true prevalence. Lower prevalence bounds were presented by VanderWeele et al. [7] for the case where both exposures are binary. We generalize the work by VanderWeele et al. [7] by considering categorical exposures with arbitrary many levels, and by presenting both lower and upper prevalence bounds. The paper is organized as follows. In ‘‘Notation and definitions’’ section we present basic notations and definitions. In ‘‘Bounds’’ section we present the prevalence bounds on sufficient-cause interaction. We demonstrate that the bounds can be made significantly narrower by making an assumption about monotone exposure effects, which is often a reasonable assumption. In ‘‘Relation to previously proposed statistical tests’’ section we show that the bounds can be used to test for the presence of sufficient-cause interaction, and we discuss how such a test relates to the tests previously proposed in the literature. In ‘‘Relative prevalence of sufficient-cause interaction’’ section we discuss the possibility to ‘standardize’ the prevalence of sufficient-cause interaction, to obtain a more meaningful measure when the outcome is rare. In ‘‘Estimation of the bounds’’ section we discuss how the bounds can be estimated using statistical regression models. In ‘‘Application’’ section we show an application of the bounds to data from the Ka¨llberg et al. [3] study of gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE.
123
Notation and definitions We let X and Z denote two categorical exposures of interest, with levels coded as integers 0; . . .; KX and 0; . . .; KZ , respectively. We let Y denote the binary (0/1) outcome of interest, and we let Yxz ðxÞ denote the potential outcome [9, 10] for subject x when exposed to levels X ¼ x and Z ¼ z. We make two definitions, which generalize the definitions given by VanderWeele [8] by allowing for categorical exposures with arbitrary many levels. Definition 1 We say that there is strong sufficient-cause interaction between exposure levels x and z for subject x if Yxz ðxÞ ¼ 1 but Yx0 z ðxÞ ¼ 0 for all x0 6¼ x and Yxz0 ðxÞ ¼ 0 for all z0 6¼ z. Definition 2 We say that there is weak sufficient-cause interaction between exposure levels x and z for subject x if Yxz ðxÞ ¼ 1 but Yx0 z ðxÞ ¼ 0 for some x0 6¼ x and Yxz0 ðxÞ ¼ 0 for some z0 6¼ z. We note that VanderWeele [8] made a distinction between (strong/weak) sufficient-cause interaction and (definite/weak) interdependence. However, he also showed that these two concepts are mathematically equivalent, so we ignore this distinction here. Strong sufficient-cause interaction between x and z implies that any change from level x and z would change the outcome from 1 to 0, whereas weak sufficient-cause interaction between x and some change from level z implies that there is some change from level x and some change from level z that would change the outcome from 1 to 0. For instance, suppose that KX ¼ KZ ¼ 2 and that the potential outcomes for a specific subject is given by Table 1. For this subject there is strong sufficient-cause interaction between x ¼ 2 and z ¼ 2, and there is weak sufficient-cause interaction between x ¼ 2 and z ¼ 2, between x ¼ 1 and z ¼ 0, and between x ¼ 1 and z ¼ 1. Strong and weak sufficient-cause interaction are equivalent when both X and Z are binary, i.e. when KX ¼ KZ ¼ 1. Whether sufficient-cause interaction is present or not depends on how the outcome is coded. To see this, suppose that we reverse the coding of the outcome in Table 1, thus obtaining the potential outcomes in Table 2. Under this reversed coding there no strong sufficient-cause interaction between any levels x and z, and there is weak sufficient-cause interaction between x ¼ 1 and z ¼ 2, between x ¼ 2 and z ¼ 1, and between x ¼ 2 and z ¼ 2. For instance, that two exposures exhibit sufficient-cause interaction with respect to death does not mean that they exhibit sufficient-cause interaction with respect to survival. The coding dependence of sufficient-cause interaction necessitates careful thinking and biological knowledge when coding the outcome and drawing conclusions about sufficient-cause interaction.
Bounds on sufficient-cause interaction Table 1 Potential outcomes for a fictitious subject
Table 2 Potential outcomes for a fictitious subject when the coding of the outcome is reversed
z
815
x 0
1
2
0
0
1
0
1
0
1
0
2
0
0
1
z
x 0
1
2
0
1
0
1
1
1
0
1
2
1
1
0
We let hxz and gxz be the proportion of subjects for which there is strong and weak sufficient-cause interaction, respectively, between x and z, i.e. the prevalence of strong/ weak sufficient-cause interaction between x and z. hxz and gxz are the target parameters. Unfortunately these parameters are non-identifiable, even if both X and Z are randomized. In the next section derive bounds on hxz and gxz . We end this section by noting that sufficient-cause interaction is a special case of causal interaction also referred to as ‘biologic’ interaction [2, 11]. Causal interaction is said to be present if there is at least one level x at which Z has an effect, and there is at least one level z at which X has an effect. Because sufficient-cause interaction is a special case of causal interaction, the prevalence of the former cannot be larger than the prevalence of the latter.
Bounds Assumption-free bounds We let pxz denote the proportion of subjects who would develop the outcome if everybody, contrary to fact, would be exposed to levels x and z; pxz Pr fYxz ðxÞ ¼ 1g. If both X and Z are randomized, then pxz is equal to the proportion of subjects who develop the outcome among those factually exposed to levels x and z; pxz ¼ Pr ðY ¼ 1jX ¼ x; Z ¼ zÞ. In observational studies, this equality does not generally hold due to confounding of the exposure-outcome association. In this section we consider the pxz ’s as known, and we express the bounds on hxz and gxz in terms of these. Below, we discuss how the pxz ’s can be consistently estimated in observational studies, under an assumption of no unmeasured confounding. It can be shown (see online supplementary material) that the following inequalities always hold.
(
)
0
P px0 z z0 6¼z pxz0 9 8 > = < pxz > 1 p x0 z hxz 0 min0 > x 6¼x;z 6¼z> ; : 1 pxz0
max
pxz
P
x0 6¼x
ð1Þ
and max
x0 6¼x;z0 6¼z
0 pxz px0 z pxz0
gxz min
8 > =
pxz
x0 6¼x ð1 px0 z Þ : > ; : P ð1 p 0 Þ > z0 6¼z
xz
ð2Þ We let lðhxz Þ and uðhxz Þ denote the lower and upper bound on hxz in Eq. (1), respectively, and we let lðgxz Þ and uðgxz Þ denote the lower and upper bound on gxz in Eq. (2), respectively. When both exposures are binary (KX ¼ KZ ¼ 1) and x ¼ z ¼ 1, the bounds in Eqs. (1) and (2) simplify to 9 8 p11 > > = < 0 max h11 ¼ g11 min 1 p01 : > > p11 p01 p10 ; : 1 p10 ð3Þ The lower bound in Eq. (3) was reported by VanderWeele et al. [7]. Bounds under monotone exposure effects Sometimes, it may be reasonable to assume that the exposure effects are monotone, in the sense that the exposures are not causative for some subjects but preventative for other. Formally we define monotonicity as
123
A. Sjo¨lander et al.
816
Yxz ðxÞ Yx0 z0 ðxÞ if
x x0 and z z0
for all x:
ð4Þ
Under the monotonicity assumption in Eq. (4) we cannot have strong sufficient-cause interaction for any other x and z than x ¼ KX and z ¼ KZ , since monotonicity rules out that Yx0 z0 ðxÞ ¼ 0 for x0 [ x and z0 [ z. However, weak sufficient-cause interaction is still possible for all x and z. We define /xx0 zz0 ¼ pxz px0 z pxz0 þ px0 z0 . Under the monotonicity assumption it can be shown (see online supplementary material) that the following inequalities hold. (
)
0
h KX KZ /KX ðKX 1ÞKZ ðKZ 1Þ ( ) pKX KZ pðKX 1ÞKZ min pKX KZ pKX ðKZ 1Þ
max
ð5Þ
and 8 >
=
/xx0 zz0 > ; /xx00 zz00 /x0 x00 z0 z00 pxz p0z gxz min : pxz px0 max
x [ x0 [ x00 ;z [ z0 [ z00 > :
ð6Þ
We let lm ðhxz Þ and um ðhxz Þ denote the lower and upper bound on hxz in Eq. (5), respectively, and we let lm ðgxz Þ and um ðgxz Þ denote the lower and upper bound on gxz in Eq. (6), respectively. When both exposures are binary (KX ¼ KZ ¼ 1) and x ¼ z ¼ 1, the bounds in Eqs. (5) and (6) simplify to max
0
p11 p01 p10 þ p00 p11 p01 min : p11 p10
h11 ¼ g11 ð7Þ
The lower bound in Eq. (7) was reported by VanderWeele et al. [7]. The monotonicity assumption can lead to a dramatic improvement of the bounds. As an example, suppose that KX ¼ KZ ¼ 1, p00 ¼ p10 ¼ 0, and p11 ¼ p01 ¼ 0:5. In this example the bounds in Eq. (3) are equal to (0,0.5), whereas the bounds in Eq. (7) are equal to (0,0). Thus, in this example the lower bound coincides with the upper bound under the monotonicity assumption, so that h ¼ 0 becomes identifiable. More generally, we can characterize the situations where the bounds coincide, when KX ¼ KZ ¼ 1. The bounds in Eq. (3) coincide if either p11 ¼ 0, p01 ¼ 1, p10 ¼ 1, p11 ¼ 1 p01 ¼ 1, or p11 ¼ 1 p10 ¼ 1. The bounds in Eq. (7) coincide if either p11 ¼ p01 , p11 ¼ p10 , p00 ¼ p01 , or p00 ¼ p10 .
123
Relation to previously proposed statistical tests The bounds can be used to test for the presence of sufficient-cause interaction; if a particular lower bound is [ 0, then the corresponding sufficient-cause interaction is present. VanderWeele and Robins [5] and VanderWeele [8] have previously proposed various statistical tests for sufficient-cause interactions, for binary exposures, and for combinations of binary and ternary exposures, respectively. All these previous tests are on the form T [ 0, where the test statistic T is some linear combination of probabilities. However, VanderWeele and Robins [5] and VanderWeele [8] made no explicit reference to bounds, and their proposed test statistics have no obvious relation to the prevalence of sufficient-cause interaction. It turns out though, that all proposed test statistics by VanderWeele and Robins [5] and VanderWeele [8] can be identified as left-hand side elements of Eqs. (1), (2), (5), and (6). Thus, all the proposed test statistics can be interpreted as lower bounds on sufficient-cause interaction. For binary exposures, this was proven by VanderWeele et al. [7]. The authors further stated that ‘Similar remarks hold concerning bounds on the prevalence of sufficient-cause interactions when the exposures are categorical or ordinal’, citing VanderWeele [8]. However, VanderWeele [8] did not mention that his proposed test statistics can be interpreted as lower bounds, and VanderWeele et al. [7] did not prove their statement. To the best of our knowledge, it has not been proven elsewhere either. The tests proposed by VanderWeele and Robins [5] and VanderWeele [8] are implied by the bounds in Eqs. (1), (2), (5), and (6). The converse does not hold though; we give here two examples where our bounds provide more information, and are thus more powerful, than the previously proposed tests. The first example is when testing for weak sufficientcause interaction between x ¼ 1 and z ¼ 2, with KX ¼ 1 and KZ ¼ 2, under assumed monotonicity. For this case, VanderWeele [8] proposed the test /1020 [ 0 (Proposition 1 in VanderWeele [8]). However, from Eq. 6 we get the lower bound lm ðg12 Þ ¼ maxð/1020 ; /1021 Þ. For instance, suppose that p12 ¼ 1, p02 ¼ p11 ¼ p10 ¼ 0:75, and p00 ¼ 0. For these figures we obtain /1020 ¼ 1 0:75 0:75 ¼ 0:5, so that the proposed test fails to detect sufficient-cause interaction, but lm ðg12 Þ ¼ /1021 ¼ 1 0:75 0:75 þ 0:75 ¼ 0:25, so that the lower bound is [ 0. The second example is when testing for weak sufficient-cause interaction between x ¼ 2 and z ¼ 2, with KX ¼ KZ ¼ 2, under assumed monotonicity. For this case, VanderWeele [8] proposed the test max ð/2121 ; /2021 ; /2120 ; /2020 Þ [ 0 (Proposition 4 in VanderWeele [8]). However, from Eq. 6 we get the lower bound lm ðg22 Þ ¼ maxð/2121 ; /2021 ; /2120 ; /2020 ; /2020 /1010 Þ. For instance, suppose that p22 ¼ p21 ¼ p12 ¼ 1, p20 ¼ p02 ¼ p11 ¼ p10 ¼ p01 ¼ 0:75, and p00 ¼ 0. For these figures we obtain
Bounds on sufficient-cause interaction
/2121 ¼ 0:25, /2021 ¼ /2120 ¼ 0, and /2020 ¼ 0:5, so that the proposed test fails to detect sufficient-cause interaction, but lm ðh22 Þ ¼ /2020 /1010 ¼ 0:5 ð0:75Þ ¼ 0:25, so that the lower bound is [ 0. We finally note that the tests proposed by VanderWeele and Robins [5] and VanderWeele [8] only apply to exposures with at most three levels, whereas the bounds in Eqs. (1), (2), (5), and (6) can be used to test for sufficient-cause interaction between exposures with arbitrary many levels.
Relative prevalence of sufficient-cause interaction The upper bounds in Eqs. (1) and (2) are close to 0 if pxz is close to 0, which implies that the prevalence of sufficientcause interaction between x and z is close to 0 as well. This is reasonable; if the outcome is rare, then the total exposure effects (on the additive scale) cannot be large, and there is no room for large (additive) interactions either. Thus, if the outcome is rare we would always expect the prevalence of sufficient-cause interaction to be close to 0. The question arises as to whether it is possible to ‘standardize’ the prevalence of sufficient-cause interaction by a meaningful number, to obtain a more informative measure of interaction when the outcome is rare. One option is to standardize by pxz , which by definition is equal to the proportion of subjects for which the outcome is equal to 1 under levels x and z. Thus, hxz =pxz is equal to the proportion of subjects for which any change from levels x and z would change the outcome from 1 to 0, among those subjects for which the outcome is equal to 1 under levels x and z. Similarly, gxz =pxz is equal to the proportion of subjects for which there is some change from levels x and z that would change the outcome from 1 to 0, among those subjects for which the outcome is equal to 1 under levels x and z. We refer to hxz =pxz and gxz =pxz as the ‘relative’ prevalence of strong and weak sufficient-cause interaction, respectively. By dividing the bounds on hxz and gxz with pxz , we obtain bounds on the relative prevalence of sufficient-cause interaction.
Estimation of the bounds To estimate the bounds we first need to estimate the pxz ’s. In observational studies, this requires appropriate control for confounding. We let C denote a set of measured covariates. We say that C is sufficient for confounding control if the potential outcome Yxz ðxÞ is conditionally independent of ðX; ZÞ, given C [4]. Under this condition, we have that Pr fYxz ðxÞ ¼ 1jCg ¼ Pr ðY ¼ 1jX ¼ x; Z ¼ z; CÞ, from which it follows that pxz ¼ EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. If C is categorical with few levels, we
817
may estimate Pr ðY ¼ 1jX ¼ x; Z ¼ z; CÞ non-parametrically for each level of C, and then average over the sample distribution of C to obtain an estimate of EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. In real studies, C may often contain many covariates, of which some may be continuous or categorical with many levels. In such cases, we may use a parametric model for PrðY ¼ 1jX; Z; CÞ, e.g. a logistic regression model. We can fit this model through maximum likelihood, and we can use the fitted model to estimate PrðY ¼ 1jX ¼ x; Z ¼ z; CÞ for each observed value of C (i.e. for each subject in the data set). Finally, we can averaged these estimates to yield an estimate of EfPrðY ¼ 1jX ¼ x; Z ¼ z; CÞg. Alternatively, we can use inverse probability weighting, which requires a regression model for the conditional distribution of X and Z, given C [7].
Application In this section we show an application of the bounds on sufficient-cause interaction to real data. Ka¨llberg et al. [3] studied the gene–gene interaction between PTPN22 R620W and HLA-DRB1 SE in the development of RA, using pooled data from three large case–control studies. The authors coded PTPN22 R620W as binary (none/any) and HLADRB1 SE as either binary (none/any allele) or ternary (none/ single/double allele). They defined the outcome as anti-CCP positive RA (a genetically distinct sub phenotype of RA which is characterized by presence of antibodies toward citrullinated peptides), and found a statistically significant (p value\0:001) additive interaction, for both the binary and the ternary coding of HLA-DRB1 SE. The authors concluded that ‘the observed interaction thus indicates the existence of a disease mechanism for anti-CCP positive RA that requires the simultaneous presence of the HLA-DRB1 SE alleles and the PTPN22 R620W allele’. We re-analyzed data from one of the case–control studies, the Swedish Epidemiological Investigation of Rheumatoid Arthritis (EIRA) study. Originally, these data comprised 1,183 cases and 793 controls, but were later extended to 1,426 cases and 868 controls; we here analyzed the extended data. The allele distribution across cases and controls is given by Table 3. In their statistical analysis, Ka¨llbert et al. [3] controlled for sex and age. Arguably, age is not a confounder in this setting, since genetic alleles are not associated with age. However, for illustrational purposes, and to facilitate a comparison with the results in Ka¨llberg et al. [3], we controlled for sex and age as well. We used a logistic regression model for RA, including main effects of PTPN22 R620W, HLA-DRB1 SE, sex, and age, as well as all possible interactions terms between PTPN22 R620W and HLA-DRB1 SE. We only
123
A. Sjo¨lander et al.
818 Table 3 EIRA case–control data on PTPN22 R620W and HLADRB1 SE PTPN22
HLA-DRB1
No. of cases
No. of controls
None
None
271
327
None
Single
509
276
Table 4 Bounds for strong sufficient-cause interaction between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, not assuming monotonicity x
z
hxz =pxz
hxz Bounds
Confidence interval
Bounds
Confidence interval
None
Double
239
67
Any
None
107
117
0
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
Any
Single
207
66
0
1
0.00,0.01
0.00,0.01
0.00,1.00
0.00,1.00
Any
Double
93
15
0
2
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
Sum = 868
1
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
1
1
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
1
2
0.00,0.02
0.00,0.04
0.00,1.00
0.00,1.00
Sum = 1,426
The case refers to anti-CCP positive Rheumatoid Arthritis
considered a ternary coding of HLA-DRB1 SE. We used weighted logistic regression, in which each case was given the weight p=p and each control was given the weight ð1 pÞ=ð1 p Þ, where p and p is the population and sample prevalence of RA, respectively. These weights account for the oversampling of cases due to the case–control sampling scheme. From the figures in Table 3 we have that p ¼ 1;426=ð1;426 þ 868Þ ¼ 0:62. The population prevalence of RA varies from population to population; we used the population prevalence of RA in United States, which is well documented and approximately equal to 0.006 [12]. In case–control settings where the disease prevalence is less well documented we would recommend to carry out a sensitivity analysis, by varying the prevalence over a range of plausible values and estimate the bounds for each value separately. From the weighted logistic regression model we obtained estimates of the pxz ’s, which were subsequently plugged into Eqs. (1), (2), (5), and (6) to obtain estimated bounds on the prevalence of strong and weak sufficient-cause interaction, with and without the monotonicity assumption, for each combination of x (level of PTPN22 R620W) and z (level of HLA-DRB1 SE). We also estimated bounds on the relative prevalence of strong and weak sufficient-cause interaction, as described above. To address sampling variability we constructed confidence intervals for the bounds as follows. For each estimated lower bound and for each estimated upper bound we calculated a 95 % bootstrap confidence interval, based on 1000 bootstrap replications, using the percentile method [13]. We constructed a confidence interval for the bounds by joining the lower 95 % confidence limit for the lower bound with the upper 95 % confidence limit for the upper bound. It is easy to show that a confidence interval constructed in this way has at least 95 % probability of covering the interval defined by the true bounds. Tables 4 and 5 display the estimated bounds for strong and weak sufficient-cause interaction, respectively, not assuming monotonicity. The results for strong and weak sufficient-cause interaction are very similar. The lower
123
Table 5 Bounds for weak sufficient-cause interaction, between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, not assuming monotonicity x
z
gxz =pxz
gxz Bounds
Confidence interval
Bounds
Confidence interval
0
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
0
1
0.00,0.01
0.00,0.01
0.00,1.00
0.00,1.00
0
2
0.00,0.01
0.00,0.02
0.00,1.00
0.00,1.00
1
0
0.00,0.00
0.00,0.00
0.00,1.00
0.00,1.00
1 1
1 2
0.00,0.01 0.01,0.02
0.00,0.02 0.00,0.04
0.15,1.00 0.26,1.00
0.00,1.00 0.00,1.00
bounds on hxz and gxz are equal to 0.00 for all allele combinations, except for the lower bound on gxz at x ¼ 1; z ¼ 2 which is equal to 0.01. However, at x ¼ 1; z ¼ 2 the lower confidence limit for gxz is equal to 0.00 as well. Thus, the presence of sufficient-cause interaction cannot be ruled in for any allele combination. The upper bound is close to 0 for all allele combinations, which means that if sufficient-cause interaction is present, then its prevalence is most likely small. This is not surprising, given that RA is a rare disease. The bounds on the hxz =pxz and gxz =pxz are quite uninformative, ranging from 0.00 to 1.00 at most levels x and z. Under monotonicity, the only allele combination at which there can be strong sufficient-cause interaction is x ¼ 1; z ¼ 2. The estimated bounds on h12 are equal to 0.00,0.01, with confidence interval equal to 0.00,0.03. The estimated bounds on h12 =p12 are equal to 0.19,0.42 with confidence interval equal to 0.00,0.69. Table 6 displays the estimated bounds for weak sufficient-cause interaction, assuming monotonicity. The upper bounds on gxz are all close to 0. However, at x ¼ 1 the lower bounds on gxz =pxz are all far from 0, as well as the corresponding lower confidence limits.
Bounds on sufficient-cause interaction
819
Table 6 Bounds for weak sufficient-cause interaction between PTPN22 R620W (X) and HLA-DRB1 SE (Z) in the development of Rheumatoid Arthritis, assuming monotonicity x
z
gxz =pxz
gxz Bounds
Confidence interval
Bounds
Confidence interval
0
0
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
0
1
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
0
2
0.00,0.00
0.00,0.00
0.00,0.00
0.00,0.00
1
0
0.01,0.01
0.00,0.03
0.40,0.42
0.19,0.72
1
1
0.01,0.01
0.00,0.03
0.40,0.42
0.19,0.69
1
2
0.01,0.01
0.00,0.03
0.40,0.42
0.20,0.70
In conclusion, under monotonicity assumption, we see evidence of sufficient-cause interaction between the risk allele of PTPN22 gene with all the alleles of HLA-DRB1. This interaction is specific, in the sense that it does not occur with the non-risk alleles of PTPN22 gene, and accounts for a substantial proportion of the observed risk of RA due to the combination of the two genes. Sjo¨lander et al. [11] estimated bounds on the prevalence of causal interaction between PTPN22 R620W and HLADRB1 SE in the development of RA, using the same data as we have used here. These authors found a statistically significant presence of causal interaction, regardless of whether they assumed monotonicity or not. Given that the prevalence of sufficient-cause interaction is generally smaller than the prevalence of causal interaction, the power to detect sufficient-cause interaction will generally be smaller than the power to detect causal interaction. It is therefore not surprising that the presence of causal interaction did turn out statistically significant, whereas the presence of sufficientcause interaction did not in most cases.
Discussion In this paper we have generalized the definitions of strong and weak sufficient-cause interaction to cover scenarios with multi-level categorical exposures. We have derived bounds on the prevalence of sufficient-cause interaction, and we have shown that the bounds can be improved upon by making an assumption about monotone exposure effects. We have discussed how the bounds relate to previously proposed statistical tests, and we have shown how the bounds can be estimated from data using regression models. Estimation of the bounds requires some programing. We have implemented an R-function, bounds.suff, which estimates the bounds using a specified logistic regression model for the outcome, fitted through maximum likelihood, to control for measured covariates. The R-function is
available at the web page http://www.meb.ki.se/personal/ arvsjo/. We describe its usage in the online supplementary material. In line with the literature [4–8] we have focused on the proportion of subjects for which there is strong or weak sufficient-cause interaction between two particular exposure levels x and z. Other possible target parameters are the proportion of subjects for which there is strong and weak sufficient-cause interaction, respectively, between any x 2 ð0; . . .; KX Þ and z 2 ð0; . . .; KZ Þ, say h and g. We can show that the bounds on hxz and gxz can be used to bound h and g as well. Without assuming monotonicity we have that maxx;z flðhxz Þg P P h xz uðhxz Þ and maxx;z flðgxz Þg g xz uðgxz Þ. Under monotonicity we have that maxx;z flm ðhxz Þg h P P xz um ðhxz Þ and maxx;z flm ðgxz Þg g xz um ðgxz Þ. It is possible that these simple bounds can be made sharper, and we recognize this as an interesting topic for future research. Although the prevalence of sufficient-cause interaction is generally not identifiable, there may be special circumstances under which it can be identified. Liao and Lee [14] and Lee [15] demonstrated that the ‘causal pie weight’ and the ‘completions potential of sufficient component causes’, which are strongly related to the prevalence of sufficientcause interaction, can both be identified under certain parametric assumptions which may sometimes be reasonable. One important direction for future research is to explore whether there are plausible assumptions under which the prevalence of sufficient-cause interaction can be identified. Acknowledgments Arvid Sjo¨lander acknowledges financial support from The Swedish Research Council (340-2012-6007). Woojoo Lee acknowledges financial support from Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF2013R1A1A1061332). Conflict of interest of interest.
The authors declare that they have no conflict
References 1. Rothman K. Causes. Am J Epidemiol. 1976;104(6):587–92. 2. Rothman K, Greenland S, Lash T. Modern epidemiology. 3rd ed. Philadelphia: Lippincott Williams & Wilkins; 2008. 3. Ka¨llberg H, Padyukov L, Plenge R, Ro¨nnelid J, Gregersen P, van der Helm-van A, et al. Gene–gene and gene–environment interactions involving HLA-DRB1, PTPN22, and smoking in two subsets of rheumatoid arthritis. Am J Hum Genet. 2007;80(5):867. 4. VanderWeele T, Robins J. The identification of synergism in the sufficient-component-cause framework. Epidemiology. 2007; 18(3):329–39. 5. VanderWeele T, Robins J. Empirical and counterfactual conditions for sufficient cause interactions. Biometrika. 2008;95(1):49–61.
123
820 6. VanderWeele T. Sufficient cause interactions and statistical interactions. Epidemiology. 2009;20(1):6–13. 7. VanderWeele T, Vansteelandt S, Robins J. Marginal structural models for sufficient cause interactions. Am J Epidemiol. 2010; 171(4):506–14. 8. Vanderweele T. Sufficient cause interactions for categorical and ordinal exposures with three levels. Biometrika. 2010;97(3): 647–59. 9. Rubin D. Estimating causal effects of treatments in randomized and nonrandomized studies. J Educ Psychol. 1974;66(5): 688–701. 10. Pearl J. Causality: models, reasoning, and inference. 2nd ed. New York: Cambridge University Press; 2009.
123
A. Sjo¨lander et al. 11. Sjo¨lander A, Lee W, Ka¨llberg H, Pawitan Y. Bounds on causal interactions for binary outcomes. Biometrics. 2014. doi:10.1111/ biom.12166. 12. Helmick C, Felson D, Lawrence R, Gabriel S, Hirsch R, Kwoh C, Liang M, Kremers H, Mayes M, Merkel P, et al. Estimates of the prevalence of arthritis and other rheumatic conditions in the united states: part I. Arthritis Rheum. 2008;58(1):15–25. 13. Efron B, Tibshirani R. An introduction to the bootstrap. New York: Chapman & Hall; 1993. 14. Liao S, Lee W. Weighing the causal pies in case–control studies. Ann Epidemiol. 2010;20(7):568–73. 15. Lee W. Completion potentials of sufficient component causes. Epidemiology. 2012;23(3):446–53.
