Modern Stochastics: Theory and Applications

Many complex diseases are influenced by a number of risk factors that interact in a complicated way. This is often quantified by means of a regression model with affection status of a given disease as binary response, whereas the risk factors and possibly some other variables are chosen as covariates. Logistic regression models have often been used to quantify main effects and strength of interaction among the risk factors with regards to disease. There are several reasons for this. The logistic transformation is first of all the canonical link of a generalized linear model with a binomially distributed response, and the parameters of this model have a straightforward multiplicative odds interpretation [20]. A second reason is that many epidemiological datasets are collected retrospectively based on outcomes rather than on covariates. In particular, it is well known that many parameters of the logistic regression model can be estimated consistently for case-control studies under suitable sampling assumptions on the cases and controls [24]. There are also additive models of joint effects and interaction. They have recently gained in popularity, since they are believed to approximate a biological system with causal mechanisms more accurately than multiplicative odds [30, 26]. The additive measures of interaction are functions of relative risks between individuals with or without exposure to the risk factors, such as the relative excess of risk due to interaction (RERI), the attributable proportion (AP) due to interaction or the synergy index (SI) [12, 28]. Although these relative risks cannot be estimated consistently for case-control data, they are well approximated by estimable odds ratios when the disease risk is small [15, 4, 3]. When relative risks are replaced by odds ratios in the expressions for RERI, AP and SI, they correspond to additive odds models of interaction.

In previous work [16] we developed a unified theory for quantifying and estimating different measures of marginal effects, joint effects and interaction among risk factors on a multiplicative, additive, additive odds or some other scale. We also described how to estimate and produce confidence intervals for these quantities from a prospective study, with data sampled based on their covariates, or from a case-control study. Traditional definitions of attributable proportion have an unbounded negative range [18, 29, 19]. In [16], we introduced a novel normalization of AP that guarantees a range between $-1$ and 1, with negative or positive values depending on whether there is synergism or antagonism between the risk factors.

In this paper we concentrate on case-control data and additive odds measures of main effects, joint effects and interaction. Our main result is to express the odds ratio of the risk factors as a sum of terms, which include their main effects and different orders of interaction, when the effect of other confounding covariates is controlled for. In this way we extend and unify some previously used measures of interaction [23, 17, 16]. In order to find confidence intervals (CIs) for the attributable proportion, synergy index and excess odds ratio due to joint effects and interaction from case-control data, we first estimate the parameters of a logistic regression model by maximum likelihood. Then we use the delta method [27, 6] for an appropriate transformation of these measures of joint effects and interaction in order to find their standard errors and CIs on the new scale, before transforming back to the original odds scale.

The paper is organized as follows: In Section 2 we define a logistic regression model that includes marginal and interaction effects for the risk factors of interest, as well as marginal effects of other confounding covariates. The additive odds measures of joint effects and interaction are defined in Section 3, and the procedures for estimating them are given in Section 4. Our methodology is illustrated in Section 5 for a multiple sclerosis (MS) dataset from Hedström et al. [14], and a concluding discussion appears in Section 6.

Let $Y\in \{0,1\}$ be a binary outcome variable, with $Y=1$ for individuals that carry a certain disease, and $Y=0$ for those that do not. Consider a large population, and let be the disease probability of a randomly chosen individual. It is assumed to be a function of $p+q$ covariates $\boldsymbol{x}=(\boldsymbol{v},\boldsymbol{z})=(v_{1},\dots ,v_{p},z_{1},\dots ,z_{q})$, of which the first p are binary risk factors $v_{1},\dots ,v_{p}$ of main interest, with $v_{j}=1$ indicating presence and $v_{j}=0$ absence of each such factor j. The other q covariates $z_{1},\dots ,z_{q}$ are not necessarily binary, and they are included in the model as possible confounders. We will parametrize the disease probability on a logit scale with $\boldsymbol{w}=(w_{1},\dots ,w_{p})$ a binary vector and $\mathbf{0}=(0,\dots ,0)$ a vector with zero components. Inequalities between vectors are interpreted componentwise, so that the first sum in (2) is taken over all nonzero vectors $\boldsymbol{w}$ such that $w_{j}\le v_{j}$, with $w_{j}=1$ for at least one factor j. This sum equals 0 when $\boldsymbol{v}=\mathbf{0}$.

The logistic model in (2) is saturated for the binary risk factors, including all orders of their interaction, whereas it is linear in the other covariates [1, 2]. In particular, if $v_{j}=1$ and all other components of $\boldsymbol{v}$ are zero, $\psi _{\boldsymbol{v}}$ quantifies the marginal effect of factor j in absence of the others, whereas $\psi _{\boldsymbol{v}}$ for $\boldsymbol{v}$ such that $|\boldsymbol{v}|=\sum _{j}v_{j}\ge 2$ quantifies interaction of order $|\boldsymbol{v}|$ among those factors j for which $v_{j}=1$.

It is assumed that $\boldsymbol{\kappa }=(\kappa _{0},\kappa _{1},\dots ,\kappa _{q})$ are nuisance parameters, whereas are the ${2}^{p}-1$ structural parameters of main interest that quantify marginal effects and interaction among the p risk factors, with $\mathbf{1}=(1,\dots ,1)$. Since they are defined on a logit scale, it is possible to estimate them from suitable case-control datasets. In the next section we will develop alternative measures of marginal effects, joint effects and interaction that are functions of $\boldsymbol{\psi }$, but expressed in terms of odds ratios.

Our purpose is to estimate and produce confidence intervals for parameters that quantify marginal effects, joint effects or interaction between a subset $J\subset \{1,\dots ,p\}$ of factors, when the levels of the remaining factors in $K=\{1,\dots ,p\}\setminus J$, and the q other confounder variables are controlled for. It is assumed that (4) only involves the structural parameters $\boldsymbol{\psi }$, not the nuisance parameters $\boldsymbol{\kappa }$. The reason is that we consider parameters (4) that are functions of odds ratios i.e. ratios between the disease odds of two subjects with identical confounding variables $\boldsymbol{z}$ but different risk exposures $\boldsymbol{v}$ and $\mathbf{0}$ respectively, for some or all of the ${2}^{p}-1$ vectors $\mathbf{0}<\boldsymbol{v}\le \mathbf{1}$. In particular, when the disease risks in (1) are small, the odds ratio in (5) approximates the relative risk well.

After possible reordering of factors we may assume without loss of generality that those in J come first, so that $\boldsymbol{v}=(\boldsymbol{v}_{J},\boldsymbol{v}_{K})$, with $\boldsymbol{v}_{J}=(v_{j};\hspace{0.1667em}j\in J)$ and $\boldsymbol{v}_{K}=(v_{j};\hspace{0.1667em}j\in K)$ being the exposure levels of the factors in J and K respectively. We will consider linear combinations of the odds ratios (5), and the following concept will be central for the rest of the paper:

When $|\boldsymbol{v}_{J}|=0$, (7) is the odds ratio when none of the risk variables in J are turned on, and when $|\boldsymbol{v}_{J}|=1$ it is the marginal odds ratio increment of one single factor. For $|\boldsymbol{v}_{J}|\ge 2$ we interpret (7) as a measure of additive odds interaction among those factors $\{j\in J;\hspace{0.1667em}v_{j}=1\}$ in J that are at risk, when all interaction terms of lower order among these factors have been removed. More explicitly, the terms of order 0, 1, 2 and 3 in (7) have the form where $\boldsymbol{e}_{j}$, $\boldsymbol{e}_{jk}$ and $\boldsymbol{e}_{jkl}$ are binary vectors of length $|J|$ with zero components, except in positions $\{j\}$, $\{j,k\}$ and $\{j,k,l\}$ respectively.

It is possible to give another interpretation of (7). Let $\boldsymbol{w}=(w_{j};\hspace{0.1667em}j\in J)$ be a vector with $0\le w_{j}\le 1$, and extend the domain of the first $|J|$ components of the odds ratio $\text{OR}_{(\boldsymbol{w},\boldsymbol{v}_{K})}$ to ${[0,1]}^{|J|}$. If $\boldsymbol{v}_{J}$ is a binary vector of length $|\boldsymbol{v}_{J}|=l>0$ with nonzero components $v_{j_{1}}=\cdots =v_{j_{l}}=1$, then ${\varDelta }^{(\boldsymbol{v}_{J},\boldsymbol{v}_{K})}\text{OR}$ is a finite difference approximation of ${\partial }^{l}\text{OR}_{(\boldsymbol{w},\boldsymbol{v}_{K})}/(\partial w_{j_{1}}\dots \partial w_{j_{l}})$, where $\boldsymbol{w}=\boldsymbol{v}_{J}/2$.

Having defined the odds ratio increment in (7), we are ready to state the main result. It tells that the odds ratio for exposure $\boldsymbol{v}_{J}$ is a sum of the odds ratio increments. This includes the baseline odds ratio when none of the factors in J are turned on, the additive marginal odds ratio increments for those $|\boldsymbol{v}_{J}|$ factors in J that are exposed, and the odds ratio interactions of order $2,3,\dots ,|\boldsymbol{v}_{J}|$ among these factors (see Appendix A for a proof):

The expansions in (7) and (9) define a combinatorial inclusion–exclusion principle. In order to see this, we associate to each $j\in J$ a set $A_{j}$ and its complement ${A_{j}^{c}}$. Then ${\varDelta }^{(\boldsymbol{v}_{J},\boldsymbol{v}_{K})}\text{OR}$ represents an intersection $(\bigcap _{j;v_{j}=1}A_{j})\cap (\bigcap _{j;v_{j}=0}{A_{j}^{c}})$, whereas $\text{OR}_{(\boldsymbol{v}_{J},\boldsymbol{v}_{K})}$ corresponds to a union $\bigcup _{j;v_{j}=1}A_{j}$, and in particular $\text{OR}_{(\mathbf{0},\boldsymbol{v}_{K})}=\varnothing $.

It is of interest to know how much of the odds ratio (5) can be explained by marginal effects and lower order interaction among the factors in J. For this reason we introduce the following concept:

It is possible to rewrite the predicted odds ratio (10) as a linear combination of lower order odds ratios $\{\text{OR}_{(\boldsymbol{w},\boldsymbol{v}_{K})};\hspace{0.1667em}0\le \boldsymbol{w}\le \boldsymbol{v}_{J},|\boldsymbol{w}|\le i\}$ for exposure vectors $\boldsymbol{w}$ that have at most i factors in J at risk (see Appendix B for a proof):

Assume from now on that $\boldsymbol{v}_{J}=\mathbf{1}$, so that all factors in J are at risk. We will quantify how much of the odds ratio $\text{OR}_{(\mathbf{1},\boldsymbol{v}_{K})}$ is left unexplained by marginal effects and lower orders of interaction among the factors in J. To this end, we introduce the following concept:

The unadjusted excess odds ratio can be interpreted as an unstandardized residual of a regression model, where only marginal effects and interaction terms up to order $i-1$ are included as independent variables in order to predict the odds ratio (5). Equivalently, it quantifies the contribution of the odds ratio expansion (9) from terms of order at least i. The special case $i=|J|$ is treated by Katsoulis and Bamia [17] when $K=\varnothing $. The definition of the unadjusted odds ratio then simplifies to ${\varDelta }^{(\mathbf{1},\boldsymbol{v}_{K})}\text{OR}$.

We will define three measures (4) of marginal, joint or interaction effects among the factors in J, and all of them are functions of $\text{UnadjEOR}$.

All three quantities in equations (14)–(16) are stratified for the risk factors in K, like confounders that are controlled at their observed levels. In contrast, there is only partial (additive) control for the remaining q covariates of $\boldsymbol{z}$. EOR and AP can be viewed as different types of standardized residuals of a regression model where marginal effects and lower order interactions among the factors in J are used to predict the odds ratio (5). The inverse synergy index ${\text{SI}}^{-1}$ is the analogue of the coefficient of determination. It quantifies how large fraction of the odds ratio increment above the baseline level $\text{OR}_{(\mathbf{0},\boldsymbol{v}_{K})}$ that is explained by the regression model.

Some special cases of formulas (14)–(16) are of particular interest. When $i=|J|=1$, there is one single risk factor ($J=\{1\}$). Equation (14) then quantifies the excess odds ratio due the marginal effect (EORM) of factor 1, when those in K are controlled for at level $\boldsymbol{v}_{K}$. Similarly, equation (15) defines the attributable proportion of risk due to the marginal effect of factor 1, whereas the synergy index is not well defined.

When $i=1$ and $|J|\ge 2$, we refer to the quantity in (14) as EORJ, the excess odds ratio due to joint (marginal and interaction) effects of all factors in J when those in K are controlled at level $\boldsymbol{v}_{K}$. In the same way, $\text{AP}$ is the attributable proportion due to joint effects among the factors in J, whereas the synergy index is undefined.

When $i=2$ we refer to the quantity in (14) as the total excess odds ratio due to all levels of interaction (TotEORI) among the factors in J, when those in K are controlled at level $\boldsymbol{v}_{K}$. $\text{AP}$ and $\text{SI}$ are similarly referred to as the attributable proportion and synergy index due to all levels of interaction among the factors in J.

Finally, when $i=|J|$ we refer to the quantity in (14) as the excess odds ratio due to the highest order $|J|$ of interaction (EORI) among the factors in J, when those in K are controlled at level $\boldsymbol{v}_{K}$. Analogously, equations (15)–(16) quantify the attributable proportion and the synergy index due to the highest order $|J|$ of interaction among the factors in J.

Since EOR, AP and SI are all functions of UnadjEOR, it suffices to specify the latter. In the subsections to follow we will do so for models with 1, 2, or 3 risk factors in J.

Assume that a case-control dataset $(\boldsymbol{x}_{1},Y_{1}),\dots ,(\boldsymbol{x}_{n},Y_{n})$ of size $n=n_{0}+n_{1}$ is available, with $n_{0}$ controls ($Y_{a}=0$) and $n_{1}$ cases ($Y_{a}=1$). Since this is a retrospective sample, we cannot estimate all parameters of the logistic regression model (2). On the other hand, it is possible to estimate the structural parameters $\boldsymbol{\psi }$ consistently when controls are frequency matched with cases. This is the unconditional logistic regression approach, whereby the prospective log likelihood is maximized over the model parameters in (2). This approach can also be used if cases and controls are matched in strata, as long as all variables used for matching are included in the model as covariates, and the number of strata does not grow with sample size. On the other hand, conditional logistic regression is more appropriate if the number of strata is large [5].

Let be the plug-in estimator of (4) obtained from the maximum likelihood estimator (17). In order to produce a confidence interval for ξ with asymptotic coverage probability $1-\alpha $ we use the delta method in conjunction with some appropriately chosen monotone increasing and differentiable transformation h. This leads to where $z_{\beta }$ is the β-quantile of a standard normal distribution and $\text{SE}$ is the standard error of $h(\hat{\xi })$. This method relies on asymptotic normality of the structural part of (17) for large samples, with $\boldsymbol{\varSigma }=\boldsymbol{I}{(\boldsymbol{\psi },\boldsymbol{\kappa })}^{-1}$ an asymptotic approximation of the covariance matrix $\text{Cov}(\hat{\boldsymbol{\psi }})$ of $\hat{\boldsymbol{\psi }}$, which equals the inverse of the Fisher information matrix $\boldsymbol{I}$ [6]. In order to find $\text{SE}$ we first approximate the asymptotic variance of $\hat{\xi }$ by ${\sigma }^{2}=\boldsymbol{D}\boldsymbol{\varSigma }{\boldsymbol{D}}^{T}$, using (20) and the first order Taylor expansion of $\xi (\boldsymbol{\psi })$, with $\boldsymbol{D}=\boldsymbol{D}(\boldsymbol{\psi })=d\xi (\boldsymbol{\psi })/d\boldsymbol{\psi }$ and T referring to vector transposition (see Appendix C for an explicit expression of $\boldsymbol{D}$). Then one chooses the function h so that the distribution of $h(\hat{\xi })$ is closer to normal than that of $\hat{\xi }$. Typically h is a variance stabilizing transformation that maps the range of ξ to the real line $(-\infty ,\infty )$. The Taylor expansion of h gives the asymptotic variance of $h(\hat{\xi })$. By taking the square root of an estimate of this variance, we find the standard error of $h(\hat{\xi })$, with ${\hat{\sigma }}^{2}=\hat{\boldsymbol{D}}\hat{\boldsymbol{\varSigma }}{\hat{\boldsymbol{D}}}^{T}$ an estimate of ${\sigma }^{2}$, $\hat{\boldsymbol{D}}=\boldsymbol{D}(\hat{\boldsymbol{\psi }})$ an estimate of $\boldsymbol{D}$, and $\hat{\boldsymbol{\varSigma }}=\boldsymbol{I}{(\hat{\boldsymbol{\psi }},\hat{\boldsymbol{\kappa }})}^{-1}$ the inverse of the observed Fisher information matrix. In this paper we will use the transformations see for instance Rothman [25] and Hössjer et al. [16]. The three functions in (21) map the corresponding indices ξ from their original ranges ($(-\infty ,\infty )$ for EOR, $(-1,1)$ for AP, $(0,\infty )$ for SI) to $(-\infty ,\infty )$.

Multiple sclerosis (MS) is a complex and inflammatory disease causing damage to the central nervous system. Its prevalence is over 0.1% in many countries, affecting large regions of the world [21]. There is solid evidence for a genetic component of the disorder, with a major contribution from variants at the human leukocyte antigen (HLA) complex. It is also 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 [7]. Several studies reveal that environmental factors, in particular smoking, impact the risk of the disease as well [22]. A Swedish study in Hedström et al. [13] demonstrated positive pairwise additive interaction between the two genetic factors, and also between smoking and each genetic factor. These results have more recently been replicated and refined in Hedström et al. [14], by merging case-control studies from several countries. Due to the size of this meta analysis, it was also possible to investigate whether third order interaction was present between the two genetic factors and smoking. In order to illustrate the methodology of this paper, we present some of the findings from the four Nordic studies (see Table 1) of Hedström et al. [14].

Apart from the two genetic factors and smoking, three other covariates (gender, age, study) were also part of the model. This gives a total of 8 covariates, encoded as The last three study covariates were only included for the meta analysis, so that Let be the odds ratio for the joint effect of all the risk factors J, when the confounding risk factors in K are fixed at level $\boldsymbol{v}_{K}$, and the remaining covariates are $\{1,\dots ,p+q\}\setminus (J\cup K)$. Table 2 lists maximum likelihood estimates (18) of these odds ratios for various choices of risk factors and confounders. The associated confidence intervals (19) are obtained using the logarithmic transformation in (21). Although MS is a common disorder, its prevalence is small enough to assume that the point estimates and confidence intervals of Table 2 are representative for the relative risks $\text{RR}_{(\mathbf{1},\boldsymbol{v}_{K},\boldsymbol{z})}/\text{RR}_{(\mathbf{0},\boldsymbol{v}_{K},\boldsymbol{z})}$ as well.

We find, for instance, that the point estimate of the marginal odds ratio (or relative risk) of having MS in the combined dataset is 3.6 for individuals with the DRB15 risk allele, compared to those that lack this allele. The corresponding marginal odds ratios for absence of the protecting A2 allele and for smoking are 1.75 and 2.0 respectively. Since the joint odds ratios for all pairs of risk factors are much larger than the corresponding marginal odds ratios, there are strong indications of two-way interactions between all pairs of risk factors. There is possibly some three-way interaction between DR15, A2- and smoking as well, since the joint OR for all three factors is higher than the pairwise odds ratios. On the other hand, the OR for the two genetic factors is only higher among smokers than among non-smokers for one study (EIMS).

The estimates of Table 2 motivate further analysis of the MS datasets in Hedström et al [14] based on the attributable proportion, excess odds ratio and synergy index for the three risk factors DR15, A2- and smoking. Table 3 gives confidence intervals for all three quantities when $J=\{1,2,3\}$ and $K=\varnothing $, for various choices of i. It confirms a strong joint effect of all three factors, since AP and EOR are both significantly different from 0, and SI is significantly different from 1. For instance, the estimate $\widehat{\text{AP}}=0.92$ for the combined Nordic data set indicates that about 92% of the odds ratio (or disease risk) for smokers with both genetic risk factors, is not present among those that lack all three factors. We also find that the total amount of second and third order interaction between DR15, A2- and smoking is strongly significant. In particular, $\widehat{\text{AP}}=0.52$ for the Nordic metapopulation indicates that about half of the disease risk is due to interaction. One the other hand, $\widehat{\text{SI}}=2.38$ tells that the disease risk increment of a smoker with both genetic risk factors, compared to an individual with none of these risk factors, is more than twice the disease risk increment due to marginal effects only. Finally, we find from the rightmost AP, EOR and SI columns of Table 3 that third order interaction between DR15, A2- and smoking is only significant for one dataset (EIMS).

We have also estimated the attributable proportion separately for males and females (data not shown). The results are in quite good agreement with the upper part of Table 3, although the values for males are slightly larger than those for females, and since the gender specific datasets are smaller, the confidence intervals are wider. Such an analysis illustrates how joint effects of and interaction among three factors ($J=\{1,2,3\}$) is affected when one controls for a fourth factor ($K=\{4\}$).

In this paper we studied how a collection J of binary risk factors jointly influence a binary outcome. Our objective was to develop a general framework for quantifying marginal effects and various orders of interaction between these factors on an additive scale, when stratifying or partially controlling for other confounding variables. This led to the key result; an expansion of the odds ratios as a sum of marginal effects and interaction terms of different orders, with a finite difference and a combinatorial inclusion–exclusion interpretation. We also showed how to use these odds ratio expansions for estimating and producing confidence intervals for the excess odds ratio, attributable proportion and synergy index from a case-control dataset. The methodology was illustrated using a Nordic meta study for multiple sclerosis. The inferential procedure relies on maximum likelihood estimates from a logistic regression model and the delta method. It has been implemented in a Matlab program that is available from the first author upon request.

Our approach makes it possible to stratify for some variables (in K) at various levels, and yet use the whole data set to estimate effect parameters of the other q confounders, for which there is only partial control. But this requires a large data set in order to estimate all ${2}^{p}+q={2}^{|J|+|K|}+q$ parameters of the logistic regression model. An alternative and simpler strategy is to use only those observations $(\boldsymbol{x}_{a},Y_{a})$ for which the variables in K are at a pre-specified level. This gives a smaller model with only ${2}^{|J|}+q$ parameters to estimate.

The delta based confidence intervals are fast to compute. This is suitable in applications where a number of different putative risk factors are sought for. We have implemented confidence intervals based on resampling as well, using the bias-corrected accelerated percentile method [9, 10, 16]. There is generally good agreement between the resampling and delta-based intervals for odds ratios and attributable proportions, as long as the interaction effects are not too strong, the order of interaction is not too high and the data set is not too small. On the other hand, the delta based confidence intervals tend to be less accurate for excess odds ratios. For this reason we recommend to report excess odds ratios with resampling based confidence intervals. As a topic of future research, it would be of interest to compare the asymptotic accuracy of the delta and resampling based confidence intervals more systematically.

Another object of further study is to develop odds ratio expansions of main effects and interactions when some of the risk factors are continuous [11] and to find analogous expansions for non-binary outcomes.