1 Introduction

Most complex diseases are influenced by a number of genetic and environmental factors, co-factors or predictors that either have a marginal effect and/or interact in some way (Maher 2008; Eichler et al. 2010; Bookman et al. 2011). The impact of these factors on disease risk is specified with a penetrance function

where Y equals 1 (0) for an individual with (without) the disease, and x=(x1,…,xp) is a vector with exposure status for p co-factors. We will assume that the factors are binary, with xi=1 or 0 depending on whether an individual is exposed to factor i or not, although extensions to continuous covariates will be discussed in Section 8. Often a few of the variables J⊂{1,…,p} are of main interest, whereas the others are confounders. In this paper we develop a unified theory for the attributable risk or attributable proportion (AP), and the population attributable fraction (PAF), when the effects of confounders are taken into account. This theory will involve both joint effects, marginal effects and interaction.

The joint effect of all non-confounders will refer to their cumulative marginal effects and interaction. The traditional definition of AP;

is the joint fraction of disease risk due the factors in J. The penetrance function θrem(x) is the disease risk of an individual whose exposure to all factors i∈J of interest has been removed, whereas the confounding variables i∉J remain at the levels specified by x. Equation (2) generalizes the definition of attributable proportion, first introduced by Levin (1953) to quantify the effect that one single risk factor, smoking, has on lung cancer. The joint effect is then equivalent to a marginal effect, since it is not possible to have interaction for one single factor. Greenland and Robins (1988) noted that sometimes there is an ambiguity as to whether θrem should include all subjects that develop the disease without exposure to the relevant factor(s), or only those for which the factor(s) has no causal role. Often the former definition of attributable proportion is used, since it is easier to estimate from data.

Since smoking is a risk factor (θ>θrem), APtrad is a number between 0 and 1. On the other hand, it is negative without any lower bound for preventive factors (θ<θrem). It is then possible to use the preventive fraction

a number between 0 and 1 that equals the fraction of reduced disease risk for a subject exposed to preventive factors, see for instance Aschengrau et al. (2003, pp. 65–66). But it is redundant to use two different quantities APtrad and PF in order to quantify the effect of one set of factors, in particular when it is not known before a study whether a predictor increases or decreases risk. The relative risk could replace eqs (2) and (3) but it has no upper bound for risk factors. Instead, we propose another normalization

of the attributable proportion, a number between –1 and 1 that equals APtrad for risk factors and −PF for preventive factors.

In order to define a version of AP that averages the effect of all confounders, while still controlling for their effect, we introduce θaddx, the penetrance function for a subject whose exposure to all factors i∈J has been added or turned on, whereas all confounding factors i∉J remain at the levels specified by x. Let =(X1,…,Xp) refer to the exposure vector of a randomly chosen subject in the population. The joint proportion of risk due to the predictors in J is

when averaging the effect of all confounders. The numerator of eq. (5) is the average causal effect (ACE) when causation can be replaced by association (Agresti 2013), and the denominator normalizes ACE to a value between –1 and 1.

It is also possible to define counterparts of eqs (2), (3) and (4) for a whole population. The population attributable fraction

quantifies the fraction of disease prevalence Eθ(X) that is jointly attributable to a subset J of risk factors. It is an important quantity for public health and a tool for deciding which possible interventions that should be prioritized, see Deubner et al. (1980), Northridge (1995), Benichou (2001) and Sjölander and Vansteelandt (2011). For preventive factors, one uses instead the population prevented fraction

i. e. the fraction by which prevalence is decreased when the population is jointly exposed to all factors in J. Analogously to eq. (4), we combine eqs (6) and (7) into a generalized version

of population attributable fraction that takes values between –1 and 1 and equals PAFtrad for risk factors and −PPF for preventive factors.

It is possible to reinterpret eqs (4), (5) and (8) to quantify epistasis or interaction among a subset J of at least two predictors, rather than their joint effect. Several methods are available to quantify interaction, see Greenland (1983), Cordell (2002), Phillips (2008), VanderWeele (2010) and VanderWeele and Knol (2014) for reviews. Some epistasis definitions involve causation, counterfactuals or potential outcomes (Rothman 1976a), and under certain conditions they can be tested statistically (VanderWeele and Robins 2008; Sjölander et al. 2014; Lekman et al. 2014). We will adopt a common and weaker definition of epistasis in terms of a model whose interaction parameter between a subset J of factors is nonzero. When attributable proportion is used to quantify the strength of this interaction (rather than the joint effect), we must not include any marginal effects. This amounts to reinterpreting θrem(x) as the disease risk when exposure to each factor in J is the same as in the definition of θ(x), but interaction among them, not their marginal effects, is removed. This requires that a model M of no interaction is defined. Statisticians have commonly used logistic regression models and multiplicative odds, but additive models of interaction have recently gained popularity, since they are believed to approximate an underlying biological reality more accurately (Rothman 2012). In particular, a collection of factors exhibits sufficient cause interaction when they are all present in at least one causal mechanism that is sufficient to bring about the disease. For binary factors and outcomes, it has been shown in VanderWeele (2009) and references therein, that under mild regularity conditions, sufficient cause interaction is closely related to additive interaction.

Given that a model M has been specified, interaction is positive or negative when there is excess or lowered risk due to interaction, depending on the sign of APtrad. For positive interaction that increases disease risk (synergism), APtrad is a number between 0 and 1 that quantifies the proportion by which disease risk is lowered when interaction among the factors in J is removed for subjects with exposure x. But for negative interaction that lowers disease risk (antagonism), APtrad has no lower bound. The alternative definition eq. (4) is identical to eq. (2) when positive, but normalized differently for negative values, as minus the proportion by which disease risk is lowered for a subject with exposure x when interaction among the factors in J is added, with a lower bound of –1. The two extreme cases when AP equals –1 and 1 correspond to interaction that completely eliminates or dominates disease risk. The quantities AP and PAF in eqs (5) and (8) are also numbers between –1 and 1, with similar interpretation as the proportion of disease risk when interaction among the predictors in J is removed or added.

As far as we know, APtrad has only been used for interaction between pairs of binary factors under additive models of no interaction, but we will apply our new normalization of AP and PAF much more generally. Our framework is a model with an arbitrary number p of factors, when the joint effect or interaction among a subset J of them is of interest. We consider a large class of models M of no interaction, such as an additive, multiplicative or disjunctive one. Results of Hosmer and Lemeshow (1992), Assman et al. (1996) and Andersson et al. (2005) are extended, whereby linear or logit parametrization of effect parameters is used to obtain point estimates and three types of confidence intervals, two Wald intervals based on the multivariate delta method and asymptotic normality of parameter estimates, and a bootstrap interval. In particular, we utilize that eq. (4) is restricted to [−1,1] and show that a Wald interval based on a logit type transformation and the bootstrap work very well for a wide range of models.

Our work was motivated by a recent study in Lekman et al. (2014), where possible genetic interaction effects of major depressive disorder were sought for in terms of pairs of single nucleotide polymorphisms (SNPs), using different multiplicative and additive models of no interaction. For some pairs of SNPs, AP estimates for the additive models were far below –1, with difficulty of interpretation. Although different normalizations of AP due to interaction have been proposed (VanderWeele 2013), the one in eq. (4) is the first that guarantees a range of [−1,1].

In order to estimate AP and PAF one typically uses cohort data, and external information is needed in order to calculate these quantities from case-control data. But since odds ratios can be estimated consistently from case-control data, an alternative strategy is to redefine the attributable proportion in terms of odds ratios rather than disease risks. We apply such an odds ratio based definition of AP, with our new normalization, to estimate interaction effects of a multiple sclerosis (MS) data set previously analyzed in Hedström et al. (2011).

The paper is organized as follows. In Section 2 we consider the case of two factors and additive models of no interaction, and then in Section 3 an arbitrary number of factors, and a generalized linear model (McCullagh and Nelder 1989) family of no interaction models that include additive, additive odds, multiplicative and disjunctive models. In Section 4 we define the analogous versions of AP in terms of odds ratios. Then we define point estimates and confidence intervals for prospective and retrospective studies in Section 5, conduct a simulation study in Section 6, analyze the MS data set in Section 7, and conclude with a discussion in Section 8. Further mathematical and implementation details are given in Appendices A-F.

2 Two predictors

3 Arbitrary number of predictors

4 Odds ratios

It is possible to estimate relative risks consistently for cohort studies or case-control studies with incident cases who developed the disease after exposure (Miettinen 1976). But many case-control studies involve prevalent cases that were diagnosed before the study began. It is not possible then to estimate relative risks consistently, unless the sampling fraction of cases and controls is known. But since odds ratios

can be estimated consistently also for the latter kind of studies (Cornfield 1951, Prentice and Pyke 1979, Skrondal 2003) they may replace relative risks, both for AP due to marginal effects eq. (34) and interaction eq. (37).

Consider a subset J⊂{1,…,p} of factors, and a non-interaction model M among them. In Appendix A.3 we define a number of odds-ratio based quantities in terms of eq. (19) and

the odds ratio for exposure x when interaction among the factors in J has been removed, with θrem(x; J,M) as in eq. (17). This includes the joint attributable proportion APOR(x; J) of the odds ratio for exposure x due to the factors in J, the corresponding joint attributable proportion APOR(J) when averaging the effect of all confounders, the attributable proportion APOR(x; J,M) of the odds ratio due to interaction among the factors in J for exposure x, and finally, the corresponding attributable proportion APOR(J,M) of the odds ratio due interaction among the factors in J, when averaging over the confounding factors. When APOR>0, it is the relative proportion of the odds ratio of disease attributed to joint effect or interaction among the variables in J. On the other hand, when APOR<0, it is minus the proportion by which joint effect of interaction among the variables in J decreases the odds of developing the disease.

Odds ratios typically approximate relative risks well for rare diseases. When this is the case, the odds ratio based versions of the attributable proportion are close to the corresponding risk based quantities of Section 3. For the models studied by Kalilani and Atashili (2006), the odds ratio approximation was accurate for diseases with baseline risks up to 0.2. But sometimes the difference is notable, and depending on disease risks and mode of interaction the odds ratio based attributable proportion may be biased (Kalilani and Atashili 2006; Zou 2008).

The attributable proportion of the odds ratio due to interaction has a particularly explicit form for the additive odds model. For instance, when there are p=2 risk factors and x=(1,1), we tacitly assume in the notation that J={1,2} (since otherwise APOR=0) and find that

is an odds ratio analogue of eqs (13)–(14). The corresponding traditional definition

can be found for instance in Rothman (2012).

5 Inference

We will outline a unified approach to estimate the attributable proportion eqs (4)–(5) or eqs (40)–(41), and the population attributable fraction [8] due to joint effects or interaction of an arbitrary subset J⊂{1,…,p} of factors. It is assumed that a saturated model Msat is fit to data. For the risk-based quantities eqs (4)–(5) or eq. (8), the parameter vector ψ=(ψx) of Msat has 2p components, one for each x. A linear parametrization of the risk function is

Another possibility is to parametrize a logit transformation of risk. This was utilized by Hosmer and Lemeshow (1992) for the traditional definition [12] of AP of interaction with p=2 risk factors and an additive model of no interaction. They noted that the disease probability for a subject with exposure x can be parametrized in two ways. Either

with y=(y1,…,yp) and y≤x interpreted componentwise, i. e. yi≤xi for i=1,…,p, and ψx an interaction parameter defined for the |x|=∑ixi nonzero risk factors of x. This is an ordinary logistic regression parametrization for a model with p predictors having 0 as their baseline levels, and interactions of all orders (Agresti 2013). If all nonzero exposure profiles x are viewed as different levels of one single predictor variable, then

with ψx the log odds of x compared to the baseline level 0. Since the two parametrizations eqs (23) and (24) only differ by a linear transformation, they give identical inference, although eq. (24) is often simpler. We can also express the odds ratio eq. (19) conveniently in terms of any of the two logistic parametrizations, either

for eq. (23), or

for eq. (24). Neither eq. (25) nor eq. (26) involve the baseline parameter ψ0, and therefore we redefine ψ to equal (ψx; x≠0) when quantities with odds ratios are estimated.

All quantities AP(x),APOR(x),AP,APOR and PAF of interest are functions

of the parameter vector ψ, with different choices of a and b, as summarized in Table 5. In order to estimate ξ, suppose we have data (xi,Yi), i=1,…,N from N subjects, with ψˆ=(ψˆx) the accompanying maximum likelihood estimator of ψˆ. If this is a cohort design or some other prospective study, ψˆ maximizes a prospective likelihood, the probability of disease given exposure. For the linear parametrization (22) we have that

with Nxy the number of subjects i with xi=x and Yi=y, and Nx=Nx0+Nx1. For either of the two logit scale parametrizations (23)–(24), we use standard logistic regression software to compute ψˆ. This is also possible for a case-control study, where sampling is on the response variable, and ψˆ maximizes a retrospective likelihood, the probability of exposure given disease, as long as we focus on quantities that only involve odds ratios (Prentice and Pyke 1979). We estimate ξ by replacing ψ with ψˆ in eq. (27), so that

For large N one often approximates the distribution of ψˆ by a multivariate normal N(ψ,Σ) with asymptotic covariance matrix Σ. The denominator of eq. (27) is not differentiable when a(ψ)=b(ψ). But since the numerator is zero at such points, ξ(ψ) is still differentiable and ξˆ approximately normally distributed Nξ,σ2 with asymptotic variance σ2=DΣDT,D=∂ξ(ψ)/∂ψ, and DT its transpose. As in Hosmer and Lemeshow (1992) we use the delta method to find a Wald type confidence interval

for ξ with asymptotic coverage probability 1−α, where zβ is the β-quantile of a standard normal distribution, and

the standard error of ξˆ, with Σˆ and Dˆ estimates of Σ and D. For the linear parametrization (22), Σˆ is a diagonal matrix with elements

since observed risks are independent binomial proportions. For the logit parametrizations, Σˆ is available from standard logistic regression software. A formula for Dˆ is derived in Appendix B.

When ξ is close to one of its two boundary values –1 and 1, the distribution of ξˆ will typically be skewed, unless the sample size N is very large. Since eq. (22) assumes a symmetric normal distribution of ξˆ, the actual coverage probability π(1−α)=P(ξ∈I) may deviate substantially from the nominal 1−α. This suggests that a conveniently chosen transformation of h(ξˆ) of ξˆ will have a distribution closer to normal. When ξ∉{−1,1}, the logit type transformation h(x)=log(1+x)/(1−x) gives a Wald confidence interval

where σ˜=h′(ξˆ)σˆ=2σ˜/{(1+ξˆ)(1−ξˆ)}. Whereas the interval in eq. (29) may extend outside (−1,1), eq. (32) is always a subset of (−1,1) as long as ξˆ∈(−1,1). Other confidence interval approaches for ξ that deal with asymmetric distributions include the methods of variance estimates recovery (Zou 2008), likelihood ratio based intervals obtained from a saturated additive odds ratio model (Richardson and Kaufman 2009), nonparametric bootstrap (Assman et al. 1996), with accelerated bias correction (Zou 2008). In Appendix D we define a bias-corrected and accelerated version of the percentile bootstrap method.

The exposure distribution q(x)=P(X=x) is assumed to be known in eq. (28) for the population based AP, PAF and APOR. In Appendix C we show how to compute standard error and confidence intervals when q is estimated from a prospective study.

6 Numerical illustration and simulation

Figures 1–2 illustrate how normalizations eqs (2) and (4) affect AP(1,1; M) of interaction for two saturated models I and II with p=2 factors. Model I has two risk factors, whereas II has one risk and one preventive factor. As expected, the normalization has a large impact for both models when AP is negative. The model of no interaction influences the value of AP a lot for model I, and to some extent for model II. Table 2 gives values of the disease risk θrem(x; J,M) when various types of interaction among the variables in J is removed. This is done for models I and II, as well as for a model III with p=3 co-factors.

Figure 1:

Plots of normalized attributable proportion AP(1,1)=AP(1,1; J,M) of interaction [37] (solid) and the corresponding traditional quantity [2] (dashed). The curves are functions of θ(1,1) for a saturated model with p=2 co-factors, whose other three risks θ(0,0)=0.05, θ(1,0)=0.25, θ(0,1)=0.4 are the same as for Model I of Table 2. Both variables are non-confounders (J={1,2}), and three different no interaction models M are assumed; additive (thick line), additive odds (upper thin line) and multiplicative (lower thin line).

Figure 2:

Plots of normalized attributable proportion AP(1,1)=AP(1,1; J,M) of interaction [37] (solid) and the corresponding traditional quantity [2] (dashed) versus θ(1,1) for a saturated model with p=2 risk factors, whose other three risks θ(0,0)=0.1, θ(1,0)=0.05, θ(0,1)=0.15 are the same as for Model II of Table 2. Both variables are non-confounders (J={1,2}), and two different no interaction models M are assumed; additive (thick line), and multiplicative (thin line). The curves for the additive odds model (not shown) are very close to those of the additive model.

Table 3 displays estimated coverage probabilities of confidence intervals for the odds ratio based attributable fraction APOR(x; J,M) of interaction. Each coverage probability is estimated from a number of randomly drawn data sets with N/2 cases and controls, generated from the two saturated models I and II with p=2 co-factors, and the model III with p=3. We use two models M of no interaction, additive odds and multiplicative, and three types of confidence intervals; the delta method (29), the logit delta method (32) and the bootstrap method (48). It is seen that the accuracy of all three confidence intervals increase with N, and overall the bootstrap method has the best performance, although the simpler logit delta method works well for a wide range of models. Appendix F provides more details about the simulation study.

7 A real data set

MS is a complex and inflammatory disease causing damage to the central nervous system. A number of family studies indicate a genetic component of the disease, with a major contribution from variants at the human leukocyte antigen (HLA) complex. It is well known that presence of allele 15 of the HLA-DRB1 gene is a risk factor, whereas allele 02 of the HLA-A gene has a protective effect (Cree 2014). Several data sets reveal that environmental factors, in particular smoking, impact the aetiology of the disease as well. A recent paper of Hedström et al. (2011) demonstrated pairwise positive additive interaction between these three factors for a Scandinavian case-control data set with 843 cases and 1209 controls. Here we will redo parts of the analysis, using our novel measures APOR of joint effects and interaction. As in Hedström et al. (2011), we use a binary coding xDRB1∗15=1 for an individual with genotype 15/X or 15/15, xA∗02=1 if the genotype is Z/Z, and xsmoking=1 for a current smoker, with X and Z referring to any allele different from HLA-DRB1*15 and HLA-A*02 respectively. For simplicity, we will not adjust for the confounding effect of age, gender, residential area and ancestry, as in Hedström et al. (2011). Table 4 summarizes results for different models, three of which have one single factor, one with both genetic factors included, and one model with all three factors present.