Bayesian estimation of the measurement of interactions in epidemiological studies

Measures of additive interaction

Consider two binary factors, A and B, with values of 0 and 1, where 0 denotes the absence of a factor and 1 denotes its presence. Four possible combinations exist: A₀B₀, A₁B₀, A₀B₁ and A₁B₁. Suppose the reference group is A₀B₀, and let RR₁₀, RR₀₁, and RR₁₁ denote the relative risks for groups A₁B₀, A₀B₁ and A₁B₁, respectively. Rothman, Greenland & Lash (2008) developed three measures of interaction on an additive scale: the relative excess risk due to interaction (RERI), ${\displaystyle RERI=R{R}_{11}-R{R}_{10}-R{R}_{01}+1,}$ the attributable proportion due to interaction (AP), ${\displaystyle AP=\frac{RERI}{R{R}_{11}}=\frac{R{R}_{11}-R{R}_{10}-R{R}_{01}+1}{R{R}_{11}},}$ and the synergy index (S) ${\displaystyle S=\frac{R{R}_{11}-1}{\left(R{R}_{10}-1\right)+\left(R{R}_{01}-1\right)}=\frac{R{R}_{11}-1}{R{R}_{10}+R{R}_{01}-2}.}$ Note that both the RERI and AP can range from - ∞ to + ∞, while S can range from 0 to + ∞. If no interaction is present on the additive scale, both the RERI and AP will be equal to 0, and S will be equal to 1; if there is more than additivity, both the RERI and AP will be greater than 0, and S will be greater than 1; if there is less than additivity, both the RERI and AP will be less than 0, and S will be less than 1.

The RR is readily available in cohort studies. However, in case-control studies, the RR can be approximated by the odds ratio (OR) if the prevalence is sufficiently rare (Appendix S1B). As suggested by Rothman, Greenland & Lash (2008), the RERI, AP and S can be estimated via multiple logistic regression with indicator variables created for categories A₁B₀, A₀B₁ and A ₁B₁. Given that (1)${\displaystyle logit\left(p\right)=log\left(\frac{p}{1-p}\right)={\beta }_0+{\beta }_{1}I\left(A_1{B}_0\right)+{\beta }_{2}I\left(A_0{B}_1\right)+{\beta }_{3}I\left(A_1{B}_1\right)+X\gamma }$ where p = P(case|A, B, X) is the corresponding probability of the case given factors A, B and X (X is a vector of the potential confounders) and I(⋅) is the indicator function; then, OR₁₀ = e^β₁, OR₀₁ = e^β₂ and OR₁₁ = e^β₃ can be substituted for the RRs in the three measures of interaction. As noted by Hosmer Jr, Lemeshow & Sturdivant (2013), the ORs can be represented as OR₁₀ = e^η₁, OR₀₁ = e^η₂ and OR₁₁ = e^{η₁+η₂+η₃} based on logistic regression with factors A and B and the product of A and B (2)${\displaystyle logit\left(p\right)={\eta }_0+{\eta }_{1}A+{\eta }_{2}B+{\eta }_{3}AB+X\gamma .}$

Because both models are saturated, the two estimations are equivalent.

Methods for calculating confidence intervals

The methods used to calculate confidence intervals fall into two categories: symmetric intervals based on a normal distribution and asymmetric intervals based on quantile estimation. Let lnS denote S on the log scale and Z denote any one of the three measures RERI, AP or lnS. As proposed by Hosmer and Lemeshow, the confidence interval of Z is obtained by assuming a normal distribution for Z. The variance in Z can be estimated using the delta method, which employs the first-order approximation of a Taylor series (Hosmer & Lemeshow, 1992). Therefore, the 95% confidence interval of the delta method for Z, which is symmetric about the point estimate, is given by $\hat{Z}\pm 1.96\times \sqrt{\hat{D}\left(Z\right)}$, where $\hat{Z}$ and $\hat{D}\left(Z\right)$ are the point estimate and the variance estimate for Z, respectively. The confidence intervals for the RERI and AP can be obtained directly, while the confidence interval for S can be obtained using an exponential function. Zou (2008) suggested a recover variance estimate of the variance for measures that considers Taylor series expansion using the multivariate delta method.

Another way to estimate confidence intervals is through bootstrapping (Efron & Tibshirani, 1994). This technique allows estimation of the sampling distribution of almost any statistic using random sampling, as do the measures of additive interaction. Let Z denote any one of the three measures RERI, AP or S. The bootstrap samples are resampled from the original sample, and then Z can be estimated in each of the bootstrap samples. When obtaining the distribution of Z, deriving estimates of standard errors and confidence intervals is straightforward. There are several methods for constructing confidence intervals from the bootstrap distribution.

The first of these methods is the normal bootstrap method. The method assumes that Z follows a normal distribution. The variance $\hat{D}\left(Z\right)$ is estimated from the bootstrap samples, and the confidence interval is calculated as $\hat{Z}\pm 1.96\times \sqrt{\hat{D}\left(Z\right)}$ Alternatively, the quantile method is the most popular method for constructing confidence intervals. Let $Z_{\alpha }^{\ast }$ denote the α percentile of the bootstrap-estimated distribution for Z. Efron & Tibshirani (1994) used percentiles of the bootstrap distribution to estimate the 1 − α confidence interval for Z as (percentile bootstrap, PB) ${\displaystyle \left(Z_{\alpha /2}^{\ast },Z_{1-\alpha /2}^{\ast }\right).}$ Davison & Hinkley (1997) proceeded in a similar way, but a different formula was used to construct the confidence interval for Z, as (basic bootstrap, BB) ${\displaystyle \left(2\hat{Z}-Z_{1-\alpha /2}^{\ast },2\hat{Z}-Z_{\alpha /2}^{\ast }\right).}$

Bayesian estimation of the additive measure

Given a model, there are three steps involved in Bayesian analysis (Bolstad & Curran, 2016): (1) determining the marginal likelihood of the data, (2) specifying prior probabilities for the parameters, and (3) applying Bayes’ theorem to predictor variables and the outcome variable via Bayesian modeling.

Let Y denote the binary outcome variable. Consider the logistic regression of the form ${\displaystyle Y:Bernoulli\left(p\right)}$ where Bernoulli(⋅) is a Bernoulli distribution and p satisfies Eq. (1) or Eq. (2). The likelihood for the logistic regression is ${\displaystyle p\left(x|\theta \right)=\prod _{i}h{\left(x;\theta \right)}^Y{\left(1-h\left(x;\theta \right)\right)}^{1-Y}}$ where θ is a parameter (such as θ = (β₀, β₁, β₂, β₃, γ) in Eq. (1)) and $h\left(x;\theta \right)=\frac{1}{1+e^{-logit\left(p\right)}}$ is the predicted probability that Y is 1.

The prior probability is an unconditional probability that is assigned before any data are accounted for. There are three main categories of prior distributions: informative priors, weakly informative priors and noninformative priors. When a family of conjugate priors exists, the conjugate priors are chosen for computational efficiency. However, there is no conjugate prior for the likelihood function in logistic regression. When Bayesian inference is performed analytically, this makes the posterior distribution difficult to calculate except in very low dimensions. However, software such as OpenBUGS (Lunn et al., 2009), JAGS (Plummer, 2003) and Stan (Carpenter et al., 2017) allows these posteriors to be computed via simulation; hence, a lack of conjugacy is not a concern. This article employs noninformative or weakly informative priors p(θ) for all model parameters.

According to Bayes’ theorem, the joint posterior distribution of the model parameters is proportional to the product of the likelihood and priors, ${\displaystyle p\left(\theta |x\right)\propto p\left(\theta \right)\cdot p\left(x|\theta \right)=p\left(\theta \right)\prod _{i}h{\left(x;\theta \right)}^Y{\left(1-h\left(x;\theta \right)\right)}^{1-Y}.}$

This equation is too complex to solve analytically; therefore, Monte Carlo methods are often used to summarize the posterior distribution. A combination of Monte Carlo and Markov chain can result in effective sampling and evaluation. Thus, this posterior computation can be accomplished easily using Markov chain Monte Carlo (MCMC) methods with a random walk Metropolis algorithm. Since RERI, AP and S are functions of θ, the values of the three measures of interaction can be calculated directly when a sample is generated from the posterior distribution p(θ|x).

The uncertainties of the RERI, AP and S can be summarized by giving a range of values on the posterior probability distribution that includes 95% of the probability, (P_2.5, P_97.5), which is called a 95% credibility interval (calculated based on an equal-tailed interval). Similarly, in a Bayesian credibility interval, the null hypothesis is rejected if the credibility interval does not encompass the parameters of the null hypothesis. For example, there is an interaction if the credibility interval of an RERI does not contain 0. Notably, the confidence interval captures the uncertainty about the interval obtained (whether it contains the true value or not), and it cannot be interpreted as a probabilistic statement about the true parameter value. However, the credible interval captures the current uncertainty in the location of the parameter value and thus can be interpreted as a probabilistic statement about the parameter.

Preventive factors

Measures of interaction were developed for risk factors, so it is not appropriate to use them directly when a preventive factor exists. The recoded method, which can be seen as a choice of reference category, is used to address this issue (Knol et al., 2011). The use of the recoded method and the use of an updated model is a cumbersome and tedious process. In fact, when the reference category changes, the new logistic regression coefficients ${\theta }^{\prime }=\left({\beta }_0^{\prime },{\beta }_1^{\prime },{\beta }_2^{\prime },{\beta }_3^{\prime }\right)$ are completely dependent on the original logistic coefficients θ = (β₀, β₁, β₂, β₃), as resetting the reference category is equivalent to subtracting the new reference category coefficient from θ (Hosmer Jr, Lemeshow & Sturdivant, 2013). The relationships between θ′and θ for the three interaction measures are shown in Appendix S1A. In the Bayesian approach, the three measures of interaction are calculated directly by setting the reference category as that whose coefficient is the smallest (set β₀ = 0 in the dummy coding scheme). Thus, it is not necessary to update the model, and the uncertainty of the parameters will be obtained from the posterior sample. The R (R Core Team, 2023) function for computing new coefficients based on the lowest risk category as a reference category is provided in Appendix S1C.