Abstract

I. Introduction

Cardiac resynchronization therapy (CRT) is indicated for patients suffering from drug-refractory congestive heart failure (CHF) associated with intraventricular dyssynchrony [1]. CRT improves hemodynamic parameters, ejection fraction or distance covered in the 6 minutes walking test [2]. Furthermore, CRT has shown to decrease hospitalizations for patients treated with the implantable devices (ID). Although the efficiency of this treatment has been proven, 20 to 30% of patients show either no improvement or worsening of their symptoms [3].

Individual follow-up of implanted patients is a key to understand the difference between responders and non-responders, and to prevent severe health degradation. Besides regular follow-up visits, during the post-operative period, an everyday follow-up is possible with the new IDs recently developed for CRT. They offer an increased storage capability of data acquired by the ID, providing information on the ID itself (e.g. event counters of pacing and sensing activities) or on the state of the patient (e.g. arrhythmias, electrograms) and on its activity [4]. Recorded data are very promising towards the home monitoring of patients, the prediction of adverse events or the reduction of hospitalizations. However, this source of information is under-exploited because data are large, multivariate, time-dependent and heterogeneous, and consequently difficult to interpret for caregivers.

The objective of the present study is to propose a methodology to process this amount of multivariate data, in order to i) evaluate and extract the information content of the time-dependent data downloaded from the pacemaker memory, ii) define synthetic indices which are easy to interpret and iii) characterize and compare different populations of patients. Given the dimensionality of the recorded data, methods of data reduction are investigated. The interest of the multidimensional analysis of the data recorded in the ID memory to objectively assess the patients’ response to the therapy and the validity of the exploratory techniques to process these data have been shown in two previous studies, using principal component analysis (PCA) [5] and multiple correspondence analysis (MCA) associated with a spatio-temporal fuzzy coding of the time-series [6]. The former method has been successfully used to differentiate a test population (patients with rate-responsive pacemakers) from a population of patients suffering from CHF by jointly exploiting a number of physiological variables and using a simple representation for the temporal dimension. Providing an appropriate adaptation of its table of analysis, MCA has been successfully applied to the analysis of the evolution of time-series across time, and is then used here as well. MCA performs a reduction of the dimensionality of the data and provides synthetic indices called “factorial axes”. A plane formed by two factorial axes is called “factorial plane”. Each patient is finally represented by trajectories on the factorial plane.

From a methodological point of view, several questions are raised: i) how to link the factorial axes with the variables acquired from the ID? ii) do patients with similar trajectories on the factorial plane have a similar clinical state, and if yes, how to cluster patients according to their evolution in the factorial plane? and iii) are the obtained clusters consistent with the clinical data available from the patients? The study addresses the clustering of trajectories with different numbers of points in the factorial plane, which implies the choice of appropriate distance measure and clustering method. This problem is related to temporal clustering (i.e. the clustering of time-series) [7], [8].

The paper is organized as follows. In section II, the clinical protocol and the ID data are presented. Then, in section III, the proposed methodology is described and the issues arising from the clustering of trajectories are presented and solutions are proposed. The proposed methodology is tested on simulated datasets and applied on the real recorded data in Section VI. The validity of the obtained clusters is evaluated in comparison with the medical records of patients participating in the protocol.

II. Data recorded by cardiac implantable devices

III. A methodology for clustering multivariate time-series

An overview of the proposed methodology is provided in Figure 1. The analysis is first performed on a reference population to determine significant factorial axes and clusters of similar evolutions in the factorial planes. In a follow-up situation, new data would regularly be retrieved from patients’ implantable cardiac devices. The results of the analysis (i.e. factorial axes and clusters) would then be used by projecting the new data on the factorial planes and assigning the evolution of the patient to an existing cluster towards the diagnosis (i.e. the patient is improving or degrading). In this paper, only the analysis is presented.

Open in a separate window

Fig. 1

Overview of the proposed methodology. The analysis is performed on the reference population and leads to the determination of clusters of similar evolutions – according to an appropriately chosen dissimilarity measure – on the factorial plane defined by the smoothed multiple correspondence analysis (SMCA). A practical application of this methodology would be the projection of subsequent follow-up data as supplementary individuals on the factorial axes defined during the analysis and the assignation of the obtained trajectories to the “closest” cluster, in the sense of the chosen distance.

The analysis consists of 2 successive steps, namely fuzzy coding of the data and multidimensional analysis with smoothed multiple correspondence analysis (SMCA), described in the following subsections.

A. Space-time fuzzy coding

A coding of the recorded data is required, as MCA has been at first conceived for categorical variables. MCA exploits disjunctive tables Z = (z_ij)_{(i,j)[1,R]×[1,C]} where z_ij is the membership value of the i^th object to the j^th modality. As a coding of the multivariate time-series, a fuzzy space-time windowing, defined by Loslever and Bouilland for characterizing and coding biomechanical temporal data [12], is proposed. Instead of an indicator matrix Z (i.e. z_ij {0,1}), the MCA analyses a fuzzy version of Z where z_ij [0,1] with the condition Σ_jJv, z_ij = 1, J_v being the set of modalities of the v^th attribute (variable). In statistics, one modality of a variable is one possible level of this variable. In classical crisp coding, continuous variables are divided into several modalities (whose number depends on the distribution of the variable). For example, the variable “age” of a population can be split into 3 modalities “age ≤ 30”, “30 < age < 60” and “age ≥ 60”.

As depicted in Figure 2, the fuzzy space-time windowing divides the time domain of each variable into N_T overlapping windows. The membership value of the q^th time sample t_q to the time window T_j is denoted μ_Tj(t_q) and falls between [0,1]. The membership values of each data point in the N_T time windows meet the condition ${\sum }_{j=1}^{N_T}{\mu }_{T_j}(t_q)=1$. The amplitudes of each variable (referred to as spatial domain) are coded in a similar manner with N_A spatial fuzzy windows verifying the same properties.

Open in a separate window

Fig. 2

Temporal and spatial (amplitude) fuzzy coding of a continuous signal.

From the membership values in the time and the spatial (amplitude) domains, the membership value of the space-time window $W_{i,j}^n$, for a given time-series (signal) TS_k and the variable V_n, is denned as:

$${\mu }_{W_{i,j}^n}(T{S}_k)=\frac{{\sum }_{q=1}^Q{\mu }_{T_j}(t_q)·{\mu }_{A_{i,n}}(V_n(t_q))}{{\sum }_{q=1}^Q{\mu }_{T_j}(t_q)}$$

(1)

where V_n(t_q) is the value taken by the n^th variable at time unit t_q, μ_Tj(t_q) is the membership value of the time window T_j for the time unit t_q, μ_Ai,n(V_n(t_q)) is the membership value of the i^th space window A_i,n for the V_n(t_q) value, and Q is the number of time units in TS_k. With this definition, ${\mu }_{W_{i,j}^n}$ (TS_k) is the weighted average of the space membership values with the time membership values as weights and verifies ${\sum }_{i=1}^{N_A}{\mu }_{W_{i,j}^n}(T{S}_k)=1$. This property is required to maintain the statistical context and to allow ${\mu }_{W_{i,j}^n}$ (TS_k) to be interpreted as the frequency of appearance of the signal in the space-time window $W_{i,j}^n$.

IV. Application of the fuzzy coding and multidimensional analysis on the recorded database

The first two steps of the proposed methodology, namely the fuzzy spatio-temporal coding and the SMCA, are applied to the time-series available in the recorded database. In this study, a fuzzy window set T = {T₁, · · ·, T_J, · · ·, T_NT}, where each T_j is 7-day long, is considered. The length of the fuzzy time windows has been chosen by considering the patients’ behaviors, being quite similar from one week to another. Each trajectory, representing a three-month length period (i.e. around 13 7-day long), links then around 13 points.

For each variable V_n, N_A = 3 modalities are considered: “Low” corresponding to the spatial fuzzy window A_1,n = [−∞, median(V_n)], “Medium” corresponding to A_2,n = [prctile(V_n, 2.5), prctile(V_n, 97.5)] and “High” corresponding to A_3,n = [median(V_n), +∞], where prctile(V_n,p) is the p^th percentile of the variable V_n. The 3 spatial fuzzy windows are denoted with the suffixes “−H” for the higher level (modality), “−M” for the medium level and “−L” for the lower level.

The elements of the proximity matrix P are defined as:

performing a temporal average of the statistical individuals and respecting the condition ${\sum }_{j=1}^L{p}_{ij}=1$.

The protocol provided 58 records for 41 patients. SMCA is thus applied to 793 statistical individuals related to these 41 patients, i.e. to an array of 793 rows and 39 columns (13 variables with N_A = 3 spatial fuzzy windows).

The 1^st and 2^nd factorial axes represent respectively 68.0% and 15.7% of the total variance, showing that the majority (83.7%) of the information is contained in the first factorial plane. So a great part of the variance of the initial data can be represented by only two factorial axes when 13 variables were initially considered. These first two factors define synthetic indices for patient follow-up but their interpretation, from physiological and functional points of view, has to be performed. Consequently, this study will focus on the clustering of trajectories on this factorial plane.

The projections of the variables and individuals of analysis on the first factorial plane of the SMCA are provided in Figures 3 and ​and4,4, respectively. The first axis is mainly defined by the lower (−L) and higher (−H) modalities of Imp_E, Dur_E, NbCardCyc_E, NbBreaths_E, ACC_E and ImpRate, which reflect the time spent in exercise and the intensity of the efforts made by the patient. Consequently, in Figure 4, the more the individuals are located to the right of the plane, the lower is the time spent in exercise and the less important are the efforts they make. The second axis is mainly defined by the medium (−M) and extreme (−H, −L) modalities of ImpMinVent_E, ImpOverAcc, Imp_R and HeartRate_E, which define the ventilation activity in terms of amplitude, frequency and flow rate, especially in rest. This axis can be interpreted as an “axis of cardiovascular efficiency”. In Figure 4 the more the individuals are located to the lower part of the plane, the less important is their ventilation in rest (i.e. their cardiovascular system is more “efficient”), independently of the daily activity duration and intensity (i.e. of the position along the first axis).

Open in a separate window

Fig. 3

Variables of analysis represented on the first plane of the Smoothed Multiple Correspondence Analysis (SMCA). For each variable, only the first (Low, −L) and the third (High, −H) levels are labelled, unlabelled squares correspond to the second levels (Medium, −M).

Open in a separate window

Fig. 4

Individuals of analysis represented on the first plane of the Smoothed Multiple Correspondence Analysis (SMCA). Each record of each patient (corresponding a three-month period) is represented by one trajectory. The first point in time for each trajectory is marked up by a circle.

As it can be seen on Figure 4, trajectories present different locations and evolutions on the factorial plane, with a high overlapping.

V. Clustering trajectories on the factorial plane

The aim of the present study is to cluster patients according to their clinical state during the follow-up period, which, in terms of methodology, corresponds to the clustering of trajectories on the factorial plane according to their location and evolution. Considering the characteristics of the trajectories, the clustering methodology has to address the following issues:

• Unsupervised: no a priori knowledge on the clusters is required.

• Similarity measure: relevance in comparing trajectories on the factorial plane having different numbers of points and possibly subjects to nonlinear spatio-temporal deformations.

• Location and evolution: both locations and evolutions of the trajectories on the factorial plane have to be taken into account, as they are both informative.

Each of the previous three points is described in the following sections.

VI. Results obtained with the clustering

Before applying the proposed clustering method, namely agglomerative hierarchical clustering with complete linkage and the dissimilarity matrix based on dynamic time warping, it is tested on two simulated datasets, described in the following sections.

A. Tests on simulated datasets

The aim of the first test is to explore the capability of the proposed similarity measure to take into account both locations and evolutions of the trajectories on the factorial plane. The second test interests in the complete methodology described in Figure 1, from the fuzzy space-time coding to the clustering.

1. Test on the similarity measure

From two trajectories selected among the trajectories obtained with the real dataset (cf. Figure 4), three prototypes of trajectories are computed, the third being obtained by inverting the time samples of the first one.

Nine trajectories are simulated from the first prototype at different locations in the factorial plane, 6 from the second prototype and 12 from the third prototype. For each trajectory, a white noise is added to obtain slightly different trajectories for the same prototype. The simulated dataset of 27 artificial trajectories is represented in Figure 5. Figure 5 shows that the simulation reproduces the characteristics of the real dataset, namely the overlapping of several trajectories of different shapes at different locations on the factorial plane.

Open in a separate window

Fig. 5

Simulated dataset of 27 trajectories obtained from three real trajectories and represented on the first factorial plane. The first point in time for each trajectory is marked up by a circle.

The clustering method proposed above is applied to the simulated dataset with three dissimilarity measures based on the DTW which is alternately computed on: i) the coordinates of the trajectories, ii) the derivatives of the trajectories, iii) the derivatives with addition of the Euclidian distance between the means of the coordinates as proposed in equation 5.

Figure 6 illustrates the resulting dendrogram for each dissimilarity measure and provides the mean silhouette value against the number of clusters. A dendrogram is a tree-like plot where each step of hierarchical clustering is represented as a fusion of two branches of the tree into a single one. The branches represent clusters obtained on each step of hierarchical clustering. This representation eases the choice of the number of clusters. The dendrogram is cut at the threshold (horizontal dashed line) producing the number of clusters corresponding to the maximum mean silhouette value.

Open in a separate window

Fig. 6

Dendrograms and mean silhouette value vs. the number of clusters, for three dissimilarity measures based on the DTW. The clustering is applied to the simulated dataset of 27 trajectories. The dendrogram is cut at the threshold whose value is obtained after the maximum of the mean silhouette.

For each of the three measures, the mean silhouette value S presents one global maximum, indicating a value for the number of clusters K to be chosen, and leads to a different clustering. For the DTW on coordinates, max_K S is obtained for K = 3 and it tends to group the trajectories close in the sense of their location on the factorial plane. For the DTW on derivatives, max_K S is obtained for K = 3, but it regroups trajectories only according to their evolution. For the DTW on derivatives with the addition of the Euclidian distance between coordinates, max_K S is obtained for K = 9 and each obtained cluster is composed of trajectories with similar evolutions and close locations. Despite the overlapping of several trajectories with different evolutions at the same location, the proposed dissimilarity measure is able to group similar trajectories. Moreover, it is suitable to distinguish between shapes (independent of any notion of time) and evolutions (with start- and end-points), as trajectories created with the first prototype and with the third prototype are assigned to two different clusters, even if they have similar shapes.

2. Test of the complete method

In this section, the noise robustness of the proposed methodology, described in Figure 1, is tested. Among the trajectories obtained with the real database (cf. Figure 4), seven trajectories with different evolutions and locations are selected. The corresponding time-series (i.e. the 13 variables of analysis) are retrieved. For each selected trajectory, a white noise is added independently to each of its 13 variables with a given signal-to-noise ratio (SNR). The original time-series and 10 realizations of the noisy time-series are coded by fuzzy space-time coding and constitute the individuals of analysis for the SMCA. The 77 resulting trajectories are then clustered by AHC with the DTW on derivatives with the addition of the Euclidian distance between coordinates as a dissimilarity measure (cf. equation 5). Figure 7 presents the mean silhouette value against the number of clusters for SNR = 3 dB, and the clusters obtained for the maximum value of the silhouette value being K = 7. Despite the noise, the correct number of clusters is determined by the internal validation indices, the trajectories being grouped according to their evolution and location.

Open in a separate window

Fig. 7

Test of the complete method on 77 trajectories computed from 7 selected trajectories of the real dataset. For each selected trajectory, a white noise is added independently to each of its 13 variables with a given signal-to-noise ratio (SNR = 3 dB), and the resulting 77 trajectories are clustered. For each cluster, individuals of analysis (in gray) and individuals of the given cluster (in black) are represented in the first plane of the SMCA. The first point in time for each trajectory is marked up by a circle.

It appears that above SNR = 0 dB the clustering is correct and is not disturbed by the noise added to the time-series and that under SNR = 0 dB, S is unable to provide the correct value of K. The method seems then robust to noise on the variables of analysis, which can be explained by two of its steps: the time fuzzy coding and the smoothing performed during the SMCA. The mean silhouette value increases when the length of the time fuzzy windows increases and is higher with the SMCA than with the MCA. The averaging performed by both the fuzzy coding and the SMCA smoothes the time-series and the trajectories, respectively, improving the performance of the clustering.

The previous two tests, performed on datasets close to the real database, have proven the validity of the proposed approach in terms of clustering methodology and dissimilarity measure relevant for the processed trajectories, and of noise robustness.

B. Performance of the clustering on the recorded database

In this section, the recorded database, constituted of 58 trajectories on the first factorial plane, is clustered with the proposed approach. Figure 8 provides the mean silhouette value S against the number of clusters K. The number of clusters is chosen at K = 10, after the maximum value of S. The resulting clusters are provided in Figure 8, where for each cluster individuals of analysis are in gray and individuals of the given cluster are in black. One can notice that trajectories within a given cluster have visually similar evolutions and close locations.

Open in a separate window

Fig. 8

Mean silhouette value S against the number of clusters K. Ten clusters are determined with the proposed methodology for the recorded database (max_K S for K = 10). For each cluster, individuals of analysis (in gray) and individuals of the given cluster (in black) are represented in the first plane of the SMCA. Each record of each patient (corresponding to a three-month period) is represented by one trajectory. The first point in time for each trajectory is marked up by a circle.

To evaluate the methodology, the resulting clusters have to be compared with the appreciation of the physicians on the evolution of the patients during each of the three-month length periods. For the present protocol, information is available on the global state of each patient, updated at the end of each 3-month period of the follow-up, and each “adverse event” is reported (date and type). In the present study, the records with adverse events other than cardiac events are discarded, as a non cardiac event like a fall or a bronchitis can have very different effects on patients with RHF. As no everyday report of patients’ clinical state is available, the difficulty resides in defining the actual evolution of a given patient: has a patient undergoing a cardiac adverse event at the beginning of his/her three-month length period and recovering rapidly after hospitalization a favorable or an unfavorable evolution? In this context, indices like the specificity or the sensitivity are difficult to compute. Consequently, the present study can only focus on the relation between the different areas of the factorial plane defined by SMCA and the health of the patients as reported at the end of each follow-up period. More data would be necessary to study the evolutions of patients on this factorial plane.

According to the analysis performed on the position of the individuals of analysis relatively to the position of the modalities of analysis (cf. Figures 3 and ​and4),4), the evolution of a patient (during one of his/her three-month length period) is a priori i) favorable if the corresponding trajectory is located on the left of the plane or evolves from the right to the left of the plane, and ii) unfavorable if the corresponding trajectory is located on the right part of the plane or evolves from the left to the right of the plane. The study of each cluster is necessary to confirm this division of the factorial plane into several areas with different meanings in terms of patient’s health.

On the left of the plane, clusters 2, 4 and 5 contain 13 trajectories in total. The medical reports indicate that all the corresponding patients were, during the given period, in a favorable state. On the right of the plane, clusters 8 and 9 contain 5 trajectories in total, associated with patients all undergoing a health deterioration during the given period. These five clusters enable the definition of two distinct areas in the factorial plane, the left one associated with a favorable state and the right one with a health deterioration.

The other five clusters are interesting as their trajectories show transitions from one area to the other. Cluster 1 comprises 27 trajectories, with small evolutions on the upper half-plane, overlapping both right and left half-planes. During the given period, patients whose trajectories belong to cluster 1 were all in good health, except one patient. According to the medical records on this patient, he was tired during the period of interest (3 months post-op) and died 1 month after the end of this period, the date of his/her adverse event is not reported. Figure 9 shows that the corresponding trajectory (solid line) is evolving from the upper-plane to the right of the plane, with a loop around the 7^th post-op week. The trajectory in cluster 3 is very specific, evolving bottom-to-top firstly and right-to-left secondly. This patient underwent an adverse event (heart failure) 28 days after the beginning of the period of interest, corresponding to the changing of direction in his/her trajectory, and recovered rapidly which explains why the trajectory evolves from right to left. Cluster 6 groups 6 trajectories evolving in the same direction (bottom-right to top), associated with patients all, except one, presenting a favorable state during the given period. This patient (bold line in Figure 9) underwent a cardiac adverse event around 2 months after the beginning of the period of interest. He was hospitalized for 1 week and recovered rapidly, as can be seen from his/her trajectory finally evolving right to left. In Figure 9, the arrow corresponds to the reported date of the adverse event, and one can notice that a changing in the direction of the trajectory has occurred around 5 weeks before the adverse event. Cluster 7 comprises 2 trajectories evolving from top-right to top-left of the plane and corresponding to one patient with a favorable evolution and one patient with an unfavorable evolution. This patient underwent a cardiac adverse event around 2 weeks before the end of the period of interest, as can be seen with the abrupt U-turn in his/her trajectory (dashed line in Figure 9). Cluster 10 groups 4 trajectories with two different evolutions, two trajectories have very small dynamics, associated with patients in good health, and the two others have large evolutions from the top to the right of the plane, associated with patients undergoing adverse events during the given period.

Open in a separate window

Fig. 9

Three specific trajectories of patients undergoing an adverse event during the period of interest. Solid line: trajectory in cluster 1, bold line: trajectory in cluster 6, dashed line: trajectory in cluster 7. Reported dates of adverse events are indicated by an arrow. The first point in time for each trajectory is marked up by a circle.

The clustering performed on the real dataset provides 10 clusters, grouping trajectories with similar evolutions and close locations, as expected. The obtained clusters seem consistent with the medical records of the patients, but no quantitative evaluation is available. Two areas in the factorial plane have been identified, the bottom-right quarter-plane is related with a health degradation and the bottom-left one with a stable clinical state. The other two areas, although being distinguished by the clustering algorithm, are not so easily identifiable.

VII. Discussion and conclusion

This study was designed to show the informative potential of acceleration and impedance data recorded in implantable devices and to evaluate the appropriateness of a multiple correspondence analysis in this particular context. A clinical protocol has been designed to provide data on patients that suffer from heart failure.

MCA has been chosen as it is a multivariate method that exhibits linear and non-linear relationships between variables. However, MCA demands to transform continuous variables into nominal ones, i.e. to code amplitude and time domains of the variables by means of crisp or fuzzy modalities. MCA does not explicitly exploit time, but the temporal dimension is introduced in this study with a basic solution that consists in representing each time sample (or time window) of a time-series by one statistical individual (a row of the table of analysis) and to perform MCA. Statistical individuals are represented by trajectories onto the factorial plane and their temporal evolution can then be exploited. Consequently, the proposed method can be useful for graphically follow the evolution of a given patient’s state by means of a synthetic representation that takes into account the most pertinent information in the data.

In the present study, it has been possible to evaluate and discuss the synthetic indices provided by the first two factorial axes of the MCA and to explain them according to functional and physiological points of view. The 58 records provided by the clinical protocol have been represented by trajectories on the first factorial plane of the MCA. The definition of an appropriate similarity measure has been discussed and has enabled the clustering of the trajectories into groups of trajectories with similar locations and evolutions on the factorial plane. The proposed distance has been validated on simulated datasets and enables the clustering of trajectories in the factorial plane according to their location and/or their evolution depending on the problem to be solved. Experiments on simulated datasets regarding the noise robustness also show the relevance of the association of fuzzy spatio-temporal coding and smoothed MCA. The proposed method is robust to white noise with SNR as low as 3 dB.

A database constituted of clinical observations from 41 patients has also been analyzed by using a data mining approach so as to characterize the data. Two areas in the factorial plane, corresponding to two large groups of patients, have been identified, the bottom-right quarter-plane being related with a health degradation and the bottom-left one with a stable clinical state. Most of the trajectories projected on this first factorial plane were correctly clustered according to their location and shape. Discussed individually, these clusters were efficient in grouping trajectories corresponding to similar patients’ clinical state. The present study has shown that patients undergoing an adverse event often present trajectories with abrupt variations. A detection of such phenomena would permit to identify patients with a critical evolution.

In the future, additional data would enable the identification and detection of typical evolutions related to health deterioration. Rules may be extracted from the location of modalities on the factorial plane, and enable the definition of thresholds on time-series to generate alarms associated with adverse events or health deterioration. These alarms could be sent from the pacemaker, via a data communication device (such as a mobile phone, PDA, etc.), to a telemonitoring center. These results are encouraging and may be useful for the definition of new selection criteria of candidate patients for CRT. Finally, the proposed methodology can be generalized to other monitoring problems.

References