G × E model definition

A standard G × E model can be represented in the following way:

(1)$${\bi y} = \beta _0 + \beta _e{\bi e} + \beta _g{\bi g} + \beta _{eg}{\bi eg} + {\bf \varepsilon}, \, $$

where yis a vector representing the n observed outcomes, β ₀ is the intercept, β _e the regression weight for the environment main effect (e), β _g the regression weight for the gene main effect (g), and β _eg is the regression weight for the product of environment and genes (eg), which represents the G × E interaction. For a fitted model, we use the following formulation:

(2)$$E\lsqb {{\bi y} \vert {\bi e}\comma \,{\bi g}} \rsqb = {\hat \beta} _0 + {\hat \beta} _e{\bi e} + {\hat \beta} _g{\bi g} + {\hat \beta} _{eg}{\bi eg}.$$

The hat (^) is used on top of parameters to represent that these are estimated parameters rather than the original ones. This notation will be used through the paper.

Simple slopes and RoS

Assuming a single binary genetic variant (e.g., short vs. long 5-HTTLPR allele), a G × E model has two environment-predicting outcome regression lines: one representing individuals with, putatively, low environmental sensitivity and another representing individuals with, putatively, high environmental sensitivity. The traditional approach for distinguishing between differential susceptibility, diathesis-stress, and vantage sensitivity is to determine at which values of the environment the lines differ significantly from each other.

To this end, one must first reparametrize a G × E model as an explicit function of the genotype:

$$\displaylines{E\lsqb {{\bi y} \vert {\bi e}\comma \,{\bi g}} \rsqb = {{\hat \beta}} _0 + {{\hat \beta}} _e{\bi e} + {{\hat \beta}} _g{\bi g} + {{\hat \beta}} _{eg}{\bi eg} \cr \qquad \qquad \hskip5pt = \lpar {{{\hat \beta}}_0 + {{\hat \beta}}_e{\bi e}} \rpar + ({{\hat \beta}} _g + {{\hat \beta}} _{eg}{\bi e}){\bi g} \cr \hskip-23pt = {{\hat w}}_0 + {{\hat w}}_1{\bi g}.} $$

We call ŵ ₀ the simple intercept and ŵ ₁ the simple slope. The simple slope represents the slope as a function of g. This means that if ŵ ₁ is significantly different from zero, then the slopes for different values of g (e.g., the slopes for g = 0 and g = 1 if g is dichotomous) are significantly different from one another.

The variance of the simple slope is

$$\displaylines{Var\lsqb {{{\hat w}}_1 \vert {\bi e}} \rsqb = Var\lpar {{{\hat \beta}}_g} \rpar + Var\lpar {{{\hat \beta}}_{eg}{\bi e}} \rpar + 2Cov\lpar {{{\hat \beta}}_g\comma \,{{\hat \beta}}_{eg}{\bi e}} \rpar \; \; \cr \hskip3.6pc = Var\lpar {{{\hat \beta}}_g} \rpar + {\bi e}^2Var\lpar {{{\hat \beta}}_{eg}} \rpar + 2{\bi e}Cov\lpar {{{\hat \beta}}_g\comma \,{{\hat \beta}}_{eg}} \rpar .} $$

Using standard regression asymptotic theory (Aiken et al., Reference Aiken, West and Reno1991), it can be shown that

$$\displaystyle{{{{\hat w}}_1} \over {\sqrt {Var\lsqb {{{\hat w}}_1 \vert {\bi e}} \rsqb }}} \sim t_{N-5},$$

where t _N-5 is a student's t distribution with N – 5 degrees of freedom (df) and N is the number of observations.

Based on the above equations, it is possible to test whether the simple slope is significantly different from zero, and thus determine if the slopes for different values of g are significantly different from one another at fixed values of e. However, this does not tell us at which values of e is ŵ ₁ significant. The Johnson–Neyman technique (Johnson & Fay, Reference Johnson and Fay1950) uses this formula, in a backward fashion, to find the range of values of e where the simple slope is significantly different from zero. This range corresponds to the RoS, and it can be defined by a lower bound L and an upper bound U. When e is larger than L and smaller than U, the slopes for different values of g do not differ significantly, whereas outside these bounds they do.

Kochanska et al. (Reference Kochanska, Kim, Barry and Philibert2011) suggested using RoS (Aiken et al., Reference Aiken, West and Reno1991; Hayes & Matthes, Reference Hayes and Matthes2009; Preacher et al., Reference Preacher, Curran and Bauer2006) to help differentiate diathesis-stress from differential susceptibility. When the RoS are within the observable ranges of values for e, the data are considered to reflect differential susceptibility. When only the lower bound L of the RoS is within the observable range, the data is considered to reflect diathesis-stress. When only the upper bound U of the RoS is within the observable range, the data is considered to reflect vantage sensitivity. If neither L nor U is within the observable range, results are inconclusive.

Confirmatory and competitive models

As shown above, the RoS approach focused entirely on the lower and upper bounds of the RoS. An alternative point of view is to focus instead on the crossover point, the point where the slopes for the different values of the environment cross. It was shown (Aiken et al., Reference Aiken, West and Reno1991; Widaman et al., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) that the crossover can be found using the following formula:

(3)$$c = -\displaystyle{{\beta _g} \over {\beta _{eg}}}.$$

Similarly to the RoS approach, if the complete 95% confidence interval of the crossover point is within the observable range of the environment, it is suggestive of differential susceptibility. To obtain a confidence interval for the crossover point, Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) suggested fitting a model where the crossover point c is an unknown parameter to be estimated. The authors showed that the standard G × E formulation (1) can be reparametrized as

(4)$${\bi y} = \beta _0 + \beta _e\lpar {{\bi e}-c} \rpar + \beta _{eg}\lpar {{\bi e}-c} \rpar {\bi g} + {\bf \varepsilon}.$$

Although β _g does not appear in this equation, Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) have shown that this parametrization is equivalent due to the inclusion of the crossover point c.

Using nonlinear regression estimation, Model 5 can be fitted and the crossover point can be estimated directly along with its confidence intervals. It is important to note that the estimate of the crossover point is not normally distributed (Lee, Lei, & Brody, Reference Lee, Lei and Brody2015; Marsaglia, Reference Marsaglia1965), suggesting that the resulting confidence interval may be biased. However, it has been shown that when the sample size is large enough (N > 500), the crossover point is approximately normally distributed (Lee et al., Reference Lee, Lei and Brody2015; Marsaglia, Reference Marsaglia1965). In addition, even when the sample size is small, the nonnormality of the crossover point does not lead to significant problems in the case of the competitive-confirmatory approach, as it does not rely exclusively on confidence intervals to distinguish between differential susceptibility and diathesis-stress; this is in contrast to the RoS approach. Instead, Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) constructed confirmatory and competitive models of differential susceptibility and diathesis-stress to determine which model best fits the data. In accordance, the purpose of the confidence interval of the crossover point is only to verify that the crossover point is within the observable range of the environment when a differential susceptibility model is deemed to be the best fit-wise. To prevent any confusion, we will use the notation c for the true crossover point, ĉ for the estimated crossover of the differential susceptibility models, c _low for the crossover point of the vantage sensitivity models, and c _high for the crossover point of the diathesis-stress models.

Although the original paper by Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) focused only on the diathesis-stress and differential susceptibility models, we extend this approach in the current effort by also presenting the vantage sensitivity models. Furthermore, we extend all concepts to situations where multiple environmental variables contribute to an environmental score, and where multiple genetic variables contribute to a genetic score. Accordingly, there are six models tested in our extension to multidimensional situations:

1. 1. weak vantage sensitivity (β _e estimated, c _low = min(e))

2. 2. strong vantage sensitivity (β _e = 0, c _low = min(e))

3. 3. weak differential susceptibility (β _e estimated, c estimated)

4. 4. strong differential susceptibility (β _e = 0, c estimated)

5. 5. weak diathesis-stress (β _e estimated, c _high = max(e))

6. 6. strong diathesis-stress (β _e = 0, c _high = max(e))

All six models are constructed using Equation 4. The models of differential susceptibility (Models 3 and 4) need to estimate an extra parameter, namely, the crossover point; therefore, Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) suggested using the Akaike information criterion (Akaike, Reference Akaike, Akaike, Parzen, Tanabe and Kitagawa1998) or the BIC (Schwarz, Reference Schwarz1978) to evaluate the quality of fit.

Latent genetic and environmental score models

Instead of considering a linear regression on a single observed genetic variant and environment, as in Equation 1, Jolicoeur-Martineau et al. (Reference Jolicoeur-Martineau, Wazana, Székely, Steiner, Fleming, Kennedy and Greenwood2017) instead considered a linear regression on latent genetic (g) and environmental (e) scores and their interaction. These latent scores are defined as the weighted sum of their corresponding variables (genetic or environmental). The goal thus becomes not only to estimate the parameters of the G × E but also to estimate the weights of the genetic and environmental variables that make up the latent variables. Jolicoeur-Martineau et al. (Reference Jolicoeur-Martineau, Wazana, Székely, Steiner, Fleming, Kennedy and Greenwood2017) defined their model in three parts: the genetic score g, the environmental score e, and the main model.

(5)$$\eqalign{{\bi g} & = \mathop \sum \limits_{\,j = 1}^k p_j {\bi g}_j \cr {\bi e} &= \sum \limits_{l = 1}^s q_l {\bi e}_l \cr & y = \beta _0 + \beta _e {\bi e} + \beta _g {\bi g} + \beta _{eg} {\bi eg} + {\bf \varepsilon}\comma} $$

where p = (p ₁, p ₂, … , p _k) is a vector of parameters for the k genetic variables g₁, g₂, … , g_k and q = (q ₁, q ₂, … , q _s) is a vector of parameters for the s environmental variables e₁, e₂, … , e_s.

Note that there are infinitely many possibilities for p and q that lead to the same model. For example, β _e, β _eg with Cp, where C is a constant, leads to the exact same model as C β _e, Cβ _eg with p:

$$\eqalign{& y = \beta _0 + \beta _e {\bi e} + \beta _g \sum \limits_{\,j = 1}^k (Cp_j) {\bi g}_j + \beta _{eg}{\bi e} \sum \limits_{\,j = 1}^k (Cp_j) {\bi g}_j + {\bf \varepsilon} \cr & = \beta _0 + \beta _e {\bi e} + C\beta _g \sum \limits_{\,j = 1}^k p_j {\bi g}_j + C\beta _{eg}{\bi e} \sum \limits_{\,j = 1}^k p_j {\bi g}_j + {\bf \varepsilon} \cr & = \beta _0 + \beta _e{\bi e} + (C\beta _g){\bi g} + (C\beta _{eg}){\bi eg} + {\bf \varepsilon}.} $$

To prevent infinite possibilities for p and q, the following restrictions are added:

(6)$$\matrix { {\mathop \sum \limits_{\,j = 1}^k \vert {\,p_j} \vert = 1\comma \,} \hfill \cr {\; \mathop \sum \limits_{l = 1}^s \vert {q_l} \vert = 1.} \hfill \cr} $$

This restriction also aids interpretation because the weights of the genetic and environmental scores represent the relative contribution of the individual genetic variants and environmental variables to their respective composites.

To estimate the weights of the genetic score, the environmental score, and the main models, an alternating optimization approach is used. We first initialize p to (1/k, 1/k, … , 1/k) and q to (1/s, 1/s, … , 1/s), that is, we assume that all genetic and environmental variables are equally weighted; this is equivalent to taking the average of all the genetic or environmental variables. Then, the approach comprises three steps: (a) estimating β ₀, β _e, β _g, β _eg assuming that p and q are known, (b) estimating p assuming that q and β ₀, β _e, β _g, β _eg are known, and (c) estimating q assuming that p and β ₀, β _e, β _g, β _eg are known (Jolicoeur-Martineau et al., Reference Jolicoeur-Martineau, Wazana, Székely, Steiner, Fleming, Kennedy and Greenwood2017). This is done in iterative steps until convergence, at each step ensuring that the parameters of the genetic and environmental scores sum to one in absolute values (constraints from Equation 6). Table 1 shows an example of alternating optimization being used to estimate a LEGIT G × E model with four genes and one environmental variable. The approach is discussed in greater details in our previous methodological paper (Jolicoeur-Martineau et al., Reference Jolicoeur-Martineau, Wazana, Székely, Steiner, Fleming, Kennedy and Greenwood2017).

Implementations in R and SAS are freely available online on CRAN (cran.r-project.org/web/packages/LEGIT) and Github (github.com/AlexiaJM/LEGIT).

Combining RoS with alternating optimization

Combining RoS with alternating optimization is a simple step, whereby one can apply RoS directly to the main model of Equation 5 assuming that p and q are known.

Combining confirmatory with alternating optimization

To combine the competitive-confirmatory approach with alternating optimization, Model 6 can be reparametrized and reformulated to include a crossover point, thus enabling the testing of the interaction, in the following way:

(7)$$\eqalign{&{\bi g} = \sum \limits_{\,j = 1}^k p_j {\bi g}_j \cr & {\bi e} = \sum \limits_{l = 1}^s q_l {\bi e}_l \cr & y = \beta _0 + \beta _e({\bi e}-c) + \beta _g {\bi g} + \beta _{eg}({\bi e}-c) {\bi g} \,+{\bf \varepsilon}, } $$

This new formulation corresponds to the competitive-confirmatory model (Belsky et al., Reference Belsky, Pluess and Widaman2013; Widaman et al., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) shown in Equation 4. Using this formulation, estimating the crossover point c would require the use of a nonlinear regression; however, by simply adding a negative intercept to the environmental score, this equation can be reparametrized in a way that each part remains a linear model, as shown below:

(8)$$\eqalign{&{\bi g} = \sum \limits_{\,j = 1}^k p_j {\bi g}_j \cr &{\bi e}^{\prime} = -c + \sum \limits_{l = 1}^s q_l {\bi e}_l \cr & y = \beta _0 + \beta _e {\bi e}^{\prime} + \beta _g {\bi g} + \beta _{eg}{\bi e}^{\prime}g +{\bf \varepsilon} .} $$

Accordingly, the alternating optimization approach can be easily modified to include a crossover point and test for the six different hypotheses mentioned above (i.e., diathesis-stress, differential susceptibility, or vantage sensitivity; weak or strong). See Appendix A for more details on how to adapt the alternating optimization algorithm for competitive-confirmatory testing and Appendix B for special considerations regarding testing for the presence of a G × E in case of multiple genes and environments.

Simulation setup

We have proposed extensions of the concepts of RoS and competitive-confirmatory approaches to multidimensional environmental and genetic scores. To test the performance of these extensions with one or multiple genes and environments, we examined how often these two approaches can correctly determine the pattern of interaction (i.e., diathesis-stress, vantage sensitivity, or differential susceptibility) under various scenarios. We simulated 100 samples for each of the six different models (representing the three hypotheses, each as weak or strong) and report the average accuracy of all models. We defined accuracy as the percentage of correctly assigned patterns of interactions. In addition, we simulated 100 samples of the model without a G × E (genes and environments only) to estimate the false positive rate (i.e., the percentage of models not assigned as having no evidence of G × E when there is no actual G × E).

We present one scenario with a single genetic variant and single environment (traditional G × E model) and one scenario with four genetic variants and three environmental variables. In both scenarios, these variables are sampled from the following distributions:

$$\eqalign{& {\bi g}_j \sim {\rm Bernouilli} \lpar {.30} \rpar \comma \, \cr & {\bi e}_l \sim \rm Beta \lpar {2\comma \,\beta} \rpar .} $$

Table 1. Example of alternating optimization being used to estimate a simple G × E model with four genetic variables and one environment variable; E[y] = β ₀ + β _gg + β _ee + β _egeg with a single e and g = p ₁g ₁ + p ₂g ₂ + p ₃g ₃ + p ₄g ₄

Given that the vast majority of environmental measures in psychological/epidemiological research are ordinal (e.g., Likert scale) and thus bounded by a minimum and maximum value, we chose a beta distribution for the environmental factors. Note that the beta distribution is bounded by 0 and 1. We present two scenarios, one with β = 2, which leads to a fully symmetric normal-like distribution, and one with β = 4 for a highly left-skewed distribution. Left-skewed variables are very frequently used as environmental factors (e.g., socioeconomic status, income, depressive symptoms, etc.), which justifies using this type of distribution. A priori, we theorized that it might be harder to confirm a diathesis-stress model with crossover at the maximum score of the environmental variable when very few individuals have values close to the maximum.

In all models, we let β ₀ = 3 and β_eg = 2. The weak models have β _e = 1, while the strong models have β _e = 0. The crossover point is set to 0 for the vantage sensitivity models and 1 for the diathesis-stress models. We present two scenarios for the choice of the crossover of the differential susceptibility models: c = .50, the easiest case given that it is right in the middle of the observable range of the environmental score, and c = .25, a more difficult case given that it is closer to the minimum possible value of the environmental quality score.

Assuming a Gaussian error term with variance σ (i.e., ε ~ Normal(0, σ)), as in the case of traditional linear regression assumptions, we set up realistic simulations by using three different choices of σ to represent small, medium, and large effect sizes. For the models with only one genetic and environmental variable, we set the variance of the error term so that the R ² was .05, .10, and .15 for small, medium, and large effect sizes, respectively. For the models with four genetic and three environmental variables, we set the variance of the error term so that the R ² was .10, .20, and .40 for small, medium, and large effect sizes, respectively. Although the values for the scenario including multiple genes and environments may seem large, they represent realistic values that have been previously observed using the alternating optimization approach (Jolicoeur-Martineau et al., Reference Jolicoeur-Martineau, Wazana, Székely, Steiner, Fleming, Kennedy and Greenwood2017), because including multiple genetic variants and environmental measures in a model tends to greatly increase its predictive power. In addition, we tested these scenarios under different sample sizes: N = 100, 250, 500, 1000, and 2000, representing very small, small, medium, large, and very large sample sizes, respectively.

To summarize, we examine the following scenarios:

1. 1. The interaction between (a) a single genetic and environmental variable or (b) three genetic and four environmental variables

2. 2. Assuming small, medium, or large effect sizes

3. 3. Under sample sizes: N = 100, 250, 500, 1000, or 2000

4. 4. Assuming symmetric Beta(2,2) or skewed Beta(2,4) distribution of the environmental measures

5. 5. Fixing the crossover of differential susceptibility at either .50 or .25

See Table 2 for a list of all the parameters used.

Table 2. Parameters used in the different simulations, where E[y] = β ₀ + β _e (e – c) + β _eg (e – c)g, for N = 100, 250, 500, 1000, 2000 and small, medium, and large effect sizes

Note: Small, medium, and large effect sizes refers to R ² = .05, .10, and .15 in the one gene and one environment case and to R ² = .10, .20, and .40 in the four genes and three environments case.

Hypothesis testing

As originally presented by Belsky et al. (Reference Belsky, Pluess and Widaman2013), the competitive-confirmatory approach is primarily concerned with the nature of the G × E, even when the G × E is not significant. Recently, Belsky and Widaman (Reference Belsky and Widaman2018) suggested using the competitive-confirmatory approach only if the F-ratio of the G × E is greater or equal to 1; this strategy prevents trying to fit ill-conditioned models with near zero interaction effect, but it may not be penalizing enough to bring the rate of false positive to 5% or lower. A natural approach to minimize the presence of false positives is to additionally test for models without a G × E term. Thus, in addition to the six G × E models of interest, we also examine the following four models: (a) Intercept only, (b) gene(s) only, (c) environment(s) only, and (d) gene(s) and environment(s) only.

If any of the four models without an interaction had the lowest BIC, we classified the interaction as “no evidence of G × E.” Otherwise, we classified the interaction as “differential susceptibility” if (a) the weak or strong differential susceptibility models had the lowest BIC and (b) the 95% interval of its estimated crossover point was within the observable bounds of the environmental score. If one of these conditions was not met, we classified the interaction based on which of the remaining four models (i.e., weak/strong vantage sensitivity and diathesis-stress models) had the lowest BIC.

For the RoS approach, we classified the pattern of interaction as reflecting “differential susceptibility” when both lower and upper bounds of the RoS were within the observable range, “vantage sensitivity” when only the upper bound of the RoS was within the observable range, “diathesis-stress” when only the lower bound of the RoS was within the observable range, or “no evidence of G × E” when neither of the bounds of the RoS were within the observable range. The lower and upper bounds were determined using the 95% confidence interval of the simple slope.

With multiple genes and environments, we found that the .05 α level (95% confidence intervals) lead to extremely high rates of false positives (80–97%). To remedy this issue, we decreased the α level until we reached a false positive rate of 15% or lower in every scenario. The resulting α level was .0001 (99.99% confidence intervals). Given that this is very conservative, it inevitably resulted in lower accuracy levels. However, this step was necessary in order to prevent the extremely high false positive rates.