A V Bulinski

O A Butkovsky

A P Shashkin

P A Yaskov

Statistical methods of SNP data analysis with applications

Keywords: Genetic data statistical analysis, multifactor dimensionality reduction, logic regression, random forests, stochastic gradient boosting, independent rule, single nucleotide polymorphisms, external factors, ischemic heart disease, myocardial infarction, supercomputer AMS 2010 classification: 92B15, 92D10, 65C20

Various statistical methods important for genetic analysis are considered and developed. Namely, we concentrate on the multifactor dimensionality reduction, logic regression, random forests and stochastic gradient boosting. These methods and their new modifications, e.g., the MDR method with "independent rule", are used to study the risk of complex diseases such as cardiovascular ones. The roles of certain combinations of single nucleotide polymorphisms and external risk factors are examined. To perform the data analysis concerning the ischemic heart disease and myocardial infarction the supercomputer SKIF "Chebyshev" of the Lomonosov Moscow State University was employed.

Introduction

The detection of genetic susceptibility to complex diseases (such as cardiovascular, oncological ones etc.) has recently drawn much attention in many leading research centers, see, e.g., [START_REF] Bock | Computational epigenetics[END_REF] and [START_REF] Szymczak | Machine Learning in Genome-Wide Association Studies[END_REF]. According to the forecast of the World Health Organization (www.who.int), in 2030 the deaths related to cardiovascular diseases will exceed 23 millions (this year about 17 millions), the oncological diseases will take the lives of more than 11 millions of our planet inhabitants and at least 2 millions of people will be the victims of the diabetes. Thus this research domain is important since one expects to provide for each person the prophylactic measures and medical treatment taking into account his/her genetic peculiarities which increase the risk of some diseases and protect from the others, see, e.g., [START_REF]Bioinformatics From Genomes to Therapies[END_REF]. Individual's DNA variations are typically described in terms of single nucleotide polymorphisms (SNP), i.e. the fragments of genetic code where a nucleotide change is possible. For more details see, e.g., [START_REF] Venter | The Sequence of the Human Genome[END_REF]. The first examples of genetically based diseases (e.g., sicklemia) were related with a single mutation. Contrariwise many hard diseases such as diabetes, Alzheimer's disease and others have a complex character as they 1 Dept. of Mathematics and Mechanics, Lomonosov Moscow State University, Moscow, 119991, Russia. 2 LPMA, Université Pierre et Marie Curie (Paris-6), 4, place Jussieu 75252, Paris Cedex 05, France. 3 Steklov Mathematical Institute, Gubkina str. 8, Moscow, 119991, Russia.

1 can be provoked by mutations in different parts of the DNA code which are responsible for the formation of certain types of proteins. Quite a number of recent studies (see, e.g., [START_REF] Dai | SHARE: an adaptive algorithm to select the most informative set of SNPs for candidate genetic association[END_REF], [START_REF] Schwender | Testing SNPs and sets of SNPs for importance in association studies[END_REF] and [START_REF] Schwender | Empirical Bayes analysis of single nucleotide polymorphisms[END_REF]) support the paradigm that the increasing risks of complex diseases can be explained by combinations of certain SNP whereas separate mutations have no dangerous effects.

Thereupon it should be mentioned that there existed a longstanding demand for statistical analysis of biological and medical data. However, only in the first part of the 20th century, due to the classical contributions by K.Pirson, R.Fisher, H.Cramér, A.Kolmogorov, N.Smirnov, A.Wald and other prominent statisticians, the essential progress was achieved both in theory and applications. The methods developed were sufficient, e.g., for investigation of the efficiency of new medicaments. The situation has changed radically at the beginning of the 2000's when the laboratory methods of DNA analysis provided the data related to the personal human code structure. The achievements in decoding of the human genome have led to formation of vast data bases in the frameworks of the International Research projects, see, e.g., GAW16 [17] and HapMap [20]. Note also that software engineering plays an important role in such studies, see, e.g., [START_REF] Karchin | Next generation tools for the annotation of human SNPs[END_REF] and [START_REF] Zhang | Software engineering the mixed model for genome wide assosiation studies on large samples[END_REF]. The cost of genomic analysis has fallen considerably in the last 10 years, allowing to collect large volumes of genetic data for genetic mapping of complex diseases. However statistical problems arising here require new methods of inference rather than classical ones. Indeed, the modern statistical models involve huge number of variables, parameters, hypotheses etc., while the sample size is usually moderate (several hundred or sometimes several thousand of observations, see, e.g., [START_REF] Kohara | Identification of Hypertension-Susceptibility Genes and Pathways by a Systemic Multiple Candidate Gene Approach: The Millennium Genome Project for Hypertension[END_REF]). The sample design is limited both by costs of analysis which are still high and by difficulties due to the sample selection. In particular, the ethnic homogeneity should be taken into account, as well as the influence of external risk factors such as obesity, smoking etc.

To perform reliable statistical inference, it is necessary to apply new powerful tools developed in high-dimensional statistics, artificial intelligence, information retrieval, econometrics etc. Some of them have been adapted and further generalized in numerous papers by biostatisticians. Among the most important SNP analysis methods are the multifactor dimensionality reduction (MDR), logic regression (LR), random forests (RF) and stochastic gradient boosting (SGB). All approaches based on these methods do not impose any strong restrictions on the dependence structure of variables under consideration (apart from independence and identical distribution of observations within certain groups). Thus a broad class of statistical models is defined and the model providing the best out-of-sample fit is selected.

If one deals with too many parameters, overfitting is likely to happen, i.e. the estimated parameters depend too much on the given sample. As a result the constructed estimators give poor prediction on new data. On the other hand, application of a very sophisticated model may not capture the studied dependence structure of various factors efficiently. However the trade-off between the model's complexity and its predictive power allows to perform reliable statistical inference via new model validation techniques (see, e.g., [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF] and [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]). The main tool of model selection is the cross-validation, see, e.g., [START_REF] Winham | A comparison of internal validation techniques for multifactor dimensionality reduction[END_REF]. Its idea is to estimate parameters by involving only a part of the sample (training sample) and afterwards use the remaining observations (test sample) to test the predictive power of the obtained estimates. Then an average over several realizations of randomly chosen training and test samples is taken, see [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF].

There are two closely connected research directions in genomic statistics. The first one is aimed at the disease risk estimation when the genetic portrait of a person is known (in turn this problem involves estimation of disease probability and classification of genetic data into high and low risk domains, see, e.g., [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]). The second trend is to identify relevant combinations of SNPs having the most significant pathogenic (or, in other way, protective) influence. Both directions are presented in this paper. Moreover, the authors propose further development of various statistical methods and apply them to study of the risks of cardiovascular diseases. For this purpose the new software concerning the employment of the mentioned statistical methods was designed and used.

Due to high-dimensionality of data many numerical procedures based on the above mentioned statistical methods are very time consuming. The authors are grateful to the Chancellor of the Lomonosov Moscow State University (MSU) Professor V.A.Sadovnichy and to the Deputy Director of the MSU Research Computing Center Professor V.V.Voevodin for the opportunity to use the supercomputer SKIF MSU "Chebyshev".

This investigation was started in the framework of the project headed by Professor V.A.Tkachuk, the Dean of the Faculty of the Fundamental Medicine of the MSU. An overview of preliminary results of the work was presented at the International conference "Postgenomic methods of analysis in biology, and laboratory and clinical medicine" in the talk by Professor A.V.Bulinski (see [START_REF] Bulinski | Stochastic methods of identification of SNP interactions[END_REF] and [START_REF] Bulinski | Investigations in the framework of the MSU Research project Postgenomic Medical Studies and Technologies[END_REF]).

Methods

We start with some definitions. Let N be the the number of patients in the sample and the vector X j = X j 1 , . . . , X j n consist of genetic (SNP) and external risk factors of j-th individual, j = 1, ..., N. Here n is the total number of factors, and X j i is the value of i-th variable (characterizing SNP or external factor) of j-th individual. These variables are also called explanatory variables or predictors. If X i stands for an SNP, we set

X i =      0, no mutation in i-th SNP, 1, heterozygous mutation, 2, homozygous mutation. (1)
We assume that the external risk factors also take no more than three values, denoted by 0, 1 and 2. For example, we can specify a presence or an absence of obesity (or hypercholesterolemia etc.) by values 1 and 0 respectively. If the external factor takes more values (e.g., blood pressure), we can divide individuals into three groups according to its values.

Further on X j 1 , . . . , X j m stand for genetic data and X j m+1 , . . . , X j n for external risk factors. Let a binary variable Y j (response variable) be equal to 1 for a case, i.e. whenever j-th individual is diseased, and to -1 otherwise (that is for a control). Set ξ = (ξ 1 , . . . , ξ N ) where ξ j = (X j , Y j ), j = 1, . . . , N.

(

Suppose ξ 1 , . . . , ξ N are i.i.d. discrete random vectors having the same law as a vector (X, Y ) and independent of this vector. Assume that X = (X 1 , . . . , X n ). All random vectors (and random variables) are considered on a probability space (Ω, F , P), E denotes the integration w.r.t. P.

The main problem is to find a function in genetic and external risk factors describing the phenotype (that is the individual healthy or sick) in the best way.

Prediction algorithms

Let X := {0, 1, 2} n denote the space of all possible values of explanatory variables. Any function f : X → {-1, 1} will be called a theoretical prediction function. Define the balanced or normalized prediction error for the theoretical prediction function f as

Err(f ) := E|Y -f (X)|ψ(Y )
where the penalty function ψ : {-1, 1} → R + . Obviously,

Err(f ) = 2ψ(-1)P(f (X) = 1, Y = -1) + 2ψ(1)P(f (X) = -1, Y = 1). ( 3
)
Clearly Err(f ) depends also on the law of (X, Y ). Following [START_REF] Velez | A balanced accuracy function for epistasis modeling in imbalanced datasets using multifactor dimensionality reduction[END_REF] and [START_REF] Winham | A comparison of internal validation techniques for multifactor dimensionality reduction[END_REF] we put

ψ(y) = 1 4P(Y = y) , y ∈ {-1, 1}, (4)
the trivial cases P(Y = -1) = 0 and P(Y = 1) = 0 are excluded. Then

Err(f ) = 1 2 P(f (X) = 1|Y = -1) + 1 2 P(f (X) = -1|Y = 1). (5)
For a balanced sample considered in [START_REF] Ritchie | Multifactor-dimensionality reduction reveals high-order iteractions among estrogen-metabolism genes in sporadic breast cancer[END_REF], P(Y = -1) = P(Y = 1) = 1/2 and Err(f

) = E|Y -f (X)|/2 is equal to the classification error P(Y = f (X)).
The reason to consider this weighted scheme is that a misclassification in a more rare class should be taken into account with a greater weight. Otherwise, if the probability of disease P(Y = 1) is small, then the trivial function f (x) ≡ -1 may have the least prediction error. The approach to calculation of the prediction error based on penalty functions is not the only one possible. Nevertheless Velez et al. [START_REF] Velez | A balanced accuracy function for epistasis modeling in imbalanced datasets using multifactor dimensionality reduction[END_REF] showed that for models with high computational costs it outperforms substantially other methods such as over-and undersampling.

It is easy to prove that the optimal theoretical prediction function minimizing the balanced prediction error is given by

f * (x) = 1, p(x) > P(Y = 1), -1, otherwise. ( 6
)
where

p(x) = P(Y = 1|X = x), x ∈ X . (7)
Then each multilocus genotype (with added external risk factors)

x ∈ X is classified as high-risk if f * (x) = 1 or low-risk if f * (x) = -1.
Since p(x) and P(Y = 1) are unknown, the immediate application of ( 6) is not possible. Thus we try to find an approximation of unknown function f * using a prediction algorithm that is a function f P A = f P A (x, ξ(S)) with values in {-1, 1} (recall that Y ∈ {-1, 1} a.s.) which depends on x ∈ X and the sample ξ(S) = {ξ j , j ∈ S} where S ⊂ {1, . . . , N}.

The simplest way is to employ formula (6) with p(x) and P(Y = 1) replaced by their statistical estimates. For example introduce

p(x, ξ(S)) = j∈S I{Y j = 1, X j = x} j∈S I{X j = x} , x ∈ X , (9)
and take

P S (Y = 1) = 1 ♯S j∈S I{Y j = 1} (10)
where I{A} stands for the indicator of an event A and ♯D denotes the cardinality of a finite set D.

Along with (10) we will consider

P S (Y = 1|X ∈ C) = j∈S I{Y j = 1, X j ∈ C} j∈S I{X j ∈ C} , C ⊂ X . (11)
Thus ( 9) is a special case of [START_REF] Chawla | Data mining for imbalanced datasets: An overview. Data mining and knowledge discovery handbook[END_REF] for C = {x} with x ∈ X . Note that more difficult way is to search for the estimators of f * using several subsamples of ξ.

Assume that we constructed a prediction algorithm f P A . Then taking in mind [START_REF] Biau | Consistency of Random Forests and Other Averaging Classifiers[END_REF] set

Err(f P A (•, ξ(S))) = 1 2 y∈{-1,1} P (f P A (X, ξ(S)) = y|Y = y) . ( 12
)
As a law of (X, Y ) is unknown one can only construct an estimate Err(f P A (•, ξ(S))) of Err(f P A (•, ξ(S))). In Section 3 we use the estimated prediction error of a prediction algorithm f P A which is based on K-fold cross-validation and has the form

Err K (f P A (•, ξ), ξ) = 1 2 y∈{-1,1} 1 K K k=1 j∈S k I f P A (X j , ξ(S k )) = y, Y j = y j∈S k I{Y j = y} (13)
where

S k = (k -1) N K + 1, . . . , k N K I{k < K} + NI{k = K} , (14)
S k = {1, . . . , N} \ S k and [a] is the integer part of a ∈ R.

A very important problem is to make sure that the prediction algorithm f P A gives statistically reliable results. The quality of an algorithm is determined by its prediction error [START_REF] Coffey | An application of conditional logistic regression and multifactor dimensionality reduction for detecting gene-gene interactions on risk of myocardial infarction: the importance of model validation[END_REF] which is unknown and therefore the inference is based on consistent estimates of this error. Clearly the high quality of an algorithm means that it captures the dependence between predictors and response variables, so the error is made more rarely than it would be if these variables were independent. Consider a null hypothesis H 0 that X and Y are independent. If they are in fact dependent, then for any significant prediction algorithm f P A an appropriate test procedure involving f P A should reject H 0 at the approximate significance level α, e.g., 5%. Intuitively, this shows that results of the algorithm could not be obtained by chance. For such a procedure, we take a permutation test (see [START_REF] Golland | Permutation tests for classification[END_REF]). Its idea is as follows.

Permutation test for a given statistic L(ξ) (we consider L(ξ) = Err(f P A (•, ξ)) is done by the following steps.

1. Generate B independent random vectors (π b 1 , . . . , π b N ), 1 ≤ b ≤ B, with the uniform distribution over all permutations Π N of 1, . . . , N.

Compute Err

K,b = Err K (f P A , ξb ), 1 ≤ b ≤ B, with ξb = X 1 , Y π b 1 , . . . , X N , Y π b N .
3. Find the Monte Carlo p-value (see, e.g., [27, p. 63]):

p = F Err K (f P A , ξ) (15)
where F = F (z) is the empirical cumulative distribution function (c.d.f.) defined by the relation

F (z) = 1 B B b=1 I Err K,b ≤ z , z ∈ R. 4. If p < α, reject H 0 , otherwise not.
According to [START_REF] Golland | Permutation tests for classification[END_REF], one ideally has to use all permutations belonging to Π N but this is impractical in view of computational costs. Thus the Monte Carlo approximations for the true p-value p = F Err K (f P A , ξ) are employed, here F is the c.d.f. of Err K,b . The upper bound for |p -p| is 1/2 √ B (see [START_REF] Golland | Permutation tests for classification[END_REF]). This could be used to determine the number B of simulations for a desired accuracy.

Note also that if the estimate of the error function for the algorithm f P A (•, ξ) is asymptotically optimal, i.e. converges in probability to the error of the optimal prediction function f * as N → ∞ (ξ depends on N), then the rule of thumb is to suspect overfitting if Err K (f P A , ξ) is close to 1/2, which is a probability limit of this error under H 0 as N → ∞.

We use complementary approaches to analyze dataset related to complex diseases. Each approach (MDR, LR and machine learning) is characterized by its own way of constructing prediction algorithms. For each method one or several prediction algorithms admitting the least estimated prediction error are found (a typical situation is that there are several ones with almost the same estimated prediction error). These prediction algorithms provide a way to determine the domains where the disease risk is high or low (depending on the value of the corresponding prediction function). It is also possible to select combinations of SNPs and external risk factors whose presence influences the liability to disease to a great extent. Some methods allow to present such combinations immediately. Others, which employ more complicated forms of dependence between explanatory and response variables, need further analysis based on modifications of permutation tests. Now we pass to the description of various statistical methods and their applications to the cardiovascular risk detection.

Multifactor dimensionality reduction

Ritchie et al. [START_REF] Ritchie | Multifactor-dimensionality reduction reveals high-order iteractions among estrogen-metabolism genes in sporadic breast cancer[END_REF] introduced multifactor dimensionality reduction (MDR) as a new method of analyzing gene-gene and gene-environment interactions. Rather soon the method became very popular. Since the first publication more than 200 papers applying MDR in genetic studies were written (see, e.g., references in [START_REF] Velez | A balanced accuracy function for epistasis modeling in imbalanced datasets using multifactor dimensionality reduction[END_REF]).

MDR is a flexible non-parametric method not depending on a particular inheritance model. We give a rigorous description of the method following ideas of [START_REF] Ritchie | Multifactor-dimensionality reduction reveals high-order iteractions among estrogen-metabolism genes in sporadic breast cancer[END_REF] and [START_REF] Velez | A balanced accuracy function for epistasis modeling in imbalanced datasets using multifactor dimensionality reduction[END_REF]. As mentioned earlier, the probability p(x) introduced in ( 7) is unknown. To find its estimate one can apply maximum likelihood approach assuming that the random variable I{Y = 1} conditionally on X = x has a Bernoulli distribution with unknown parameter p(x). Then we come to [START_REF] Bulinski | Stochastic methods of identification of SNP interactions[END_REF].

A direct calculation of estimate in [START_REF] Bulinski | Stochastic methods of identification of SNP interactions[END_REF] with exhaustive search over all possible values of x is highly inefficient, since the number of different values of x grows exponentially with number of risk factors. Moreover, such a search leads to overfitting. Instead, it is usually supposed that p(x) non-trivially depends not on all, but certain variables x i . That is, there exist l ∈ N, l < n, and (k * 1 , . . . , k * l ), where 1 ≤ k * 1 < . . . < k * l ≤ n, such that for each x = (x 1 , . . . , x n ) ∈ X , the following relation holds:

p(x) = P(Y = 1|X k * 1 = x k * 1 , ..., X k * l = x k * l ). ( 16
)
In other words only few factors influence the disease and the others can be neglected. A minimal combination of factors (X k * 1 , . . . , X k * l ) in formula ( 16) is called the most significant. Clearly it is the most significant combination which has the least prediction error. Indeed, if we consider any other combination of pairwise different indices k 1 , . . . , k r and set

f k 1 ,...,kr (x) = 1, P(Y = 1|X k 1 = x k 1 , ..., X kr = x kr ) > P(Y = 1), -1, otherwise, then we obviously have Err f k * 1 ,...,k * l ≤ Err(f k 1 ,...,kr ) (17)
where Err(f ) is calculated according to [START_REF] Biau | Consistency of Random Forests and Other Averaging Classifiers[END_REF].

To choose the most significant combination, exhaustive search over all possible combinations of factors is applied. For each {k 1 , . . . , k r } ⊂ {1, . . . , n} and any x ∈ X consider

C k 1 ,...,kr (x) = {u = (u 1 , . . . , u n ) ∈ X : u k i = x k i , i = 1, . . . , r}
and for S appearing in (8) define a prediction algorithm (cf. ( 6)) by

f k 1 ,...,kr (x, ξ(S)) := 1, P S (Y = 1|X ∈ C k 1 ,...,kr (x)) > P S (Y = 1), -1, otherwise, (18)
here we use formulas [START_REF] Bulinski | Investigations in the framework of the MSU Research project Postgenomic Medical Studies and Technologies[END_REF] and [START_REF] Chawla | Data mining for imbalanced datasets: An overview. Data mining and knowledge discovery handbook[END_REF]. It is easy to show that Err K ( f k 1 ,...,kr , ξ) → Err(f k 1 ,...,kr ) in probability as N → ∞ (ξ depends on N). Consequently, (17) implies that, for any ε > 0 and all N large enough, with probability close to 1 one has

Err K ( f k *
Hence, it is natural to pick one or a few combinations of factors with the smallest empirical prediction errors as an approximation for the most significant combination. The last step in MDR is to determine statistical significance of the results. Here we test a null hypothesis of independence between X and Y i.e. between risk factors X and a disease Y . This can be done via the permutation test described in Section 2.1.

MDR method with "independent rule". We propose multifactor dimensionality reduction with "independent rule" (MDRIR) method to improve the estimate of probability p(x). This approach is motivated by Park [START_REF] Park | Independent rule in classification of multivariate binary data[END_REF], who deals with classification of large array of binary data. The principal difficulty with employment of formula ( 9) is that the number of observations in numerator and denominator of the formula might be small even for large N (see, e.g., [START_REF] Lee | Log-linear model-based multifactor dimensionality reduction method to detect gene-gene interactions[END_REF]). This can lead to inaccurate estimates and finally to a wrong prediction algorithm. Moreover, for some samples the denominator of ( 9) can equal zero.

The Bayes formula implies that

p(x) = P(X = x|Y = 1)P(Y = 1) P(X = x|Y = 1)P(Y = 1) + P(X = x|Y = -1)P(Y = -1) , (19)
where the trivial cases P(Y = -1) = 0 and P(Y = 1) = 0 are excluded. Substituting ( 19) into ( 6) we obtain the following expression for prediction function:

f * (x) = 1, P(X = x|Y = 1) > P(X = x|Y = -1), -1, otherwise. (20)
As in standard MDR method described above, we will assume that formula (16) holds. It was proved in [START_REF] Park | Independent rule in classification of multivariate binary data[END_REF] that for a broad class of models (e.g., Bahadur model [START_REF] Bahadur | A representation of the joint distribution of responses to n dichotomous items[END_REF], logit model [START_REF] Cox | The analysis of multivariate binary data[END_REF]) the conditional probability P (X k 1 = x 1 , . . . , X kr = x r |Y = y ), where y = ±1, can be estimated in the following way:

P S (X k 1 = x 1 , . . . , X kr = x r |Y = y ) = r i=1 P S (X k i = x i |Y = y ) , (21)
here (cf. ( 11))

P S (X k i = x |Y = y ) = j∈S I{X j k i = x, Y j = y} j∈S I{Y j = y} . (22)
Combining ( 16), ( 20), ( 21) and ( 22) we find the desired estimate of f * (x).

A number of observations in numerator and denominator of ( 22) increases considerably comparing with [START_REF] Golland | Permutation tests for classification[END_REF]. It allows to estimate the conditional probability more precisely whenever the estimate introduced in ( 21) is reasonable.

Thus, as opposed to standard MDR method, MDRIR uses alternative estimates of conditional probabilities. All other steps (prediction algorithm construction, prediction error calculation) remain the same. Let us mention that as far as we know this modification of MDR has not been applied before. It is based on a combination of the original MDR method [START_REF] Ritchie | Multifactor-dimensionality reduction reveals high-order iteractions among estrogen-metabolism genes in sporadic breast cancer[END_REF] and the ideas of [START_REF] Park | Independent rule in classification of multivariate binary data[END_REF].

Logic regression

The logic regression (LR) was proposed in [START_REF] Ruczinski | it Logic regression[END_REF]. Further generalizations are given in [START_REF] Fritsch | Comparing logic regression based methods for identifying SNP interactions[END_REF], [START_REF] Kooperberg | Logic regression for analysis of the association between genetic variation in the reninangiotensin system and myocardial infarction or stroke[END_REF], [START_REF] Schwender | Empirical Bayes analysis of single nucleotide polymorphisms[END_REF], [START_REF] Schwender | Identification of SNP interactions using logic regression[END_REF] and other works. LR is based on the classical binary logistic regression [START_REF] Hosmer | Applied Logistic Regression[END_REF] and exhaustive search for relevant predictor combinations. The main difficulty is to organize the search efficiently. The LR method was applied to identification of the most significant SNP combinations in [START_REF] Albrechtsen | A Bayesian multilocus association method: allowing for higher-order interaction in association studies[END_REF], [START_REF] Nunkesser | Detecting highorder interactions of single nucleotide polymorphisms using genetic programming[END_REF] and [START_REF] Schwender | Empirical Bayes analysis of single nucleotide polymorphisms[END_REF]. Note that for genetic analysis it is convenient to use explanatory variables taking 3 values. Thus we employ ternary variables, whereas the authors of the above-mentioned papers employ binary ones.

Let p(x) be the conditional probability of a disease defined in [START_REF] Bondy | Graph Theory[END_REF]. We suppose that trivial situations when p(x) ∈ {0, 1} do not occur and omit them from the consideration. To estimate p(x) we pass now to the logistic transform

q(x) = λ (p(x)) ( 23
)
where λ(z) = ln(z/(1 -z)), z ∈ (0, 1), is the inverse logistic function. The logistic function itself equals to Λ(t) = (1 + e -t ) -1 , t ∈ R. Note that we are going to estimate the unknown disease probability with the help on linear statistics with appropriately selected coefficients. Therefore it is natural to avoid restrictions on possible values of the function estimated. Thus the logistic transform is convenient, because p(x) ∈ (0, 1) for x ∈ X while q(x) can take all real values. Consider a class G of all real-valued functions in ternary variables x 1 , . . . , x n . We call a model of the dependence between the disease and explanatory variables any subclass

M ⊂ G. Set ψ(y, ξ(S)) = 1 4 P S (Y = y) , y ∈ {-1, 1},
here P S (Y = y) was introduced in [START_REF] Bulinski | Investigations in the framework of the MSU Research project Postgenomic Medical Studies and Technologies[END_REF]. Define the normalized smoothed score function

L(h, ξ(S)) = 1 ♯S j∈S φ(-Y j h(X j )) ψ(Y j , ξ(S)) ( 24
)
where S is introduced in (8), φ(t) = log 2 (1 + e t ) for t ∈ R, and h ∈ M. In contrast to previous works our version of LR scheme involves normalization (cf. ( 3)), i.e. taking the observations with weights dependent on the proportion of cases and controls in subsample ξ(S).

An easy computation yields that arg min h∈M L(h, ξ(S)) equals to arg max

h∈M 1 ♯S j∈S ln Λ(h(X j )) I{Y j = 1} 2 P S (Y = 1) + ln(1 -Λ(h(X j ))) I{Y j = -1} 2 P S (Y = -1) .
That is, minimizing the score function is equivalent to the search of normalized maximal likelihood estimate of q. Note that estimating the disease probability in this setup is closely connected with the problem of data classification, i.e. predicting the disease by the value of x ∈ X . Recall that in standard classification problem instead of the score function [START_REF] Kohara | Identification of Hypertension-Susceptibility Genes and Pathways by a Systemic Multiple Candidate Gene Approach: The Millennium Genome Project for Hypertension[END_REF] one uses the following normalized estimate of the error probability

L(h, ξ(S)) = 1 ♯S j∈S I{Y j h(X j ) < 0} ψ(Y j , ξ(S)
(x) = (x 1 × x 2 ) × (x 3 + x 4 ).
In fact the optimal choice of h for these problems coincide if the underlying model M is correctly specified (i.e. q ∈ M), see [START_REF] Biau | Consistency of Random Forests and Other Averaging Classifiers[END_REF]. However the usage of score function L has an important advantage over L because one has to evaluate the minimum of a smooth function.

A wide and easy to handle class of models is obtained by taking functions linear in variables x 1 , . . . , x n or in their products. In turn these functions admit a convenient representation by elementary polynomials. Recall that an elementary polynomial (EP) is a function T in ternary variables x 1 , . . . , x n belonging to {0, 1, 2} which can be represented as a finite sum of products x u 1 1 . . . x un n where u 1 , . . . , u n ∈ Z + . The addition and multiplication of ternary variables is considered by modulo 3. Any EP can be represented as a binary tree4 in which knots (vertices which are not leaves) contain either addition or multiplication sign, and each leaf corresponds to a variable. Figure 1 provides an example of a binary tree. Different trees may correspond to the same EP, thus this relation is not one-to-one. However, it does not influence our problem, so we regain the notation T for a tree. A finite set of trees F = (T 1 , . . . , T s ) is called a forest. For a tree T, its complexity C(T ) is the number of leaves. The complexity C(F ) of a forest F is the maximal complexity of trees constituting F.

It is clear that if g ∈ G then there exists s ≥ 1 such that g has the following form:

g(x 1 , . . . , x n ) = β 0 + s i=1 β i T i (x 1 , . . . , x n ), (25)
here β 0 , β 1 , . . . , β s ∈ R and T 1 , . . . , T s are EP.

Let us say that function g belongs to a class G r (s), where s, r ∈ N, if there exist a decomposition (25) of g such that all trees T i (i = 1, . . . , s) have complexity less or equal r. We identify a function g ∈ G r (s) with pair (F, β) where F is the corresponding forest and β = (β 0 , . . . , β s ) is the vector of coefficients in [START_REF] Kooperberg | Logic regression for analysis of the association between genetic variation in the reninangiotensin system and myocardial infarction or stroke[END_REF].

Minimization of L(h, ξ(S)) defined by ( 24) over all functions h ∈ M ⊂ G r (s) is done in two alternating steps. First we find the optimal value of β while F is fixed (which is the minimization of a smooth function in several variables) and then we search for the best F . Here one uses stochastic algorithms, since the number of such forests increase rapidly when the complexity r grows. For s ∈ N, a forest F = (T 1 , . . . , T s ) and a subsample ξ(S) (see [START_REF] Breiman | Using random forest to learn imbalanced data[END_REF]) consider a prediction algorithm f F LR setting

f F LR (x, ξ(S)) = 1, h(x) > 0, -1, otherwise,
where h = (F, β) and

β = arg min β L β 0 + s j=1 β j T j (•), ξ(S) . (26)
Define also the normalized prediction error of a forest F = (T 1 , . . . , T s ) as

ϕ(F ) = Err K (f F LR (•, ξ), ξ).
A subgraph B of a tree T is called a branch if it is itself a binary tree (i.e. it can be obtained by selecting one vertex of T together with its offspring). Sum and product signs standing in a knot of a tree are called operations, thus * stands for sum or product. Following [START_REF] Ruczinski | it Logic regression[END_REF], call the tree T a neighbor of T if it is obtained from T via one and only one of the following transformations.

1. Changing one variable to another in a leaf of the tree T (variable change).

2. Replacing an operation in a knot of a tree T with another one, i.e. sum to product or vice versa (operator change).

3. Changing a branch of two leaves to one of these leaves (deleting a leaf).

4. Changing a leaf to a branch of two leaves, one of which contains the same variable as in initial leaf (splitting a leaf).

5. Replacing a branch B 1 * B 2 with the branch B 1 (branch pruning).

6. Changing a branch B to a branch x j * B (branch growing), here x j is a variable.

Figure 2 depicts results of these operations applied to the tree T of Figure 1. We say that forests F and F are neighbors if they can be written as F = {T 1 , T 2 , . . . , T s } and F = { T 1 , T 2 , . . . , T s } where T 1 and T 1 are neighbors. The neighborhood relation defines a finite connected graph on all forests of equal size s with complexity not exceeding r. To each vertex F of this graph we assign a number ϕ(F ). To find the global minimum of a function defined on a finite graph we employ the simulated annealing method (see, e.g., [START_REF] Hajek | Cooling schedules for optimal annealing[END_REF], [START_REF] Nikolaev | Simulated Annealing. Handbook of Metaheuristics[END_REF] and [START_REF] Rubenthaler | Fast simulated annealing in R d with an application to maximum likelihood estimation in state-space models[END_REF]). This method constructs some specified Markov process which takes values in the graph vertices and converges with high probability to the global minimum of the function. To avoid stalling at a local minimal point the process is allowed to pass with some small probability to a point F having greater value of ϕ(F ) than current one. We propose a new modification of this method in which the output is the forest corresponding to the minimal value of a function ϕ(F ) over all (randomly) visited points. x 1 = 2

x 3 = 2 x 4 > 0 1(yes) 0(no) 1(yes) 0(no) 1(yes) 0(no) Y = 1 Y = -1 Y = 1 Y = -1

Machine learning methods

Let us describe (see, e.g., [START_REF] Szymczak | Machine Learning in Genome-Wide Association Studies[END_REF]) two among the most popular machine learning methods random forests (RF) and stochastic gradient boosting (SGB). They belong to ensemble methods which combine multiple predictions from a certain base algorithm to obtain better predictive power (i.e. less estimated prediction error). We use classification and regression trees (CART) for a base learning algorithm because it showed good performance in a number of studies (see [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]).

Classification tree T is a binary tree having the following structure. Any leaf of T contains either 1 or -1 and for any vertex P in T (including leaves) there exists a subset A P of the explanatory variable space X such that the following properties hold: 1. A P = X if P is the root of T ; 2. if vertices P ′ and P ′′ are children for P , then A P ′ ∪ A P ′′ = A P and A P ′ ∩ A P ′′ = ∅.

In particular, subsets corresponding to the leaves form the partition of X . To obtain a prediction of Y given a certain value x ∈ X of the random vector X, one should go along the path which starts from the root and ends in some leaf turning at each parent vertex P to that child P ′ for which A P ′ contains x. At the end of the x-specific path, one gets either 1 or -1 which serves as a prediction of Y . Figure 3 provides an example of a classification tree. Namely, the partition of X is formed by values of boolean functions standing in parent vertices. For each x starting from the root of the tree we calculate the value of a boolean function and move along the edge marked with the value obtained (1 or 0). The left child of the root corresponds to the subset {x ∈ X : x 1 = 2}, while the right one to its complement in X . Next, the leftmost leaf stands for a subset {x ∈ X : x 1 = 2, x 3 = 2}, and if X falls in this subset, we predict that Y = 1; the rightmost leaf stands for a subset {x ∈ X : x 1 < 2, x 4 = 0}, and if X takes values in this subset, we predict Y = -1.

Classification tree could be constructed via CART algorithm (if it is the case, we will call it CART). The algorithm proceeds iteratively. That is, on the l-th step of the algorithm (l = 1, 2, . . .), each element A of the current partition A l (A 1 = X ) of the set X is divided into two disjoint parts

A + (i, t) = {(x 1 , . . . , x n ) ∈ A : x i ≤ t} and A -(i, t) = {(x 1 , . . . , x n ) ∈ A : x i > t} minimizing the sum G(A + (i, t)) + G(A -(i, t)) over i = 1, . . . , n and t ∈ {0, 1}. Here the empirical Gini index G(C) = 2 P S (Y = 1|X ∈ C) 1 -P S (Y = 1|X ∈ C)
with C ⊂ X and P S (Y = 1|X ∈ C) (see [START_REF] Chawla | Data mining for imbalanced datasets: An overview. Data mining and knowledge discovery handbook[END_REF]) measures the heterogeneity of the subsample {j ∈ S : X j ∈ C} w.r.t. response variable Y . Any uninformative partition with min

(i,t) G(A + (i, t)) + G(A -(i, t)) > G(A),
is not made.

The algorithm stops whenever a number of leaves D reaches some critical level which is chosen via some data-dependent criteria (see [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF], page 308). For a subsample ξ(S) of ξ, each CART defines a prediction algorithm

f (x, ξ(S)) = D d=1 a d (ξ(S))I{x ∈ A d (ξ(S))} (27)
where {A 1 (ξ(S)), . . . , A D (ξ(S))} is the partition of X corresponding to the leaves,

a d (ξ(S)) = 1, ♯{j∈ S : Y j = 1, X j ∈ A d (ξ(S))} > ♯{j∈ S : Y j = -1, X j ∈ A d (ξ(S))}; -1, otherwise.
RF is a non-parametric method of estimating conditional probability p(x). It was successfully applied to genetics data in a number of papers (see references in [START_REF] Szymczak | Machine Learning in Genome-Wide Association Studies[END_REF]). It could be briefly described as follows (see chapter 15 in [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] for details). Generate B bootstrap samples from the initial sample where one could choose B = max{[N log N], 1000} according to [START_REF] Olive | The number of samples for resampling algorithms[END_REF]. For b-th bootstrap sample (1 ≤ b ≤ B) construct a CART prediction algorithm f b : X × (X × {-1, 1}) N → {-1, 1} defined according to [START_REF] Lehmann | Testing Statistical Hypotheses[END_REF] and take

p RF (x, ξ(S)) = B -1 B b=1 f b (x, ξ(S)) + 1 /2
as an estimate of p(x).

It is shown in [START_REF] Biau | Consistency of Random Forests and Other Averaging Classifiers[END_REF] that generally RF method gives consistent estimates of p(x) only if the number of partitions used in CART grows slower than the sample size. A final prediction algorithm f RF (x, ξ(S)) is constructed from the estimate p RF (x, ξ(S)) similarly to [START_REF] Bock | Computational epigenetics[END_REF], i.e.

f RF (x, ξ(S)) = 1, p RF (x, ξ(S)) > P S (Y = 1), -1, otherwise.
The distinctive features of this method are low computational costs and the ability to extract relevant predictors when the number of irrelevant ones is large (see [START_REF] Biau | Analysis of a Random Forests model[END_REF]).

SGB is another non-parametric method of estimating conditional probability p(x). This method is used in a number of procedures for studying genetics data (see, e.g., [START_REF] Wan | MegaSNPHunter: a learning approach to detect disease predisposition SNPs and high level interactions in genome wide association study[END_REF]). SGB method can be described as follows ([15]).

f m = f m-1 + ρ D d=1 w m d (ξ m (S))I A m d (ξm(S)) ,
here ρ is the memory relaxation parameter.

4. If m < M, go to Step 3, otherwise determine a final estimate

p SGB (x, ξ(S)) = 1 1 + exp{-2f M (x, ξ(S), {ξ m (S)} M m=1 )} .
This algorithm is to be run for several times with different parameters D, M, ρ and η. Then their optimal values could be chosen via cross-validation (see section 16.3.1 in [21]). Small values of η (= 0.1, 0.05, 0.0225 etc.) help to get accurate estimates for relatively noisy data.

Standard RF and SGB work poorly for unbalanced samples. One needs either to balance given datasets (as in [START_REF] Chawla | Data mining for imbalanced datasets: An overview. Data mining and knowledge discovery handbook[END_REF]) before these methods are applied or use special modifications of RF ( [START_REF] Breiman | Using random forest to learn imbalanced data[END_REF]) and SGB ( [START_REF] Liu | Exploratory Undersampling for Class-Imbalance Learning[END_REF]). To avoid overfitting, permutation test needs to be done.

A common problem of all machine learning methods is a complicated functional form of the final probability estimate p(x, ξ) (w.r.t. x). In genetic studies, one wants to pick up all relevant combinations of SNP and risk factors, based on a biological pathway causing the disease. Therefore, the final estimate p(x, ξ) is to be analyzed. We describe one of the possible methods of such analysis within RF framework and called conditional variable importance measure (CVIM). One could determine CVIM for each predictor X i in X and range all X i in terms of this measure. Following [START_REF] Strobl | Conditional variable importance for random forests[END_REF], CVIM of predictor X i given certain subvector Z i of X is calculated as follows (supposing Z i takes values z i1 , . . . , z im(i) ).

1. Construct a vector (l 1 , . . . , l N ), randomly permuting 1, . . . , N in each subset

A ik = {j : Z j i = z ik }, k = 1, . . . , m(i).
CVIM b = 1 |C b | j∈C b I{Y j = f b (X j , ξ b )} - 1 |C b | j∈C b I{Y j = f b (X l j , ξ b )} where C b = {j ∈ {1, . . . , N} : (X j , Y j ) / ∈ ξ b }.
3. Compute the final CVIM using the formula

CV IM = B -1 B b=1 CVIM b . (28)
Any permutation (l 1 , . . . , l N ) in the CVIM algorithm destroys dependence between X i and (Y, Z -i ) where Z -i consists of all components of X which are not in Z i . At the same time it preserves initial empirical distribution of (X i , Z i ) calculated for the sample ξ. After that the average loss of correctly classified Y is calculated. If it is relatively large w.r.t. CVIM of other predictors, then X i plays important role in classification and vice versa.

For Z i , one could take all components X k (k = i) such that the hypothesis of the independence between X k and X i is not rejected at some significance level (e.g., 5%). Note also that CVIM-like algorithm could be used to range pairs of SNP and risk factors w.r.t. the level of association to the disease. This will be done elsewhere.

Applications: risks of IHD and MI

We employ here the various statistical methods described above to analyze the influence of genetic and external factors on risks of ischemic heart disease (IHD) and myocardial infarction (MI) using the data for 454 individuals (333 cases, 121 controls) and 333 individuals (165 cases, 168 controls) respectively. These data contain values of seven SNPs (PAI-1, GpIa, GpIIIa, FXIII, FVII, IL-6, Cx37), as well as four external risk factors, namely, obesity (Ob), arterial hypertension (AH), smoking (Sm) and hypercholesterolemia (HC). The age of all individuals in case and control groups ranges from 35 to 55 years, which reduces its influence on the risk analysis. For each of considered methods, K-fold crossvalidation is used with K = 6. As shown in [START_REF] Motsinger | The effect of reduction in cross-validation intervals on the performance of multifactor dimensionality reduction[END_REF] and [START_REF] Winham | A comparison of internal validation techniques for multifactor dimensionality reduction[END_REF], the standard choice of partition number of cross-validation from 6 to 10 does not change the prediction error significantly. We take K = 6 as the sample sizes do not exceed 500. The supercomputer SKIF MSU "Chebyshev" was involved to perform computations. All applied methods have prediction error less than 0.25, so predictions constructed have significant predictive power. Indeed, in [START_REF] Coffey | An application of conditional logistic regression and multifactor dimensionality reduction for detecting gene-gene interactions on risk of myocardial infarction: the importance of model validation[END_REF] and [START_REF] Tsai | Reninangiotensin system gene polymorphisms and coronary artery disease in a large angiographic cohort: Detection of high order gene-gene interaction[END_REF] the interplay between genotype characteristics and MI development was also studied, with estimated prediction errors 0.30-0.40. Further on we write prediction error instead of estimated prediction error.

MDR and MDRIR methods

Ischemic heart disease

Table 1 contains (estimated) prediction errors of the most significant combinations obtained by MDR analysis of ischemic heart disease data. At Figure 4 a plot of empirical distribution function of prediction error is given when the disease is not linked with explanatory variables. We use here the simulated samples ξ b introduced in Section 2.1, with b = 1, . . . , B where B = 100. One can see that out of these 100 simulations, the corresponding prediction error was not less than 0.42. Note that Monte Carlo p-value (15) of all three combinations is less than 0.01 (since their prediction errors are much less than 0.42), which is usually considered as a good performance.

Factors

Prediction error GpIa, FXIII, AH, HC 0.231 Cx37, AH, HC 0.238 GpIa,Cx37, AH, HC 0.241 Table 1: The most significant combinations obtained by MDR analysis for IHD data.

Table 2 contains the results of MDRIR method, which are similar to results of MDR method. However, it is worth mentioning that MDRIR method allows to identify additional combinations with prediction error around 0.24.

Factors

Prediction error FXIII, FVII, AH, HC 0.240 FXIII, AH, HC 0.242 GpIa, Cx37, AH, HC 0.247 Table 2: The most significant combinations obtained by MDRIR analysis of IHD data.

It follows from Tables 1 and2 that hypertension and hypercholesterolemia are the most important external risk factors. Indeed, these two factors appear in every of 6 combinations. To perform a more precise analysis of influence of SNPs on IHD provoking we analyze gene-gene interactions. We used two different strategies. Namely, we applied MDR method to a subgroup of individuals who are not subject to any of the external risk factors (i.e. to non-smokers without obesity and without hypercholesterolemia, 51 cases and 97 controls). Another strategy is to apply MDR method to the whole sample, but to take into account only genetic factors rather than all factors. It turned out that both methods yield similar results. Combination of SNPs GpIa and Cx37 has the biggest influence on IHD. Prediction error is about 0.28-0.34, and smaller error corresponds to a risk-free sample. Moreover it follows from Tables 2 and3 that prediction error significantly dropped after additional exogenous factors were taken into account (the error is 0.247 if additional external factors are taken into account and 0.343 if not).

Thus based on ischemic heart disease data with the help of Tables 1-3 we can make the following conclusions. Combination of two SNPs (GpIa and Cx37) and two external factors (hypertension and hypercholesterolemia) has the biggest influence on IHD. Also FXIII gives additional predictive power if AH and HC are taken into account.

Myocardial infarction

Prediction errors of the most significant combinations obtained by MDR analysis of MI data are presented in Table 4. Figure 5 contains the plot of empirical c.d.f. of prediction error if disease is not linked with risk factors. This curve shows that for all 100 simulations of ξ b the estimated prediction error was not less that 0.38. Note that Monte Carlo p-value of all combinations is less than 0.01.

Factors

Prediction error GpIIIa, FXIII, Cx37, AH 0.343 GpIIIa, FXIII, FVII, Cx37 0.347 Cx37, Sm 0.356 MDRIR analysis of the same dataset gives a clearer picture (see Table 5). Apparently, combination of smoking and SNP Cx37 is the most significant. These two factors appear in all combinations in Table 5. Involving any additional factors only increases prediction error.

The explicit form of the prediction algorithm based on Cx37 and Sm shows that these factors interact nonlinearly. Smoking as well as Cx37 homozygote leads to the disease. However wild-type allele can protect from consequences of smoking, because combination of smoking and Cx37 wild-type is a protective one (i.e. value of prediction algorithm of this combination is -1).

Logic regression

We performed several research procedures both for IHD and MI data, with different restrictions imposed on the statistical model. To describe these models set

(X 1 , . . . , X n ) = (Z 1 , . . . , Z m , R 1 , . . . , R k )
where variables Z = (Z 1 , . . . , Z m ) stand for SNP values (PAI-1, GpIa, GpIIIa, FXIII, FVII, IL-6, Cx37 respectively) and R = (R 1 , . . . , R k ) denote external risk factors (Ob, AH, Sm, HC), m = 7, k = 4. We consider four different models in order to analyze both total influence of genetic and external factors and losses in predictive force appearing when some factors are excluded. In our applications we will take s = 3, as search over larger forests for samples with modest sizes can give very complicated and unreliable results.

where β i ∈ R and T i are polynomials identified with trees. Thus we allow the interaction of genes and external factors in order to find significant gene-environment interactions. However we impose additional restrictions to avoid too complex combinations of external risk factors. We do not tackle here effects of interactions where several external factors are involved. Namely, we consider only the trees satisfying the following two conditions.

1. If there is a leaf containing external factor variable then the root of that leaf contains product operator.

2. Moreover, another branch growing from the same root is also a leaf and contains a genetic (SNP) variable.

Results obtained by RF and SGB methods

The given datasets were unbalanced w.r.t. response variable and we first applied the resampling technique to them. That is, enlargement of the smaller of two groups case-control in the sample by additional bootstrap observations till the final proportion case:control would be 1:1. We also employed modifications of RF by [START_REF] Breiman | Using random forest to learn imbalanced data[END_REF] and SGB by [START_REF] Liu | Exploratory Undersampling for Class-Imbalance Learning[END_REF] for unbalanced samples, but those worked poorly for permutation tests and we do not give their results here. Note that due to the resampling techniques the following effect arise. Some observations in small groups (case or control) appear in the new sample more frequently than other ones. Therefore, we took the average over 1000 iterations. Results of RF and SGB methods are given in Table 10. It shows that RF and SGB methods give statistically reliable results (prediction error in the permutation test is close to 50%). Moreover, additional SNP information improves predicting ability on 11% and 13% (SGB). It seems that SGB method is better fitted to IHD data than RF.

Ischemic heart disease

Data

Computing CVIM for each X i , we constructed Z i as follows. We included in Z i all predictors X j , j = i, for which χ 2 -criteria rejected hypothesis of independence between X j and X i at 5% significance level. Since the genetic information has second order effect on prediction of Y comparing to the risk factors, we ran the program 1000 times and then took the average CVIM to get a reliable estimate. An error over different runs of the program was around 0.01. The results are given in Table 11.

AH HC Cx37 Ob FXIII Sm GpIa FVII PAI-1 GpIIIa IL-6 8. [START_REF] Bulinski | Stochastic methods of identification of SNP interactions[END_REF] Thus, the most relevant predictors for IHD are AH, HC and Cx37.

Myocardial infarction

Results of RF and SGB methods are given in the following table. CVIM was calculated according to [START_REF]Bioinformatics From Genomes to Therapies[END_REF] and is given below. Thus, the most relevant predictors for MI are Cx37, Sm and AH.

Conclusions and final remarks

Let us briefly summarize the main results obtained. The analysis of IHD dataset showed that two external risk factors out of four considered (AH and HC) have a strong connection with the disease risk (the error of classification based on external factors only is 0.25-0.26 with p-value less than 0.01). Also, the classification based on SNPs only gives a relatively low error of 0.28. Moreover, the most influential SNPs are Cx37 and GpIa (FXIII also enters the analysis only when AH and HC are present). Prediction error decreases to 0.13 if both SNP information and external risk factors are taken into account. Note that excluding any of the 5 remaining SNPs (all except for two most influential) from data increases the error by 0.01-0.02 approximately. So, while the most influential data are responsible for the situation within a large part of population, there are smaller parts where other SNPs come to effect and provide a more efficient prognosis ("small subgroups effect"). The MI dataset gave the following results. The most significant factors of MI risk are the Cx37 SNP (more precisely, homozygous mutation) and smoking with a considerable gene-environment interaction present. The smallest prediction error of methods applied was 0.33-0.35 (with p-value less than 0.01). The classification based on external factors only yields a much greater error of 0.42. Thus genetic data improves the prognosis quality essentially. While two factors are of great importance, other SNPs considered actually do not improve the prognosis essentially, i.e. no small groups effect is observed.

The conclusions given above are based on several complementary methods of modern statistical analysis. These new data mining methods allow to analyze other datasets as well. The study can be continued with larger datasets, in particular, involving new SNP data.

Figure 1 :

 1 Figure 1: A tree T representing a function T (x) = (x 1 × x 2 ) × (x 3 + x 4 ).

Figure 2 :

 2 Figure 2: Neighbors of the tree T.

Figure 3 :

 3 Figure 3: CART representing the prediction ((x 1 = 2) and (x 3 = 2)) or (x 4 > 0).

1 .

 1 Pass on the input of the algorithm 5 initial parameters D, M ∈ N and ρ, η ∈ (0, 1).2. Put m = 0, ξ 0 (S) = ξ(S) and f 0 (x, ξ 0 (S))

3. Increase m by 1

 1 j f m-1 (X j , {ξ l (S)} m-1 l=0 )} . Choose a random subset in ξ m (S) = {(X j , Ȳ j m )} j∈S with [η♯S] elements. Construct CART prediction algorithm (with D leaves) D d=1 a m d (ξ m (S))I{x ∈ A m d (ξ m (S))} on the chosen subset. Compute weight coefficients w m d (ξ m (S)) = j∈J Ȳ j m j∈J | Ȳ j m |(2 -| Ȳ j m |) , d = 1, . . . , D,where the random set J = {j : X j ∈ A m d (ξ m (S))}, and put

2 .

 2 Generate B bootstrap samples ξ b = (X jb , Y jb ), j = 1, . . . , N , b = 1, . . . , B. For each of these samples, construct a classifier f b (x, ξ b ) and calculate

Figure 4 :

 4 Figure 4: The empirical distribution function of the estimated error for the case when disease risk and predictors are independent (permutation test), IHD dataset.

Figure 5 :

 5 Figure 5: The empirical c.d.f. of the estimated error for the case when disease risk and predictors are independent (permutation test), MI dataset.

Model 1 .

 1 We consider the class M (see Section 2.3) consisting of the functions h having a formh(Z, R) = β 0 + s v=1 β v T v (Z 1 , . . . , Z m ) + k v=1 β s+v R vwhere the coefficients β i ∈ R and T i are polynomials identified with trees. In other words we require that external factors are present only in trees consisting of one variable.Model 2. Now we assume that any function h ∈ M has the representationh(Z, R) = β 0 + s v=1 β v T v (Z 1 , . . . , Z m , R 1 , . . . , R k )

  ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	× . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	× . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	x 1	x 2	x 3	x 4

. . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	× . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	+ . . . . . . . . . . . + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .		× . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	x 1	x 5	x 3	x 4		x 1		x 2	x 3	x 4
		Variable change				Operator change
	x 1	+ . . . . . . . . . . . . x 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Deleting a leaf × x 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .		x 1	× . . . x 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + x 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
									x 6	x 7
						Splitting a leaf
	× . . . . . . . . . . . . x 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Branch pruning x 1	x 1	× . . . x 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × x 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
									x 3	x 4
						Branch growing

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

Table 3

 3 contains the most significant combinations of SNPs and their prediction errors.

	Method	Genetic factors Prediction
			Error
	MDR on a subgroup of individuals
	who do not have any risk factors	GpIa, Cx37	0.281
	MDR method on the whole group
	taking into account only genetic factors	GpIa, Cx37	0.343

Table 3 :

 3 Comparison of the most significant SNP combinations obtained by two different ways of MDR analysis of IHD data.

Table 4 :

 4 The most significant combinations obtained by MDR analysis of MI dataset.

Table 5 :

 5 The most significant combinations obtained by MDRIR analysis of MI dataset.

	Factors	Prediction error
	Cx37, Sm	0.351
	GpIIIa, Cx37, Sm	0.353
	GpIIIa, Cx37, Sm, HC	0.355

Table 10 :

 10 Prediction error/prediction error in permutation test calculated via cross-validation for IHD dataset with employment of RF and SGB methods.

		RF	SGB
	with SNP	0.20/0.454 0.134/0.473
	without SNP 0.23/0.51 0.261/0.503

Table 11 :

 11 Predictors are ranged in terms of their CVIM for IHD dataset.

	5.3	5.1 0.56 0.53 0.11 0.1	0.07	0.03	0.02	0.01

Table 12 :

 12 Table 12 shows that RF and SGB methods give statistically reliable estimates (prediction error in the permutation test is close to 50%). Moreover, additional SNP information improves predicting ability on 10%. Prediction error/prediction error in permutation test calculated via cross-validation for MI dataset with employment of RF and SGB methods.

		Data			RF		SGB
		with SNP data	0.36/0.497 0.399/0.53
		without SNP data 0.473/0.527 0.482/0.562
	Cx37 Sm AH GpIIIa FVII FXIII HC GpIa Ob IL-6 PAI-1
	7.5	2 1.86	0.03	0.02	≈0	≈0	≈0 ≈0 ≈0	≈0

Table 13 :

 13 Predictors are ranged in terms of their CVIM for MI dataset.

,...,k * l , ξ) < Err K ( f k 1 ,...,kr , ξ) + ε.

For the basic concepts of the graph theory see, e.g.,[START_REF] Bondy | Graph Theory[END_REF].

This algorithm is not to be confused with prediction algorithms.

Acknowledgments

The authors are grateful to Professor V.A.Tkachuk, Associate Professor L.M.Samokhodskaya and MD. A.V.Balatsky for providing the data concerning the complex diseases. On the basis of this preprint the joint paper will be prepared with medical interpretation of the obtained results.

A.V.Bulinski also would like to thank Professors I.Kourkova and G.Pagès for invitation to LPMA of the University Pierre and Marie Curie, he is grateful to all the members of LPMA for hospitality.

The work is partially supported by RFBR grant 10-01-00397.

Models 3 and 4 have additional restrictions that polynomials T v (v = 1, . . . , s) in [START_REF] Liu | Exploratory Undersampling for Class-Imbalance Learning[END_REF] depend only on external factors and only on SNPs respectively. These models are considered to compare their results with ones obtained with all information taken into account, in order to demonstrate the importance of genetic (resp. external) data for risk analysis.

Ischemic heart disease

We have the following results. Note that prediction error in Model 1 is only about 0.19. For the same model we performed also fast simulated annealing search of the optimal forest which is much more time-efficient, and a reasonable error of 0.23 was obtained. Model 3 application shows that external factors play an important role in IHD genesis, as classification based on external factors only gives the error less than 0.23, while usage of SNPs only (Model 4) lets the error grow to 0.34.

Model 1 gave the minimal prediction error. For the optimal forest (T 1 , . . . , R 4 ) the function h(Z, R) given before formula [START_REF] Lee | Log-linear model-based multifactor dimensionality reduction method to detect gene-gene interactions[END_REF] with S = {1, . . . , N} is provided by the expression -0.597T 1 -0.354T 2 + 0.521T 3 -0.444R 1 + 1.311R 2 -0.146R 3 + 2.331R 4 -0.226 [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF] where 6 T 1 = (

The external factors 2 and 4 (i.e. AH and HC) are the most influential since the coefficients at them are the greatest ones (1.311 and 2.331). As was shown above, MDR yields the same conclusion. If the gene-environment interactions are allowed (Model 2), no considerable increase in predictive force has been detected. However we list the pairs of SNPs and external factors present in the best forest: Z 7 and R 2 , Z 7 and R 1 , Z 7 and R 4 , Z 5 and R 1 . It is seen that Cx37 SNP is of substantial importance as it appears in combination with all risk factors except for smoking.

As formula [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF] is hard to interpret, we select the most significant SNPs via a variant of permutation test. Consider a random rearrangement of the column with first SNP in IHD dataset. Calculate the prediction error using these new simulated data and the same function h as before. The analogous procedure is done for other columns (containing the 6 The sums and products are modulo 3.

21

values of other SNPs) and the errors found are given in Table 7. It is seen that the error increases considerably when the values of GpIa and Cx37 are permutated. The statement that they are the main sources of risk agrees with what was obtained above by MDR method.

Prediciton error for

Myocardial infarction

For the MI dataset, under the same notations that above, the following results for our four models were obtained. To comment the Table 8 we should first underline that external risk factors play less important role compared with IHD risk: if they are used without genetic information, the error increases by 0.09, see Models 1 and 3 (while the same increase for IHD was 0.03). The function h(Z, R) defined before formula [START_REF] Lee | Log-linear model-based multifactor dimensionality reduction method to detect gene-gene interactions[END_REF] with S = {1, . . . , N} is is equal to

where

Thus the first tree has the greatest weight (coefficient equals -1.144), the second tree (i.e. Cx37 SNP) is on the second place, and external factors are less important.

As for IHD we performed a permutation test to compare the significance of different SNPs. Its results are presented in Table 9.

Prediciton error for Model 1

Cx37

GpIIIa IL-6 FXIII FVII PAI-1 GpIa 0.30526 0.44420 0.35345 0.33998 0.32761 0.32427 0.31918 0.30526