Artificial neural networks modeling gene-environment interaction

Simulation study

Case-control data sets are generated using a two step design. First, underlying populations are simulated with a controlled prevalence of 10% and an overall sample size of five million observations. These populations carry the information of two marginally independent and randomly drawn factors – one biallelic locus and one continuous environmental factor – and a case-control status. The minor allele frequency is 30% to ensure sufficient cell frequencies in the final case-control data sets and it is assumed that the Hardy-Weinberg equilibrium holds. The environmental factor follows a continuous uniform distribution on the interval [0,100]. Depending on the genotype G and the environmental factor U, the case-control status is allocated through eight given theoretic risk models as introduced in the next subsection. Considering each theoretic risk model in a high and a low risk scenario, this results in sixteen underlying populations. As the second step, 100 case-control data sets are randomly drawn from all underlying populations for each analysis. Thus, for each analysis, mean values over 100 data sets are considered in sixteen situations. Three different sample sizes of 2,000 subjects (1,000 cases + 1,000 controls), 1,000 subjects (500 cases + 500 controls), and 400 subjects (200 cases + 200 controls) are used.

Artificial neural networks and logistic regression models are fitted to the data, i.e. separately to all 100 case-control data sets for each situation. A multilayer perceptron (MLP, see e.g. [15]) is chosen as neural network. It is briefly described in the Appendix. For neural networks, the genotype information is coded co-dominant, i.e. the genotype takes possible values 0, 1, and 2 representing the number of mutated alleles. The environmental factor is included in the analyses as continuous variable. For all data sets, six different network topologies, from zero up to five hidden neurons, are trained to avoid an overfitting of the data. For training purposes, the data set is always used as a whole. Each training process is replicated five times each with randomly initialized starting weights drawn from a standard normal distribution to enhance the chance that the training process stops within a global instead of a local minimum. The best trained neural network for each data set, i.e. the best network topology and the best repetition, is selected based on the Bayesian Information Criterion (BIC, [16]), which takes the number of parameters into account and penalizes additional parameters. Thus in each situation, 100 best neural networks predict the underlying risk model and the mean prediction can be used to evaluate the model fit (see below).

For comparison purposes, logistic regression models are fitted to the same data sets. The genotype is coded co-dominant counting the number of risk alleles and using two dichotomous design variables, one representing the heterozygous and one representing the homozygous mutated genotype. Five different models are used: the null model, three main effect models – containing only one or both main effects – and the full model – containing both main effects and one or two interaction terms depending on the genotype coding. For both coding approaches, the best model is selected based on BIC.

where g = 0,1,2 denotes the genotype and f(g,u ^′ ) refers to the theoretic risk model of the case-control data set and ${\widehat{f}}^{\left(k\right)}(g,u^{\prime })$ to the prediction of the kth case-control data set. The smaller $\sum _{g{u}^{\prime }}{E}_{g{u}^{\prime }}$ is, the better the mean model fit of neural networks or logistic regression models is since the estimated risk model and the theoretic risk model coincide for $\sum _{g{u}^{\prime }}{E}_{g{u}^{\prime }}=0$. To take variation into account, pointwise prediction intervals are calculated as empirical 95% intervals. In particular, for all u ^′ = 0, 0.1, 0.2,…,100 and g = 0,1,2 a prediction interval is determined as the interval $\left[\widehat{f}{(g,u^{\prime })}_{\left(3\right)};\widehat{f}{(g,u^{\prime })}_{\left(98\right)}\right]$, where $\widehat{f}{(g,u^{\prime })}_{\left(3\right)}$ and $\widehat{f}{(g,u^{\prime })}_{\left(98\right)}$ denote the 3rd ordered and the 98th ordered prediction, respectively.

Data generation and all analyses are done using R [17]. The package for training the MLP was implemented by our group and is published on CRAN [18].

Theoretic risk models

Two different types of theoretic risk models for gene-environment interactions are used, namely the models introduced by Amato et al. [14] and models mainly representing a masking effect of the involved locus as defined below. For all risk models, the kind of functional relationship between the penetrance and the environmental factor depends on the genotype information, i.e. the curve shape is in general different depending on the three genotypes. The relationship is defined on a population level, i.e. the penetrance function F :{0,1,2} × [0,100] → [0,1] with F(g u) = P(Y = 1|G = g U = u), where Y ∈ {0,1} denotes the case-control status, G ∈ {0,1,2} the genotype, and U ∈ [0,100] the environmental factor, only holds in the corresponding underlying population and has to be converted to f(g u) if a case-control data set is analyzed [10].

Risk models by Amato et al

Amato et al. [14] introduced four different risk models for analyzing gene-environment interactions: a genetic model, an environmental model, an additive model and an interaction model that are characterized by the following penetrance function

$\begin{array}{l}F(g,u)=\frac{1}{1+\text{exp}\{{\alpha }_g+{\beta }_g·u\}},\\ g=0,1,2;u\in \left[0,100\right].\end{array}$

The four models are defined as follows:

• the genetic model: α₁ ≤ α₂ ≤ α₃ and β₁ = β₂ = β₃ = 0,

• the environmental model: α₁ = α₂ = α₃ and β₁ = β₂ = β₃ ≠ 0,

• the additive model: α₁ ≤ α₂ ≤ α₃ and β₁ = β₂ = β₃ ≠ 0,

• the interaction model: α₁ = α₂ = α₃ and β₁ ≤ β₂ ≤ β₃.

To be able to fix the prevalence K of disease, we introduce an upper bound z that is determined such that the prevalence is equal to K = 0.1:

$F(g,u)=\frac{z}{1+\text{exp}\{{\alpha }_g+{\beta }_g·u\}},$

The values of α _g ,β _g , g = 0,1,2, and z used in this paper for two risk scenarios are given in Table 1. Figure 1 shows the theoretic risk models of a case-control data set f(g,u) for the high risk scenario.

Risk model		Risk scenario	Constant values αg, βg ( g = 0,1,2)	Constant values c, z
	Genetic model	High risk	${\alpha }_0=\frac{2}{3}·{\alpha }_1,{\alpha }_1=2.5,{\alpha }_2=\frac{4}{3}·{\alpha }_1$	z = 0.886
			β0 = β1 = β2 = 0
		Low risk	${\alpha }_0=\frac{2}{3}·{\alpha }_1,{\alpha }_1=1.25,{\alpha }_2=\frac{4}{3}·{\alpha }_1$	z = 0.390
			β0 = β1 = β2 = 0
	Environmental model	High risk	α0 = α1 = α2 = 7.5,	z = 0.200
			β0 = β1 = β2 = −0.15,
		Low risk	α0 = α1 = α2 = 3.75,	z = 0.200
Risk models by Amato et al. [14]			β0 = β1 = β2 = −0.075,
	Additive model	High risk	${\alpha }_0=\frac{2}{3}·{\alpha }_1,{\alpha }_1=7.5,{\alpha }_2=\frac{4}{3}·{\alpha }_1$,	z = 0.177
			β0 = β1 = β2 = −0.15,
		Low risk	${\alpha }_0=\frac{2}{3}·{\alpha }_1,{\alpha }_1=3.75,{\alpha }_2=\frac{4}{3}·{\alpha }_1$,	z = 0.178
			β0 = β1 = β2 = −0.075,
	Interaction model	High risk	α0 = α1 = α2 = 7.5,	z = 0.171
			β0 = 2 · β1, β1 = −0.15, β2 = 0.5·β1,
		Low risk	α0 = α1 = α2 = 3.75,	z = 0.169
			β0 = 2 · β1, β1 = −0.075, β2 = 0.5 · β1,
	Model 1	High risk (r = 0.150)		c = 0.05, z = 0.254
		Low risk (r = 0.075)
Risk model representing a masking effect of the genetic factor	Model 2	High risk (r = 0.150)		c = 0.05, z = 0.286
		Low risk (r = 0.075)
	Model 3	High risk (r = 0.150)		c = 0.075, z = 0.631
		Low risk (r = 0.075)
	Model 4	High risk (r = 0.150)		c = 0.075, z = 0.964
		Low risk (r = 0.075)

Constant values α _g , β _g (g = 0,1,2) c, and z used to determine the penetrance functions (minor allele frequency 30%).

Risk models representing a masking effect of the genetic factor

In addition, we define four theoretic risk models representing four types of gene-environment interactions where the gene mainly has a masking effect. The kind of functional relationship between the environmental factor and the penetrance again depends on the genotype information. The four theoretic risk models are described in detail in the following:

1. 1.
   The structure of the first risk model is given by the following penetrance function F :{0,1,2} × [0,100] → [0,1]

   $F(g,u)=\left\{\begin{array}{cc}\frac{z-c}{1+\text{exp}(-r(u-50\left)\right)}+c& \text{if}g=0\\ c& \text{if}g=1.\\ c& \text{if}g=2\end{array}\right.$
    
2. 2.
   The second risk model is defined by

   $F(g,u)=\left\{\begin{array}{cc}\frac{z}{1+\text{exp}(-r(u-50\left)\right)}& \text{if}g=0\\ c& \text{if}g=1.\\ 2c& \text{if}g=2\end{array}\right.$
    
3. 3.
   In the third risk model, the penetrance function is given by

   $F(g,u)=\left\{\begin{array}{cc}c& \text{if}g=0\\ c& \text{if}g=1.\\ \frac{z-c}{1+\text{exp}(-r(u-50\left)\right)}+c& \text{if}g=2\end{array}\right.$
    
4. 4.
   For the fourth risk model, the penetrance function is determined as follows:

   $F(g,u)=\left\{\begin{array}{cc}\frac{1}{2}c& \text{if}g=0\\ c& \text{if}g=1.\\ \frac{z-2c}{1+exp(-r(u-50\left)\right)}+2c& \text{if}g=2\end{array}\right.$
    

In each of these four models, r denotes the risk increase, c a baseline risk, and z an upper bound for the penetrance function. A risk increase of r = 0.150 indicates the high risk and a risk increase of r = 0.075 the low risk scenario, respectively. The baseline risk c and the upper bound z are again determined such that the prevalence of disease is equal to 10% for each situation. The values used in this paper are given in Table 1. Figure 2 shows the theoretic risk models of a case-control data set f(g,u) for the high risk scenario. Note that the gene has a main effect on the disease in risk models 2 and 4 only and that risk model 3 differs from risk model 1 only by different cell frequencies for the three genotype classes.

Real data application

To study the performance of a neural network in a real life situation, we applied this approach to a cross-sectional study dealing with a lifestyle induced complex disease. This application should serve as an example for the general practicability of our approach without describing the study from a subject point of view. The common effect of an SNP and a continuous environmental factor on a binary outcome is investigated while controlling for the effect of one binary confounder. The data set includes 138 cases and 1599 controls. As in the simulation study, neural networks with up to five hidden neurons are trained each five times with randomly initialized weights drawn from of a standard normal distribution and the best neural network is chosen based on BIC. The analysis is done once using the whole data set and once stratified by the confounding factor. For the stratified analysis, 95% bootstrap percentile intervals are calculated using 100 bootstrap replications [19].

Risk models by Amato et al.

A graphical comparison of the general modeling ability for neural networks and logistic regression models is shown in Figures 3 and 4 for the large sample size and the risk models introduced by Amato et al. [14]. In general, neural networks have a very good model fit compared to logistic regression models. Especially if the environmental factor has an effect (environmental model, additive model, interaction model), neural networks much better predict the underlying relationship between the genetic and the environmental factor with narrow prediction intervals that always include the true theoretic risk model. On the contrary, logistic regression models are only able to satisfactorily predict the genetic model. The sigmoid effects in the case that the environmental factor has an effect are not well represented in any situation and none of the prediction intervals include the theoretic risk model. This is true for both, logistic regression models with a co-dominant coding or for those using design variables for the genetic factor.

If the sample size decreases, the modeling ability becomes worse for neural networks as well as for logistic regression models (see Table 2). However, neural networks still show the best model fit if the environmental factor has an effect. The prediction intervals include the true underlying risk model in all but two situations (interaction model, n = 500 + 500, low risk scenario and interaction model, n = 200 + 200, high risk scenario, data not shown).

Models representing a masking effect of the genetic factor

The general modeling ability for the risk models representing a masking effect of the genetic factor is shown in Figures 5 and 6 for the large sample size. Neural networks have a very good model fit if the gene has a masking effect only (risk model 1 and 3). If the gene has an own main effect, the prediction is less accurate and the variance is much larger. Nevertheless, the prediction intervals include the true theoretic risk models in all situations. Logistic regression models with a co-dominant coding for the genotype are not able to capture the underlying structure of the gene-environment interaction. Constant or non-linear effects are not detected. If design variables are used for coding the genotype, the model fit is much better. However, the theoretic risk model is not included in the prediction interval in many situations. Additionally, the constant effects are only detected by combining the single predictions with either positive or negative slopes to an average prediction with zero slope. This is reflected by the wider ends of the prediction intervals. Moreover, the sigmoid effect is again not well represented in many of the investigated situations. This is especially true for the first and the second risk model.

With decreasing sample sizes, the model fit again becomes worse and the variance increases (data not shown). If the sample size is 500 + 500 subjects, neural networks again have the best model fit for the first three risk models in the high risk scenario. In the low risk scenario, this is only true for the first and the third risk model. A sample size of just 200 + 200 subjects leads to a considerably worse model fit of neural networks. In this situation, logistic regression models with design variables coding the genotype have the best model fit for the second and fourth risk model in both risk scenarios. Neural networks still have the best model fit if the gene has a masking effect only.

Real data application

The results for the real data application are shown in Figure 7 (whole data set) and Figure 8 (stratified analysis). A neural network without any hidden neuron is chosen as best neural network. It shows that the environmental factor in general has a decreasing effect on the disease risk and that the number of mutated alleles defines the slope of this risk decrease. If one or two mutated alleles are present, the risk is lower and the risk decrease is weaker as for the homozygous wildtype genotype. We also see a strong influence of the included binary confounding factor.

Artificial neural networks

The general idea of a multilayer perceptron (MLP) is to approximate functional relationships between co-variables and response variable(s). It consists of neurons and synapses that are organized as a weighted directed graph. The neurons are arranged in layers and subsequent layers are usually fully connected by synapses. Each synapse is attached by a weight indicating the effect of this synapse. A positive weight indicates an amplifying, a negative weight a repressing effect. Neural networks have to be trained using a learning algorithm to adjust the synaptic weights according to given data. The learning algorithm minimizes the deviation of predicted output and given response variable measured by an error function.

Data passes the MLP as signals. This process starts at the input layer consisting of all co-variables and a constant neuron and it stops at the output layer consisting of the response variable(s). Hidden neurons can be included between the input and output layer in several layers to increase the modeling flexibility. These hidden layers are not directly observable and cannot be controlled by data. See Figure 10 for an MLP with one hidden layer consisting of three hidden neurons that models the functional relationship between the locus G and the environmental factor U as co-variables and the case-control status Y as response variable.

where w₀, w _j , and w _ij , i = 0,…,n, j = 1,…,m, denote the weights including intercepts, x = (x₀x₁,…,x _n ) ^T the vector of all co-variables including a constant neuron x₀ and σ the activation function that maps the output of each neuron to a given range. MLPs are a direct extension of generalized linear models (GLM, [24]) and an MLP without hidden layer is algebraically equivalent to a generalized linear model with σ as inverse link function. In this case, trained weights and estimated regression coefficients coincide.

To train neural networks according to the case-control data sets, resilient backpropagation [25] as learning algorithm with cross entropy as error function and logistic function as activation function is used.