ABSTRACT

We present a methodology for constructing a short-term event risk score in heart failure patients from an ensemble predictor, using bootstrap samples, two different classification rules, logistic regression and linear discriminant analysis for mixed data, continuous or categorical, and random selection of explanatory variables to build individual predictors. We define a measure of the importance of each variable in the score and an event risk measure by an odds-ratio. Moreover, we establish a property of linear discriminant analysis for mixed data. This methodology is applied to EPHESUS trial patients on whom biological, clinical and medical history variables were measured.

Keywords:

Ensemble Predictor, Linear Discriminant Analysis, Logistic Regression, Mixed Data, Scoring, Supervised Classification

1. Introduction

In this study, we focus on the problem of constructing a short-term event risk score in heart failure patients based on observations of biological, clinical and medical history variables.

Numerous event risk scores in heart failure patients have been proposed in recent years, but one aspect is particularly important to consider in the construction of a score and in the relevance of the results obtained. This concerns the choice of classification models whose conditions of use may be restrictive. The most currently used classification models in these studies are logistic regression and Cox proportional hazard model. Quoting for example the Seattle Heart Failure Model (SHFM) risk score [1] and the Seattle Post Myocardial Infarction Model (SPIM) risk score [2] which allow respectively predicting survival in chronic and post-infarction heart failure patients:

• SHFM risk score was derived in a cohort of 1153 patients with ejection fraction < 30% and New York Heart Association (NYHA) class III to IV and validated in 5 other cohorts of patients with similar characteristics. Area under ROC curve (AUC) at 1 year was 0.725 in resubstitution and ranged from 0.679 to 0.810 in the 5 validation cohorts.

• SPIM risk score was derived in a cohort of 6632 patients from the Eplerenone Post-Acute Myocardial Infarction Heart Failure Efficacy and Survival Study (EPHESUS) trial [3] and validated on a cohort of 5477 patients. AUC at 1 year was 0.742 in derivation and 0.774 in validation.

These two risk scores were developed using Cox proportional hazard model and characteristics available at baseline as explanatory variables. Overall, there are several limitations to using these risk scores. They were constructed using only data available at baseline. However, as many studies include inclusion criteria based on clinical or biological parameters measured at baseline, it is possible that some variables are not present in the score due to these inclusion criteria. For example, patients were included in the EPHESUS trial only if their potassium level at baseline was less than 5 mmol/L. This is a reason why potassium is not present in the SPIM score although this is an important parameter which moreover may evolve considerably over time. Concerning the model, the Cox proportional risk model assumes the proportionality of risks, an important condition not always obtained and verified.

In this study, we used a new approach:

• we develop a methodology for constructing a short-term event (death or hospitalization) risk score, taking into account the most recent values of the parameters and therefore the closest values of an event, in order to generate alerts and eventually immediately modify drug prescription; using EPHESUS trial data, we could only construct a score at 1 month in order not to have too few patients with event in the learning sample; but with the same methodology, a score could be constructed at a closer time;

• we use an ensemble predictor, that is more stable than a predictor built on a single learning sample, using bootstrap samples; this allows an internal validation of the score using AUC out-of-bag (OOB); moreover, we use two classification methods, logistic regression and linear discrimination analysis, and, in order to avoid overlearning, for each predictor we use a random selection of explanatory variables, after testing other methods of selection that did not give better results, the number of drawn variables being optimized after testing all possible choices;

• furthermore, our method of construction can be adapted to data streams: when patient data arrives continuously, the coefficients of variables in the score function can be updated online.

In the next section, we present how we defined the learning sample using the available data from EPHESUS trial and the list of explanatory variables used. In the third section, we state a property of linear discriminant analysis (LDA) for mixed data, continuous or categorical. In the fourth section, after presenting the methodology used to build a risk score and to reduce its variation scale from 0 to 100, we define a measure of the importance of variables or groups of correlated variables in the score and a measure of the event risk by an odds-ratio. In the fifth section, we describe the results obtained by applying our methodology to our data. The paper ends with a conclusion.

2. Data

The database at our disposal was EPHESUS, a clinical trial that included 6632 patients with heart failure (HF) after acute myocardial infarction (MI) complicated by left ventricular systolic dysfunction (left ventricular ejection fraction < 40%) [3] . All patients were randomly assigned to treatment with eplerenone 25 mg/day or placebo.

In this trial, each patient was regularly monitored, with visits at the inclusion in the study (baseline), 1 month after inclusion, 3 months later, then every 3 months until the end of follow-up. At each visit, biological, clinical parameters or medical history were observed. In addition, all adverse events (deaths, hospitalizations, diseases) that occurred during follow-up were collected.

To define the learning sample used to construct the short-term event risk score, we made the following working hypothesis: based on biological, clinical measurements or medical history on a patient at a fixed time, we sought to assess the risk that this patient has a short-term HF event. The individuals considered are couples (patient-month) without taking into account the link between several couples (patient-month) concerning the same patient. Therefore, it was assumed that the short-term future of a patient depends only on his current measures.

Firstly, we did a full review of the database in order to:

• identify the biological and clinical variables that were regularly measured at each visit,

• determine the medical history data that we could update from information collected during the follow-up.

We were thus able to define a set of 27 explanatory variables whose list is presented in Figure 1. Estimated plasma volume derived from Strauss formula (ePVS) was defined in [4] . Estimated glomerular filtration rate (eGFR) was assessed using three formulas [5] [6] [7] . The different types of hospitalization were defined in supplementary material of [3] .

Then, we defined the response variable as the occurrence of a composite short-term HF event (death or hospitalization for progression of HF). In order to have enough events, we defined the short term as being equal to 30 days. Patient-months with a follow-up of less than 30 days and no short-term HF event during this incomplete follow-up period, were not taken into account.

There were finally 21,382 patient-months from 5937 different patients whose 317 with short-term HF event and 21,065 with no short-term event.

3. Property of Linear Discriminant Analysis of Mixed Data

Denote A' the transposed of a matrix A.

In case of mixed data, categorical and continuous, a classical method to perform a discriminant analysis is:

1) perform a preliminary factorial analysis according to the nature of the data, such as multiple correspondence factorial analysis (MCFA) [8] for categorical data, multiple factorial analysis (MFA) [9] for groups of variables, mixed data factorial analysis (MDFA) [10] , ... ;

2) after defining a convenient distance, perform a discriminant analysis from the set of values of principal components, or factors.

See for example the DISQUAL (DIScrimination on QUALitative variables) method of Saporta [11] , which performs MCFA, then LDA or quadratic discriminant analysis (QDA).

Denote as usual T the total inertia matrix of a dataset partitioned in classes, W and B respectively its intraclass and interclass inertia matrix.

We show hereafter that when performing LDA with metrics T⁻¹ or W⁻¹, it is not necessary to perform a preliminary factorial analysis and LDA can be directly performed from the raw mixed data.

Metrics W⁻¹ will be used in the following but can be replaced by T⁻¹.

Let $I=\left\{1,2,\cdots ,n\right\}$ a set of n individuals, partitioned in q disjoint classes $I_1\mathrm{,}\cdots \mathrm{,}{I}_q$ . Denote $n_k=card\left(I_k\right)$ , $p_{ki}$ the weight of i^th individual of class $I_k$ $\left(i=1,\cdots ,n_k;k=1,\cdots ,q\right)$ and $P_k=\stackrel{n_k}{\sum _{i=1}}{p}_{ki}$ the weight of $I_k$ , with $\underset{k=1}{\sum ^q}{P}_k=1$ . p quantitative variables or indicators of modalities of categorical variables, denoted $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ , are observed on these individuals. Suppose that there exists no affine relation between these variables, especially for each categorical variable an indicator is removed.

For $j=1,\cdots ,p$ , denote $x_{ki}^j$ the value of $x^j$ for i^th individual of class $I_k$ . Denote $x_{ki}$ the vector ${\left(x_{ki}^1\cdots x_{ki}^p\right)}^{\prime }$ and $g_k$ the barycenter of the elements $x_{ki}$ for $i\in I_k$ :

$g_k=\frac{1}{P_k}\sum _{i\in I_k}{p}_{ki}{x}_{ki}.$ (1)

Intraclass inertia $\left(p\mathrm{,}p\right)$ matrix W is supposed invertible:

$W=\sum _{k=1}^q\sum _{i=1}^{n_k}{p}_{ki}\left(x_{ki}-g_k\right){\left(x_{ki}-g_k\right)}^{\prime }.$ (2)

A currently used distance in LDA $d_{W^{-1}}\left(a\mathrm{,}b\right)$ between two points a and b in ${ℝ}^p$ is such that:

$d_{W^{-1}}^2\left(a,b\right)={\left(a-b\right)}^{\prime }{W}^{-1}\left(a-b\right).$ (3)

Suppose we want to classify an individual knowing the vector a of values of $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ . Principle of LDA is to classify it in $I_k$ such that $d_{W^{-1}}^2\left(a\mathrm{,}{g}_k\right)$ is minimal.

Consider now new variables $y^1\mathrm{,}\cdots \mathrm{,}{y}^m$ affine combinations of $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ , with $m\ge p$ , such that:

$y_{ki}=A{x}_{ki}+\beta ,$ (4)

with $y_{ki}={\left(y_{ki}^1\cdots y_{ki}^m\right)}^{\prime }$ , A a $\left(m\mathrm{,}p\right)$ matrix of rank p and $\beta$ a vector in ${ℝ}^m$ .

Denote $h_k$ the barycenter of vectors $y_{ki}$ in ${ℝ}^m$ for $i\in I_k$ :

$h_k=\frac{1}{P_k}\sum _{i\in I_k}{p}_{ki}{y}_{ki}=\frac{1}{P_k}\sum _{i\in I_k}{p}_{ki}\left(A{x}_{ki}+\beta \right)=A{g}_k+\beta ,$ (5)

$y_{ki}-h_k=A\left(x_{ki}-g_k\right).$ (6)

Let Z the intraclass inertia $\left(m\mathrm{,}m\right)$ matrix of $\left\{y_{ki},i=1,\cdots ,n_k;k=1,\cdots ,q\right\}$ :

$Z=\sum _{k=1}^q\sum _{i\in I_k}{p}_{ki}\left(y_{ki}-h_k\right){\left(y_{ki}-h_k\right)}^{\prime }=AW{A}^{\prime }.$ (7)

The rank of Z is equal to the rank of A, $p\le m$ . For $m>p$ , the $\left(m\mathrm{,}m\right)$ matrix Z is not invertible. Then use in this case the pseudoinverse (or Moore-Penrose inverse) of Z, denoted Z⁺, which is equal to the inverse of Z when $m=p$ , to define the pseudodistance denoted $d_{Z^{+}}$ in ${ℝ}^m$ . The denomination pseudodistance is used because Z⁺ is not positive definite. Remind the definition of a pseudoinverse and two theorems [12] .

Definition Let A a $\left(k\mathrm{,}l\right)$ matrix of rank r. The pseudo-inverse of A is the unique $\left(l\mathrm{,}k\right)$ matrix A⁺ such that:

1) $A{A}^{+}A=A$ ,

2) $A^{+}A{A}^{+}=A^{+}$ ,

3) ${\left(A{A}^{+}\right)}^{\prime }=A{A}^{+}$ ,

4) ${\left(A^{+}A\right)}^{\prime }=A^{+}A.$

Theorem 1 Maximal rank decomposition

Let A a $\left(k\mathrm{,}l\right)$ matrix of rank r. Then there exist two full-rank (r) matrices, F of dimension $\left(k\mathrm{,}r\right)$ and G of dimension $\left(r\mathrm{,}l\right)$ ( $rg\left(F\right)=rg\left(G\right)=r$ ) such that $A=FG$ .

Theorem 2 Expression of A⁺

Let $A=FG$ a full-rank decomposition of A. Then $A^{+}=G^{\prime }{\left(F^{\prime }A{G}^{\prime }\right)}^{-1}{F}^{\prime }$ .

Prove now:

Proposition 1 $d_{Z+}^2\left(Aa+\beta ,Ab+\beta \right)=d_{W^{-1}}^2\left(a,b\right)$ .

Proof. $Z=\left(AW\right)A^{\prime }$ . AW and A' are of full-rank p. Applying theorem 2 yields:

$Z^{+}=A{\left({\left(AW\right)}^{\prime }AW{A}^{\prime }A\right)}^{-1}{\left(AW\right)}^{\prime }$ (8)

$=A{\left(A^{\prime }A\right)}^{-1}{\left(W{A}^{\prime }AW\right)}^{-1}{\left(AW\right)}^{\prime }$ (9)

$=A{\left(A^{\prime }A\right)}^{-1}{W}^{-1}{\left(A^{\prime }A\right)}^{-1}{A}^{\prime }.$ (10)

$A^{\prime }{Z}^{+}A=W^{-1}.$ (11)

Note that, when $m=p$ , A is invertible and $Z^{+}={\left(AW{A}^{\prime }\right)}^{-1}=Z^{-1}$ .

$d_{Z^{+}}^2\left(Aa+\beta ,Ab+\beta \right)={\left(A\left(a-b\right)\right)}^{\prime }{Z}^{+}\left(A\left(a-b\right)\right)={\left(a-b\right)}^{\prime }{W}^{-1}\left(a-b\right).$ □

Thus:

Proposition 2 Let A a $\left(m\mathrm{,}p\right)$ matrix, $m>p$ , of rank p and for $k=1,\cdots ,q$ , $i=1,\cdots ,n_k$ , $y_{ki}=A{x}_{ki}+\beta$ . The results of LDA of the dataset $\left\{x_{ki},k=1,\cdots ,q,i=1,\cdots ,n_k\right\}$ with the metrics W⁻¹ on ${ℝ}^p$ are the same as those of LDA of the dataset $\left\{y_{ki},k=1,\cdots ,q,i=1,\cdots ,n_k\right\}$ with the pseudometrics $Z^{+}={\left(AW{A}^{\prime }\right)}^{+}$ .

Applications

Denote $x_i^j$ the value of the variable $x^j$ for individual i belonging to I, $i=1,\cdots ,n$ , $j=1,\cdots ,p$ and $x_i={\left(x_i^1\cdots x_i^p\right)}^{\prime }$ the vector of values of $\left(x^1\mathrm{,}\cdots \mathrm{,}{x}^p\right)$ for individual i. Denote $p_i$ the weight of individual i, such that $\stackrel{n}{\sum _{i=1}}{p}_i=1$ . To perform a factorial analysis of the dataset $\left\{x_i,i=1,\cdots ,n\right\}$ , the difference between two individuals i and i' is measured by a distance $d\left(i\mathrm{,}{i}^{\prime }\right)$ defined on ${ℝ}^p$ associated to a metrics M, such that

$d^2\left(i,i^{\prime }\right)={\left(x_i-x_{i^{\prime }}\right)}^{\prime }M\left(x_i-x_{i^{\prime }}\right).$ (12)

Denote X the $\left(n\mathrm{,}p\right)$ matrix whose element $\left(i\mathrm{,}j\right)$ is $x_i^j$ . Denote D the diagonal $\left(n\mathrm{,}n\right)$ matrix whose element $\left(i\mathrm{,}i\right)$ is $p_i$ .

Perform a factorial analysis of $\left(X\mathrm{,}M\mathrm{,}D\right)$ , for instance principal component analysis (PCA) for continuous variables or MCFA for categorical variables or MDFA for mixed data. Suppose X of rank p. Denote $u_j={\left(u_j^1\cdots u_j^p\right)}^{\prime }$ a unit vector of the j^th principal axis. Denote $c^j=XM{u}_j={\left(c_1^j\cdots c_n^j\right)}^{\prime }$ the j^th principal component. Denote U the $\left(p\mathrm{,}p\right)$ matrix $\left(u_1\cdots u_p\right)$ and C the $\left(n\mathrm{,}p\right)$ matrix $\left(c^1\cdots c^p\right)=XMU$ ; as $u_1\mathrm{,}\cdots \mathrm{,}{u}_p$ are M-orthonormal, $U^{\prime }MU=I$ and

$C=XMU\iff X=C{U}^{\prime }\iff \text{for}i=1,\cdots ,n,x_i=U{c}_i$ (13)

$\iff \text{for}i=1,\cdots ,n,c_i=U^{\prime }M{x}_i.$ (14)

Using the metrics of intraclass inertia matrix inverse, LDA from C is equivalent to LDA from X.

Suppose now that the variable $x^{p+1}=1-x^p$ is introduced; when $x^p$ is the indicator of a modality of a binary variable, $x^{p+1}$ is the indicator of the other modality. Then:

$\left(\begin{array}{c}{x}_i^1\\ ⋮\\ x_i^p\\ x_i^{p+1}\end{array}\right)=\left(\begin{array}{ccc}{u}_1^1& \cdots & u_p^1\\ ⋮& \ddots & ⋮\\ u_1^p& \cdots & u_p^p\\ -u_1^p& \cdots & -u_p^p\end{array}\right)\left(\begin{array}{c}{c}_i^1\\ ⋮\\ c_i^p\end{array}\right)+\left(\begin{array}{c}0\\ ⋮\\ 0\\ 1\end{array}\right)$ (15)

Denote X₁ the $\left(n\mathrm{,}p+1\right)$ matrix whose element $\left(i\mathrm{,}j\right)$ is $x_i^j$ . LDA from C with the metrics of intraclass inertia matrix inverse is equivalent to LDA from X₁ with the metrics of intraclass inertia matrix pseudoinverse.

For instance:

1) If $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ are continuous variables, LDA from X is equivalent to LDA from C obtained by PCA, such as normed PCA, or generalized canonical correlation analysis (gCCA) [13] and MFA which can be interpreted as PCA with specific metrics.

2) If $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ are indicators of modalities of categorical variables, and if MCFA is performed to obtain C, LDA from C with the metrics of intraclass inertia matrix inverse is equivalent to LDA from X with the metrics of intraclass inertia matrix pseudoinverse.

3) Likewise, if $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ are continuous variables or indicators of modalities of categorical variables, and if MDFA [10] is performed to obtain C, LDA from C with the metrics of intraclass inertia matrix inverse is equivalent to LDA from X with the metrics of intraclass inertia matrix pseudoinverse. In this case, other metrics can also be used, such as that of Friedman [14] or that of Gower [15] .

4. Methodology for Constructing a Score

4.1. Ensemble Methods

Consider the problem of predicting an outcome variable y, continuous (in the case of regression) or categorical (in the case of classification) from observable explanatory variables $x^1\mathrm{,}\cdots \mathrm{,}{x}^p$ , continuous or categorical.

The principle of an ensemble method [16] [17] is to build a collection of N predictors and then aggregate the N predictions obtained using:

• in regression: the average of predictions $\widehat{y_i}$ ;

• in classification: the rule of the majority vote or the average of the estimations of a posteriori class probabilities.

The ensemble predictor is expected to be better than each of the individual predictors. For this purpose [16] :

• each single predictor must be relatively good,

• single predictors must be sufficiently different from each other.

To build a set of predictors, we can:

• use different classifiers,

• and/or use different samples (e.g. by bootstrapping, boosting, randomizing outputs) [17] [18] [19] ,

• and/or use different methods of variables selection (e.g. ascending, stepwise, shrinkage, random) [20] [21] [22] [23] ,

• and/or in general, introduce randomness into the construction of predictors (e.g. in random forests [24] , randomly select a fixed number of variables at each node of a classification or regression tree).

In Random Generalized Linear Model (RGLM) [25] , at each iteration,

• a bootstrap sample is drawn,

• a fixed number of variables are randomly selected,

• the selected variables are rank-ordered according to their individual association with the outcome variable y and only the top ranking variables are retained,

• an ascending selection of variables is made using Akaike information criterion (AIC) [26] or Bayesian information criterion (BIC) [27] .

Tufféry [28] wrote that logistic models built from bootstrap samples are too similar for their aggregation to really differ from the base model built on the entire sample. This is in agreement with an assertion by Genuer and Poggi [16] . However, Tufféry suggests the use of a method called “random forest of logistic models” introducing an additional randomness: at each iteration,

• a bootstrap sample is drawn,

• variables are randomly selected,

• an ascending variables selection is performed using AIC [26] or BIC [27] criteria.

Note that this method is in fact a particular case of RGLM method.

Present now the method used in this study to check the stability of the predictor obtained on the entire learning sample.

4.2. Method of Construction of an Ensemble Predictor

The steps of the method for constructing an ensemble predictor are presented in the form of a tree (Figure 2).

At first step, n₁ classifiers are chosen.

At second step, n₂ bootstrap samples are drawn and are the same for each classifier.

At third step, for each classifier and each bootstrap sample, n₃ modalities of random selection of variables are chosen, a modality being defined either by a number of randomly drawn variables or by a number of predefined groups of correlated variables, which are randomly drawn, inside each of which a variable is randomly drawn.

At fourth step, for each classifier, each bootstrap sample and each modality of random selection of variables, one method of selection of variables is chosen, a stepwise or a shrinkage (LASSO, ridge or elastic net) method.

This yields a set of $n_1\times n_2\times n_3$ predictors, which are aggregated to obtain an ensemble predictor.

4.3. Choices Made

To assess accuracy of the ensemble predictor, the percentage of well-classified is currently used. But this criteria is not always convenient, especially in the

present case of unbalanced classes. We decided to use AUC. AUC in resubstitution being usually too optimistic, we used AUC OOB [29] : for each patient, consider the set of predictors built on the bootstrap samples that do not contain this patient, i.e. for which this patient is “out-of bag”, then aggregate the corresponding predictions to obtain an OOB prediction.

Two classifiers were used: logistic regression and LDA with metrics W⁻¹. Other classifiers were tested but not retained because of their less good results, such as random forest-random input (RF-RI) [24] or QDA. The k-nearest neighbors method (k-NN) was not tested, because it was not adapted to this study due to the presence of very unbalanced classes with a too small class size.

1000 bootstrap samples were randomly drawn.

Three modalities of random selection were retained, firstly a random draw of a fixed number of variables, secondly and thirdly a random draw of a fixed number of predefined groups of correlated variables followed by a random draw of one variable inside each drawn group. The number of variables or of groups drawn was determined by optimization of AUC OOB.

Fourth step did not improve prediction accuracy and was not retained.

4.4. Construction of an Ensemble Score

Denote n the total number of patient-months and p the number of variables. Denote $x_i^j$ the value of variable $x^j$ for patient-month i, $i=1,\cdots ,n$ , $j=1,\cdots ,p$ . Each patient-month i is represented by a vector $x_i={\left(x_i^1\mathrm{<mi>\; </mi>}\cdots \mathrm{<mi>\; </mi>}{x}_i^p\right)}^{\prime }$ in ${ℝ}^p$ .

4.4.1. Aggregation of Predictors

In the case of two classes ${\Omega }_1$ and ${\Omega }_0$ , whose barycenters are respectively denoted $g_1$ and $g_0$ , Fisher linear discriminant function

$S_1\left(x\right)={\left(x-\frac{g_1+g_0}{2}\right)}^{\prime }{W}^{-1}\left(g_1-g_0\right)={{\alpha }^{\prime }}_{1}x+{\beta }_1$ (16)

can be used as score function. For logistic regression, the following score function can be used:

$S_2\left(x\right)=\ln \frac{P\left({\Omega }_1|X=x\right)}{P\left({\Omega }_0|X=x\right)}={{\alpha }^{\prime }}_{2}x+{\beta }_2.$ (17)

Remind that, in the case of a multinormal model with homoscedasticity (covariance matrices within classes are equal), when $P\left({\Omega }_1\right)=P\left({\Omega }_0\right)$ , logistic model is equivalent to LDA [17] ; indeed:

$S_2\left(x\right)=\ln \frac{P\left({\Omega }_1|X=x\right)}{P\left({\Omega }_0|X=x\right)}=\ln \frac{P\left({\Omega }_1\right)}{P\left({\Omega }_0\right)}+S_1\left(x\right)=S_1\left(x\right).$ (18)

So we used the following method to aggregate the obtained predictors:

1) the score functions obtained by LDA are aggregated by averaging; denote now S₁ the averaged score;

2) likewise the score functions obtained by logistic regression are aggregated by averaging; denote S₂ the averaged score;

3) a combination of the two scores, $\lambda S_1+\left(1-\lambda \right)S_2$ is defined, $0\le \lambda \le 1$ ; a value of $\lambda$ that maximizes AUC OOB is retained; denote S₀ the optimal score obtained by this method.

If s is an optimal cut-off, the ensemble classifier is defined by:

If $S_0\left(x\right)>s$ , x is classified in ${\Omega }_1$ ; (19)

if not, x is classified in ${\Omega }_0$ . (20)

4.4.2. Definition of a Score from 0 to 100

The variation scale of the score function $S_0\left(x\right)$ was reduced from 0 to 100 using the following method. Denote:

$S_0\left(x\right)={{\alpha }^{\prime }}_{0}x+{\beta }_0=\sum _{j=1}^p{\alpha }_0^j{x}^j+{\beta }_0.$ (21)

Denote for $j=1,\cdots ,p$ :

$P_j=\left|{\alpha }_0^j\right|\left(\underset{1\le i\le n}{\max }{x}_i^j-\underset{1\le i\le n}{\min }{x}_i^j\right)$ (22)

and

$P=\sum _{j=1}^p{P}_j=\sum _{j=1}^p\left|{\alpha }_0^j\right|\left(\underset{1\le i\le n}{\max }{x}_i^j-\underset{1\le i\le n}{\min }{x}_i^j\right).$ (23)

Let $m^j$ the minimal value of the variable $x^j$ if ${\alpha }_0^j>0$ , or its maximal value if ${\alpha }_0^j<0$ .

Denote $S\left(x\right)$ the “normalized” score function, with values from 0 to 100, defined by:

$S\left(x\right)=\frac{100}{P}\sum _{j=1}^p{\alpha }_0^j\left(x^j-m^j\right)$ (24)

$=100\sum _{j=1}^p\frac{{\alpha }_0^j\left(x^j-m^j\right)}{{\sum }_{k=1}^p\left|{\alpha }_0^k\right|\left({\max }_{1\le i\le n}{x}_i^k-{\min }_{1\le i\le n}{x}_i^k\right)}$ (25)

$={\alpha }^{\prime }x+\beta$ , with $\left(\begin{array}{c}\beta \\ {\alpha }_1\\ ⋮\\ {\alpha }_p\end{array}\right)=\left(\begin{array}{c}-\frac{100}{P}{\sum }_{j=1}^p{\alpha }_0^j{m}^j\\ 100\frac{{\alpha }_0^1}{P}\\ ⋮\\ 100\frac{{\alpha }_0^p}{P}\end{array}\right).$ (26)

4.4.3. Measure of Variables Importance

Explanatory variables are not expressed in the same unit. To assess their importance in the score, we used “standardized” coefficients, multiplying the coefficient of each variable in the score by its standard deviation. These coefficients are those associated with standardized variables and are directly comparable. For all variables, the absolute values of their standardized coefficient, from the greatest to the lowest, were plotted on a graph. The same type of plot was used for groups of correlated variables, whose importance is assessed by the sum of absolute values of their standardized coefficients.

4.4.4. Risk Measure by an Odds-Ratio

Define a risk measure associated to a score s by an odds-ratio $O{R}_1\left(s\right)$ :

$O{R}_1\left(s\right)=\frac{P\left(Y=1|S>s\right)}{P\left(Y=0|S>s\right)}\frac{P\left(Y=0\right)}{P\left(Y=1\right)}=\frac{P\left(S>s|Y=1\right)}{P\left(S>s|Y=0\right)}=\frac{Se\left(s\right)}{1-Sp\left(s\right)}.$ (27)

An estimation of $O{R}_1\left(s\right)$ , also denoted $O{R}_1\left(s\right)$ , is $\frac{n_1}{n_0}\times \frac{N_0}{N_1}$ with $n_k=\#\left\{S>s\right\}\cap \left\{Y=k\right\}$ and $N_k=\#\left\{Y=k\right\}$ , $k=0,1$ .

Note that:

• $O{R}_1\left(s\right)$ decreases when $Se\left(s\right)$ decreases and $Sp\left(s\right)$ is constant. In practice, the decrease will be much smaller when there are many observations;

• $O{R}_1\left(s\right)$ is not defined when $Sp\left(s\right)$ is equal to 1.

For these reasons, the following definition can also be used:

$O{R}_2\left(s\right)=\underset{t\le s:\mathrm{<mi>\; </mi>}O{R}_1\left(t\right)<\infty }{\max }O{R}_1\left(t\right).$ (28)

Note that $O{R}_1$ is the slope y/x of the line joining the origin to the point $\left(x\mathrm{,}y\right)$ of the ROC curve. In the case of an “ideal” ROC curve, supposed continuous above the diagonal line, assuming that there is no vertical segment in the curve, this slope increases from point $\left(\mathrm{1,1}\right)$ , corresponding to the minimal value of score, to point $\left(\mathrm{0,0}\right)$ , corresponding to its maximal value; the case of a vertical segment (Se decreases, Sp is constant), occurring when the score of a patient with event is between those of two patients without event, is particularly visible in the case of a small number of patients and also justifies the definition of $O{R}_2$ , whose curve fits that of $O{R}_1$ .

For very high score values, when n₀ or n₁ are too small, the estimation of $O{R}_1$ is no longer reliable. A reliability interval of the score could be defined, depending on the values of n₀ and n₁.

5. Results

5.1. Pre-Processing of Variables

5.1.1. Winsorization

To avoid problems related to the presence of outliers or extreme data, all continuous variables were winsorized using the 1^st percentile and the 99^th percentile of each variable as limit values [30] . We chose this solution because of the large imbalance of the classes (317 patients with event against 21,065 with no event, so there is a ratio of about 1 to 66). The elimination of extreme data would have led to decrease the number of patients with event.

5.1.2. Transformation of Variables

Among qualitative variables, two are ordinal: the NYHA class with 4 modalities and the number of myocardial infarction (no. MI) with 5 modalities. In order to preserve the ordinal nature of these variables, we chose to use an ordinal encoding. For NYHA, we therefore associated 3 binary variables: NYHA ≥ 2, NYHA ≥ 3 and NYHA ≥ 4. In the same way, for the no. MI, we considered 4 binary variables: no. MI ≥ 2, no. MI ≥ 3, no. MI ≥ 4 and no. MI ≥ 5.

On the other hand, continuous variables were transformed in the context of logistic regression. For each continuous variable, a linearity test was performed using the method of restricted cubic splines with 3 knots [31] . A cubic spline restricted with 3 knots is composed of a linear component and a cubic component. Linearity testing is to test, under the univariable logistic model, the nullity of the coefficient associated with the cubic component. To do this, we used the likelihood ratio test. The results of linearity tests are given in Table 1 (p-value 1).

At 5% level, linearity was rejected for 9 of 16 continuous variables. For each of these 9 variables, we represented graphically the relationship between the logit (natural logarithm of the ratio probability of event/probability of non-event) and the variable. An example of graphical representation is given for potassium: we observe a quadratic relationship between the logit and the potassium (Figure 3). In agreement with the relationship observed, we applied a simple, monotonous or quadratic transformation function to each of the 9 variables. The transformation function applied to each variable is given in Table 1.

For hematocrit and the three variables of eGFR, the relationship is clearly monotonous. So we considered some simple monotonic transformation functions as $f\left(x\right)=x^a$ with $a\in \left\{-\mathrm{2,}-\mathrm{1,}-\mathrm{0.5,0.5,1,2}\right\}$ or $f\left(x\right)=\ln \left(x\right)$ , then we retained for each variable the transformation for which the likelihood under univariable logistic model was maximal (minimal p-value).

For other variables not checking linearity, namely potassium, the three blood pressure measures (systolic, diastolic and mean), and heart rate, the relationship between the logit and the variable was rather quadratic. We therefore applied a quadratic transformation function ${\left(X-k^{\mathrm{<mo>\; *\; </mo>}}\right)}^2$ with $k^{\ast }$ an optimal value determined by maximizing likelihood under univariable logistic model. To compare, we also used the criterion of maximal AUC to determine an optimal value. These results are presented in Table 2. Notice that the optimal values determined by the two methods are the same for systolic BP, diastolic BP and heart rate and are very close for potassium and mean BP.

Also note that the transformation applied to potassium allows to take into account both hypokalemia and hyperkalemia, two different clinical situations pooled here that may increase the risk of death and/or hospitalization measured by the score.

To verify that the transformation of the variables was good, a linearity test for each transformed variable was performed according to the previously detailed principle. All tests are not significant at the 5% level (see Table 1, p-value 2).

5.2. Ensemble Score

5.2.1. Ensemble Score by Logistic Regression

As a first step, we applied our methodology with the following parameters:

• use of a single classification rule, logistic regression ( $n_1=1$ ),

• draw of 1000 bootstrap samples ( $n_2=1000$ ),

• random selection of variables according to a single modality ( $n_3=1$ ).

Three modalities for the random selection of variables were defined:

• 1^st modality: random draw of m variables among 32,

• 2^nd modality: random draw of m groups among 18, then one variable from each drawn group,

• 3^rd modality: random draw of m groups among 24, then one variable from each drawn group.

The groups of variables considered for each modality are presented in Table 3. For modalities 2 and 3, we formed groups of variables based on correlations between variables. For the second modality, we gathered for example in the same group hemoglobin, hematocrit and ePVS because of their high correlations. For the third modality, the same groups were used, except for the two variables linked to hospitalization for HF, the four variables linked to the no. MI and the three variables related to the NYHA class, for which each binary variable was considered as a single group.

For each modality, an ensemble score was built for all possible values of m and the one that gave maximal AUC OOB was selected. In Table 4 are reported the results obtained for each modality with the optimal m. The best result was obtained for the third modality, with AUC OOB equal to 0.8634.

The ensemble score by logistic regression, denoted $S_2\left(x\right)$ , obtained by averaging the three ensemble scores that we constructed, gave slightly better results, with AUC OOB of 0.8649.

5.2.2. Ensemble Score by LDA for Mixed Data

The same methodology was used by simply replacing the classification rule (logistic regression) by LDA for mixed data and keeping the same other settings.

Again, for each modality, we searched the optimal m parameter. The obtained results are presented in Table 5.

As for logistic regression, the best results were obtained for the third modality, with AUC OOB equal to 0.8638.

The ensemble score by LDA, denoted $S_1\left(x\right)$ , yielded better results with AUC OOB equal to 0.8654.

5.2.3. Ensemble Score Obtained by Synthesis of Logistic Regression and LDA

The final ensemble score denoted $S_0\left(x\right)$ , obtained by synthesis of the two ensemble scores $S_1\left(x\right)$ and $S_2\left(x\right)$ presented previously, provided the best results with AUC equal to 0.8733 in resubstitution and 0.8667 in OOB.

This ensemble score corresponds to the one obtained by applying our methodology with the following parameters:

• two classification rules are used, logistic regression and LDA for mixed data ( $n_1=2$ ),

• 1000 bootstrap samples are drawn ( $n_2=1000$ ),

• m variables are randomly selected according to three modalities ( $n_3=3$ ).

The scale of variation of the score function $S_0\left(x\right)$ was reduced from 0 to 100 according to the procedure described previously. We denote this “normalized” score $S\left(x\right)$ .

In Table 6, we present the “raw” and “standardized” coefficients associated with each of the variables in the score function $S_0\left(x\right)$ and the “normalized” score function $S\left(x\right)$ .

5.2.4. Importance of Variables in the Score

To have a global view of the importance of the variables in the “normalized” score, we represented on a graph the absolute value of standardized coefficient associated with each variable, from the largest value to the smallest (see Figure 4). Note that the most important variables are heart rate, NYHA class ≥ 3 and history of hospitalization for HF in the previous month. On the other hand, variables such as weight, no. MI ≥ 5 or BMI do not play a large part in the presence of others.

The same type of graph was made to represent the importance of the groups of variables in configuration 2 defined by the sum of the absolute values of the “standardized” coefficients associated with the variables of the group, from the largest sum to the smallest (see Figure 4). Note that the two most influential groups are “NYHA” (NYHA ≥ 2, NYHA ≥ 3 and NYHA ≥ 4) and “History of hospitalization for HF” (hospitalization for HF in the previous month and hospitalization for HF during life). Three important groups follow: “Hematology” (ePVS, hemoglobin, hematocrit), “Heart rate” and “Renal function” (creatinine and three formulas of eGFR). The least important groups of variables are “Obesity” (weight, BMI) and “Gender”.

5.2.5. Risk Measure by an Odds-Ratio

We represented the variation of $n_0$ , $n_1$ , $Se\left(s\right)$ , $1-Sp\left(s\right)$ , $O{R}_1\left(s\right)$ and $O{R}_2\left(s\right)$ according to the score s (Table 7). For score values $s>49.1933$ , $n_1$ is less than or equal to 30. Thus, beyond this threshold value 49.1933, $O{R}_1$ is no longer very reliable. We therefore defined as reliability interval of the $O{R}_1$ and $O{R}_2$ functions $\left[\mathrm{0\; <mo>\; ;\; </mo>\; 49.1933}\right]$ .

We represented the variation of odds-ratio $O{R}_1$ and $O{R}_2$ in this reliability interval (Figure 5). By reading the graph, for a patient with a score of 40 for

example, $\frac{P\left(Y=1|S>40\right)}{P\left(Y=0|S>40\right)}$ is about 15 times higher than $\frac{P\left(Y=1\right)}{P\left(Y=0\right)}$ .

6. Conclusions and Perspectives

In this article, we presented a new methodology for constructing a short-term event risk score in heart failure patients, based on an ensemble predictor built using two classification rules (logistic regression and LDA for mixed data), 1000 bootstrap samples and three modalities of random selection of variables. This score was normalized on a scale from 0 to 100. AUC OOB is equal to 0.8667. Note

that an important variable such as potassium that does not appear in other scores (as SPIM risk score) is taken into account in this score.

Moreover, we defined a measure of the importance of each variable and each group of variables in the score and defined an event risk measure by an odds-ratio.

Due to the nature of the data available (data obtained from the EPHESUS study), we had to define the short term to 30 days in order to have enough patients with HF event. It would be better to have data of patients with shorter intervals, in order to have data the closest possible of an event and eventually improve the quality of the score. When such data will be available, it will be interesting to apply the same methodology to construct a new score.

Furthermore, we proved a property of linear discriminant analysis for mixed data.

Finally, this methodology can be adapted to the case of a data stream. Suppose that new data for heart failure patients arrives continuously. Data can be allocated to bootstrap samples using Poisson bootstrap [32] . The coefficients of each variable in each predictor based on logistic regression or binary linear discriminant analysis can be updated online using a stochastic gradient algorithm. Such algorithms are presented in [33] for binary LDA and [34] for logistic regression; they use online standardized data in order to avoid a numerical explosion in the presence of extreme values. Thus the ensemble score obtained by averaging can be updated online. To the best of our knowledge, it is the first time that this problematics is studied in this context.

Acknowledgements

Results incorporated in this article received funding from the Investments for the Future program under grant agreement No ANR-15-RHU-0004.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Duarte, K., Monnez, J.-M. and Albuisson, E. (2018) Methodology for Constructing a Short-Term Event Risk Score in Heart Failure Patients. Applied Mathematics, 9, 954-974. https://doi.org/10.4236/am.2018.98065

References