Identification of an additive interaction using parameter regularization and model selection in epidemiology

Estimating the additive interaction using logistic regression

In this method, we let y denote the binary outcome variable labeled by 0 and 1. Z is a vector of the potential confounders which is adjusted in the model and whose dimension, p, is dependent on the number of potential confounders.

Using the method proposed by Rothman (1986), the estimates of RERI, AP, and S can be estimated from the output of a multiple logistic regression

(1) $$logit(p)=\log ({\displaystyle \frac{p}{1-p})={\beta }_0+A{\beta }_1+B{\beta }_2+AB{\beta }_3+\mathrm{Z}\gamma }$$where $p=P(y=1|A,B,z)$ is the corresponding probability of outcome variable given the factors A, B and Z, and ${\beta }_i(i=0,1,2,3)$ and $\gamma \in R^p$ is the parameter vector corresponding to the two risk factors and confounders respectively. if ${\hat{\beta }}_i(i=0,1,2,3)$ are used to stand for the estimated logistic regression coefficients of ${\beta }_i(i=0,1,2,3)$, then the RERI, AP, and S can be estimated as

$$\widehat{RERI}=e^{{\hat{\beta }}_1+{\hat{\beta }}_2+{\hat{\beta }}_3}-e^{{\hat{\beta }}_1}-e^{{\hat{\beta }}_2}+1$$

$$\widehat{AP}={\displaystyle \frac{e^{{\hat{\beta }}_1+{\hat{\beta }}_2+{\hat{\beta }}_3}-e^{{\hat{\beta }}_1}-e^{{\hat{\beta }}_2}+1}{e^{{\hat{\beta }}_1+{\hat{\beta }}_2+{\hat{\beta }}_3}}}$$and

$$\hat{S}={\displaystyle \frac{e^{{\hat{\beta }}_1+{\hat{\beta }}_2+{\hat{\beta }}_3}-1}{e^{{\hat{\beta }}_1}+e^{{\hat{\beta }}_2}-2}.}$$

The confidence interval estimator can be calculated based on normal distribution (RERI and AP) or log-normal (S) assumption. To derive the necessary estimator of the variance for these three measures, we used the delta method which based on the first order approximation of a Taylor Series expansion (Hosmer & Lemeshow, 1992). The variance and covariance of ${\hat{\beta }}_i(i=1,2,3)$ are easily found from the logistic regression outcome in most packages.

As Rothman (1986) suggested before running the logistic regression, the formula can be simplified by combining the two risk factors into one variable with four levels. Specifically, given the two risk factors, A and B, the new variable is recoded as $A_i{B}_j(i,j=0,1)$, where the subscript indicates the corresponding value of the variable. Then, the new variable is seen as a categorical variable. Using dummy variable encoding and taking the $A_0{B}_0$ as reference level, the new variable is incorporated into the logistic regression model:

(2) $$logit(p)={\theta }_0+{\theta }_{1}I(A_1{B}_0)+{\theta }_{2}I(A_0{B}_1)+{\theta }_{3}I(A_1{B}_1)+z\gamma$$where $I(\cdot )$ is the indicator function. Let ${\hat{\theta }}_i(i=0,1,2,3)$ denotes the estimated coefficients of ${\theta }_i(i=0,1,2,3)$. In Eqs. (1) and (2) both models are saturated and the two estimations are equivalent. Furthermore, it is clear that ${\hat{\theta }}_1={\hat{\beta }}_1$, ${\hat{\theta }}_2={\hat{\beta }}_2$ and ${\hat{\theta }}_3={\hat{\beta }}_1+{\hat{\beta }}_2+{\hat{\beta }}_3$. In terms of $\hat{\theta }$’s, the RERI, AP, and S can be rewritten as

$$\widehat{RERI}=e^{{\hat{\theta }}_3}-e^{{\hat{\theta }}_1}-e^{{\hat{\theta }}_2}+1$$

$$\widehat{AP}={\displaystyle \frac{e^{{\hat{\theta }}_3}-e^{{\hat{\theta }}_1}-e^{{\hat{\theta }}_2}+1}{e^{{\hat{\theta }}_3}}}$$and

$$\hat{S}={\displaystyle \frac{e^{{\hat{\theta }}_3}-1}{e^{{\hat{\theta }}_1}+e^{{\hat{\theta }}_2}-2},}$$which yield the same point estimates given above. Based on the variances and covariance of ${\hat{\theta }}_i(i=1,2,3)$. The confidence intervals of these measures are calculated as described above, using the delta method.

When the confidence interval of RERI and AP contain a 0, or of S contains 1, this indicates a lack of the additive interaction between two risk factors. If these values are not present, there is an additive interaction.

Detect additive interaction based on regularized estimation of logistic regression

To establish notation, consider the logistic regression with sample size $n$, and assume that $\{y_i,A_i,B_i,z_i\}(i=1,2,\cdots ,n)$ with $n$ observations which have identical but independent distribution. To make estimation and expression simple, Eq. (2) is adopted and rewritten. Let $X_i=({\tilde{X}}_i,z_i)$, where ${\tilde{X}}_i$ is the ith row of the design matrix corresponding to the new combined variable, which is composed of zeros and ones denoting the combination of levels for two risk factors A and B. Logistic regression models show the conditional probability $p(X_i)=P(Y=1|X_i)$ by

$$logit(p(X_i))=\eta (X_i)$$with

$$\eta (X_i)={\theta }_0+{\tilde{X}}_i{\theta }_{ab}+z_i\gamma ={\theta }_0+X_i\tilde{\theta }$$where ${\theta }_0$ is the intercept and ${\theta }_{ab}=({\theta }_1,{\theta }_2,{\theta }_3{)}^T\in R^3$, $\tilde{\theta }=({\theta }_{ab}^T,{\gamma }^T{)}^T\in R^{3+p}$ is the parameter vector corresponding to the risk factors and confounders. For simplicity, let $\theta \in R^{4+p}$ denote the whole parameter vector, i.e., $\theta =({\theta }_0,{\tilde{\theta }}^T{)}^T.$

The standard maximum likelihood estimation for the logistic regression estimates the coefficients $\theta$ by solving the following equation

(3) $$\hat{\theta }=\arg mi{n}_{\theta }\{-l(\theta )\}$$where $l(\cdot )$ is the log-likelihood function:

$$l(\theta )=\sum _{i=1}^n{y}_i\eta (X_i)-\log [1+\mathrm{e}\mathrm{x}\mathrm{p}\{\eta ({\mathrm{X}}_i)\}].$$

The LASSO reduced the relative excess risk due to the additive interaction (Tibshirani, 2011) and a penalty terms can then be added to Eq. (3):

(4) $$\hat{\theta }=\arg mi{n}_{\theta }\{-l(\theta )\}\begin{array}{cc}& s.t.\begin{array}{cc}& |e^{{\theta }_3}-e^{{\theta }_1}-e^{{\theta }_2}+1|\end{array}\end{array}\le t$$where $t\ge 0$ is the tuning parameter which controls the amount of penalization. Solving $\hat{\theta }$ in Eq. (4) is equivalent to the “Lagrange” version of the problem

(5) $$\hat{\theta }=\arg mi{n}_{\theta }\{-l(\theta )+\lambda |e^{{\theta }_3}-e^{{\theta }_1}-e^{{\theta }_2}+1|\}$$where $\lambda \ge 0$. There is a one-to-one correspondence between $t$ and $\lambda$, whose values are regularly chosen by a model selection procedure such as K cross-validation (Fushiki, 2011), Akaike information criterion (AIC) or Bayesian information criterion (BIC) (Burnham & Anderson, 2004). If $\hat{\theta }(\lambda )$ minimizes Eq. (5), then it also solves Eq. (4) with $t=|e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1|$.

It should be noted that $e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1=0$ when $t$ equals 0 or $\lambda$ is large enough. This means that IC = 0, thus all three additive measures will equal 0 and it indicates the lack of additive interaction. Therefore, if there is no additive interaction, it hopes that $\lambda$ would be large in model selection. When $e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1\ne 0$, the RERI, AP, and S can be calculated and the additive interaction can be inferred that there is an additive scale interaction without requiring additional information, such as the confidence interval.

Detect additive interaction based on model selection

Although the regularized estimation of logistic regression can determine if there is an additive scale interaction, it depends on the choice of tuning parameter $t$ or $\lambda$. Whatever cross-validation or information criterion, it fixes the value through trial and error. In detail, given a collection of tuning parameter values in advance, for example $\lambda =\{{\lambda }_1,\cdots ,{\lambda }_l\}$. For each value of $\lambda$, Eq. (5) is solved. After $\lambda =\{{\lambda }_1,\cdots ,{\lambda }_l\}$ is executed, the best ${\lambda }^{\ast }$ can be determined based on K cross-validation or information criterion.

Model selection is key for determining whether there is an additive interaction (IC = 0). Specifically, two models are fitted; one (model $M_1$) is the logistic regression without any constraint, which is solving the Eq. (3); another (model $M_2$) is the logistic regression with the constraint $e^{{\hat{\theta }}_3}-e^{{\hat{\theta }}_1}-e^{{\hat{\theta }}_2}+1=0$, which is equivalent to solving Eq. (5) with large $\lambda$. In contrast to regularized estimation, it only fits two logistic models.

In regard to model selection, let $l(\hat{\theta })$ be the maximum value of the likelihood function for the model. For a parametric statistical model and within a likelihood-based inferential procedure, the fit of model can be assessed by its fitted global deviance (GDEV), defined as $GDEV=-2l(\hat{\theta })$. $M_1$ and $M_2$ may be evaluated with a fitted global deviance $GDE{V}_1$ and $GDE{V}_2$, and degrees of freedom $d{f}_1$ and $d{f}_2$, respectively. The two models, $M_1$ and $M_2$, are nested and may be compared using the likelihood ratio test or the generalized Akaike information criterion (GAIC) (Lv & Liu, 2014; Peers, 1971). The test statistic of the likelihood ratio is

$$\mathrm{\Lambda }=GDE{V}_2-GDE{V}_1$$which has an asymptotic ${\chi }_d^2$ distribution under the null hypothesis so the model $M_2$ is recommended, where $d=d{f}_2-d{f}_1$. The GAIC is obtained by adding a penalty $\kappa$ for each degree of freedom $df$ used in the model to the fitted deviance as follows:

(6) $$GAIC(\kappa )=GDEV+\kappa \times df.$$

Then the model with the smallest value of the criterion $GAIC(\kappa )$ is selected. AIC and BIC are used most often. $\kappa$ corresponds with $\kappa =2$ and $\kappa =\ln (n)$ in Eq. (6), respectively:

$$AIC=-2\cdot l(\hat{\theta })+2df,$$and

$$BIC=-2\cdot l(\hat{\theta })+\ln (n)df.$$

The AIC typically leads to overfitting in model selection, while the BIC leads to underfitting (McLachlan & Peel, 2000). HQ is another special case of $GAIC(\kappa )$ to consider (Burnham & Anderson, 2002),

$$HQ=-2\cdot l(\hat{\theta })+2\ln (\ln (n))\cdot df.$$

When $n\ge 15$, $\kappa =2\ln (\ln (n))$ is between 2 and $\ln (n)$, it is considered to be the compromise of AIC and BIC in the penalty terms.

There is only one constraint equation in model $M_2$, so the difference of the degree freedom of the two models equals one. Therefore, when the difference of the GDEV between the two models, $\ln (n)$ and $2\ln (\ln (n))$, is greater than 2 for AIC, BIC, and HQ, respectively, the logistic regression without constraint is chosen and there is an additive interaction between the two risk factors A and B.

Simulation

We conducted simulation studies whose additive interactions were present or absent based on a priori knowledge to assess the performance of detecting the additive interaction based on regularized estimation of logistic regression or model selection. For simplicity, only two risk factors A and B were considered, regardless of the potential confounders, Z, in all simulation scenarios. The data were generated from the model

$$Y\sim Bernoulli(p)$$

(7) $$logit(p)={\theta }_0+{\theta }_{1}I(A_1{B}_0)+{\theta }_{2}I(A_0{B}_1)+{\theta }_{3}I(A_1{B}_1)$$where $Bernoulli(\cdot )$ is Bernoulli distribution, and $\theta =({\theta }_0,{\theta }_1,{\theta }_2,{\theta }_3{)}^T\in R^4$ is the coefficients with intercept.

Regularized estimation

According to Zou, Hastie & Tibshirani (2007), the effective degree of freedom of the LASSO was the number of nonzero coefficients. However, unless the condition $e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1=0$ was met, the degree of the model $df$ would not change. This would make the information criterion fail to determine the $\lambda$. Thus, the value of $\lambda$ was tuned using generalization error which calculated GDEV on new data and can be seen as simplified of the cross-validation (Golub, Heath & Wahba, 1979).

Two datasets were generated with $n=600$ according to Eq. (7), where the coefficients were set to be $\theta =(0,1,1,1.49{)}^T$ and $\theta =(0,1,1,1.90{)}^T$, corresponding to models without and with interaction, respectively. The AP is 1/3 in interaction scenario, which can be regarded as weak synergy (Hosmer & Lemeshow, 1992). For each dataset, it was split into the training set and testing set with 400 and 200 cases, respectively. The coefficients were estimated using the training set, but the generalization error was calculated on the testing set.

Figure 1 showed the relation between the coefficients and tuning parameter $\lambda$, along with the penalty term $|e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1|$. The penalty term $|e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1|$ gradually converged to 0 as the tuning parameter $\lambda$ increased. However, the primary concern was whether there was an additive scale interaction rather than the path of the coefficients. Based on the tuning parameter $\lambda$ selected by generalization error, the final models matched the a priori information in two scenarios. Specifically, in the scenario in which there is no interaction, the penalty term $e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1=0$ is met, which indicated there was not additive interaction. While, in the scenario with an interaction, the penalty term $e^{{\hat{\theta }}_3(\lambda )}-e^{{\hat{\theta }}_1(\lambda )}-e^{{\hat{\theta }}_2(\lambda )}+1=0$ disobeyed, which indicated there was additive scale interaction.

Model selection

In order to assess the performance of detecting additive scale interactions based on model selection, a simulation study was conducted to compare the classical method, such as the delta and bootstrap methods. Model selection used AIC, BIC, and HQ. Both the delta and bootstrap methods employed the three indices: RERI, AP, and S, as denoted by RERI_D, AP_D, S_D, RERI_B, AP_B, and S_B.

For the simulation study, seven scenarios were concerned (Table 1). With the exception of scenarios S2 and S7, all scenarios were derived from the study by Assmann et al. (1996). Due to RERI = AP = 0 and S = 1 in scenario S4, there was no additive interaction. The model interaction occurred in all other models. Thus, those scenarios determined which models had strong synergy, weak synergy, no interaction, weak antagonism, or strong antagonism.

In each scenario, 1,000 samples were generated based on Eq. (7) with sample sizes 400, 600 and 1,000, respectively. The parameters $\theta$ were the logarithm of the ORs. For example, in scenario S1, parameters $\theta$ were set in Eq. (7) as

$$\theta =[0,\log (4),\log (5),\log (20){]}^T=(0,1.386,1.609,2.996{)}^T.$$

For each sample, the interaction detected model selection with $GAIC(\kappa )$ and the confidence intervals of RERI, AP, and S with delta and bootstrap methods. Then the methods were compared based on the percentage of times with no additive interaction.

For the bootstrap method (Efron, 1982), 300 bootstrap samples were tested separately within $Y=0$ and $Y=1$. Therefore, the number of $Y=0$ and $Y=1$ in the bootstrap samples were identical to the original sample.

The simulation results for models with no interaction are summarized in Table 2. The percentage was calculated as the number of times that model correctly identified in scenario S4 (without interaction), but it was the rate of times that model misidentified the other scenarios (with interaction). Thus, the correct model was chosen as that which was one minus the error rate in scenarios with interaction. Scenario S4 provided an assessment of the Type I error rate when testing for additive interactions using different methods, while the other scenarios provided an evaluation of the power of these methods when testing for additive interactions with varying effect sizes. Clearly, among all methods, the percentage became smaller and smaller as the sample size increased and/or the interaction effect was augmented. This was due to the fact that it was more likely to get significant results the larger the sample size and/or the effect size (Papoulis, 1990) in favor of the model with interaction. Additionally, when the interaction effect was weak it was due to a lack of power to detect the interaction.

The bootstrap method is recommended for use in various applied fields (Jalali et al., 2022; Miwakeichi & Galka, 2023; Kaity et al., 2023) as well as for interaction detection (Assmann et al., 1996). Whatever measures were used, the percentages of no interaction were very close in various scenarios. Therefore, the bootstrap method was used as a benchmark. In contrast, however, there was a lot of variation from the delta method, especially when the sample size was relatively small. For example, in scenario S2 for sample size of 400, the percentages of no interaction corresponding to RERI_D, AP_D, and S_D were 0.762, 0.199, and 0.345, respectively, and its standard deviation was 0.292. Using bootstrap, the percentages of no interaction corresponding to RERI_B, AP_B. and S_B were 0.302, the standard deviation almost equaled to 0. This indicated that there was often conflict among the results in detecting interaction via RERI_D, AP_D, and S_D. However, based on model selection, due to $e^{{\hat{\theta }}_3}-e^{{\hat{\theta }}_1}-e^{{\hat{\theta }}_2}+1=0$, there were consistent results regardless of which measure used.

Using the bootstrap method as the criterion, the S_D had the smallest deviation and RERI_D had the largest deviation in almost all scenarios. Specifically, S_D was close to S_B, but RERI and AP were not. In contrast to the bootstrap method, RERI_D was favored without interaction, regardless of whether it was synergy or antagonism interaction. However, AP_D favored the model with a synergy interaction and without an antagonistic interaction.

Obviously, when $\kappa$ increased, the $GAIC(\kappa )$ also increased and tended to favor the model without interaction. Since $\kappa =2$, $2\ln (\ln (n))$ and $\ln (n)$ for AIC, HQ, and BIC, respectively, for each scenario. The percentage of no interactions were rising, as shown in Table 2. It is noted that the results of HQ were almost the same as the results of the bootstrap method. The difference between the results of HQ and the bootstrap method were less than 2% in all scenarios. In contrast, there was a smaller percentage of no interaction for AIC, but larger for BIC. This meant that the AIC tended to overestimate the interaction, while the BIC tended to underestimate the interaction (McLachlan & Peel, 2000).

Real data

The data came from a case-control study of oral cancer, kindly supplied by Rothman & Keller (1972). The two factors were smoking and alcohol use. There were a total of 458 participants who were male veterans under the age of 60s. The study investigated the effects of smoking and alcohol use on oral cancer. The distribution of data is shown in Table 3. It was used to illustrate their methods on how to calculate the confidence interval for measures of interaction by several authors (Hosmer & Lemeshow, 1992; Zou, 2008; Richardson & Kaufman, 2009; Chu, Nie & Cole, 2011). The confidence intervals were different among their methods. However, their results consistently showed that there was not an additive scale interaction.

Using the interaction detection based on model selection, two logistic models were employed. $M_u$ and $M_c$ indicated which model had the constraint $e^{{\hat{\theta }}_3}-e^{{\hat{\theta }}_1}-e^{{\hat{\theta }}_2}+1=0$. The result is summarized in Table 4. The GDEV of model $M_u$ and $M_c$ were 605.930 and 607.727, respectively. Thus, the difference between them was 1.803, which was less than two. The model without an interaction was favored, regardless of which criterion was used (AIC, BIC, or HQ). Therefore, consensus had to be reached with other authors. Remarkably, the coefficients of model $M_c$ without an additive interaction met the following equation:

$$\mathrm{I}\mathrm{C}=e^{2.306}-e^{1.882}-e^{1.497}+1=0,$$which would make the interpretation of the result simple. Additionally, we validated our algorithm in two other cases: a hypertension case (Zou, 2008) and a congenital heart disease case (Nie et al., 2016); the results were consistent with the authors’ analysis (Tables S1–S4).