Clustering de series temporales con Python

La proliferación y la ubicuidad de los datos con dependencia temporal en un amplio abanico de disciplinas genera un interés sustancial en el análisis y la extracción de series temporales. La agrupación o clustering es uno de los métodos de extracción de datos más populares, no sólo por su poder de exploración, sino también como paso previo al procesamiento de otras técnicas.

tslearn es una librería de Machine Learning de Python para series temporales que ofrece herramientas para el preprocesamiento y la extracción de características, así como modelos dedicados para clustering, clasificación y regresión. En esta rápida prueba de concepto, vamos a realizar una tarea de clustering de un dataset de una sola variable bien conocido.

Importando librerías

Para empezar, importamos las librerías que vamos a usar:

import tslearn
import numpy as np
import time
import matplotlib.pyplot as plt

Con TimeSeriesScalerMeanVariance procesaremos las time series para que tengan una media nula y una desviación típica igual a 1.

from tslearn.preprocessing import TimeSeriesScalerMeanVariance

Obteniendo un dataset

Afortunadamente, tslearn nos ofrece un gran número de series temporales con las que podemos empezar a trastear. Podemos acceder a los UCR/UEA time series datasets de la siguiente manera.

from tslearn.datasets import UCR_UEA_datasets

En el siguiente enlace se puede consultar la extensa lista de datasets de series temporales a la que podemos acceder. Como se puede ver, todos estos datasets están convenientemente etiquetados, por lo que podremos comprobar la bondad de nuestro modelos una vez que hayamos terminado de clusterizar las secuencias.

En este artículo vamos a trabajar con el dataset ECG5000. Este dataset contiene 20 horas de registros de electrocardiogramas (ECGs). En total se cuentan con 5,000 pares de ECGs y etiquetas. Fue publicado en «Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals. Circulation 101(23)«.

Es muy simple acceder a las particiones de entrenamiento y testeo del dataset que nos interesa.

X_test, y_test, X_train, y_train = UCR_UEA_datasets().load_dataset("ECG5000")

Las etiquetas del dataset tienen el siguiente significado:

• 1 = N – Normal
• 2 = SVEB – Supraventricular ectopic beats
• 3 = VEB – Ventricular ectopic beats
• 4 = Fusion beat
• 5 = Unknown beat

Entendiendo el dataset

Lo primero que hay que comprobar es la distribución de etiquetas en el dataset.

labels, counts = np.unique(y_train, return_counts=True)
_labels = ['Normal', 'SVEB', 'VEB', 'Fusion', 'Unknown']
plt.figure(figsize=(8, 6))
plt.bar(_labels, counts, align='center')
plt.gca().set_xticks(_labels)
plt.grid()
plt.ylabel('No. samples')
plt.xlabel('Beat type')
plt.show()

Como se puede observar, el dataset está muy desbalanceado. Hay muy pocas muestras de las clases y_train > 2 en comparación con los latidos normales o SVEB. Por ello, para simplificar el análisis, vamos a partir el dataset en latidos normales (y_train | y_test = 0) y latidos anómalos (y_train | y_test = 1)

y_test[y_test == 1] = 0
y_test[y_test > 1] = 1

y_train[y_train == 1] = 0
y_train[y_train > 1] = 1

El número de muestras pertenecientes a cada clase sigue estando desbalanceado, aunque no de manera tan extrema como antes. Esto nos permitirá clusterizar en dos grupos con mejores resultados.

labels, counts = np.unique(y_train, return_counts=True)
_labels = ['Normal', 'Anomalous']
plt.figure(figsize=(8, 6))
plt.bar(_labels, counts, align='center')
plt.gca().set_xticks(_labels)
plt.grid()
plt.ylabel('No. samples')
plt.xlabel('Beat type')
plt.show()

Preprocesado del dataset y visualización

Debemos escalar nuestras series temporales para poder obtener buenos resultados de clustering.

X_train = TimeSeriesScalerMeanVariance().fit_transform(X_train)
X_test = TimeSeriesScalerMeanVariance().fit_transform(X_test)

Como vemos, contamos con 4500 secuencias de 140 elementos de una sola variable para entrenar y 500 secuencias de igual longitud para testear.

print(X_train.shape)
print(X_test.shape)

(4500, 140, 1)
(500, 140, 1)

size = X_train.shape[1]
n_classes = len(set(y_train))

plt.figure(figsize=(10, 10))
for i, cl in enumerate(set(y_train)):
    plt.subplot(n_classes, 1, i + 1)
    for ts in X_train[y_train == cl]:
        if cl == 0:
            plt.plot(ts.ravel(), color="green" , alpha=.15)
        else:
            plt.plot(ts.ravel(), color="red" , alpha=.15)
    plt.xlim(0, size - 1)
    plt.grid()
plt.suptitle("Training time series")
plt.show()

plt.figure(figsize=(10, 10))
for i, cl in enumerate(set(y_train)):
    plt.subplot(n_classes, 1, i + 1)
    for ts in X_test[y_test == cl]:
        if cl == 0:
            plt.plot(ts.ravel(), color="green" , alpha=.15)
        else:
            plt.plot(ts.ravel(), color="red" , alpha=.15)
    plt.xlim(0, size - 1)
    plt.grid()
plt.suptitle("Test dataset time series")
plt.show()

A simple vista, únicamente en base a la forma, es sencillo distinguir entre los latidos de un individuo sano y los de un paciente con una patología (sin especificar).

Clustering de secuencias

Vamos a probar dos algoritmos de clustering: K-Means y K-Shape. Para el caso de K-Means, vamos a centrarnos únicamente en la métrica «euclidean». Hay otras dos alternativas para este hiperparámetro: «dtw» y «softdtw». Se deja como ejercicio al lector, si éste lo desea. Aviso: los tiempos de ejecución en estos dos casos son ostensiblemente más largos.

from tslearn.clustering import TimeSeriesKMeans, KShape, silhouette_score

Time Series K-Means para dos clusters

start_time = time.time()

km_euc = TimeSeriesKMeans(n_clusters=2, 
                          verbose=2, 
                          n_init=10, 
                          metric="euclidean").fit(X_train)
labels_euc = km_euc.labels_

print("--- %s seconds ---" % (time.time() - start_time))

...
--- 9.683111190795898 seconds ---

Time Series K-Shape para dos clusters

start_time = time.time()

ks = KShape(n_clusters=2, 
            verbose=True,
            n_init=10).fit(X_train)
labels_kshape = ks.labels_

print("--- %s seconds ---" % (time.time() - start_time))

...
--- 91.33036160469055 seconds ---

Las silhouette scores de ambas soluciones son muy similares:

silhouette_score(X_train, labels_euc, metric="euclidean") 

0.45911902366303975

silhouette_score(X_train, labels_kshape, metric="euclidean") 

0.4550914094152607

Resultados de los agrupamientos

Dado que contamos con las etiquetas del dataset, podemos comprobar cómo de bien hemos agrupado los electrocardiogramas entre pacientes sanos o con patologías.

from sklearn.metrics import accuracy_score

print('Training dataset accuracy for TimeSeriesKMeans clustering with euclidean metric: ', accuracy_score(y_train, labels_euc))

print('Training dataset accuracy for KShape clustering: ', accuracy_score(y_train, labels_kshape))

Training dataset accuracy for TimeSeriesKMeans clustering with euclidean metric:  0.9411
Training dataset accuracy for KShape clustering:  0.9391

print('Test dataset accuracy for TimeSeriesKMeans clustering with euclidean metric: ', accuracy_score(y_test, km_euc.predict(X_test)))

print('Test dataset accuracy for KShape clustering: ', accuracy_score(y_test, ks.predict(X_test)))

Test dataset accuracy for TimeSeriesKMeans clustering with euclidean metric:  0.958
Test dataset accuracy for KShape clustering:  0.954

La clasificación parece satisfactoria. Como vemos, tenemos una precisión cercana al 96% en el test dataset. Vamos a analizar, con matrices de confusión, los errores de Tipo I y Tipo II en los que nuestros modelos de clustering está incurriendo. Estos conceptos ya se trataron en los siguientes posts: Video Post #11: Tipos de error en Machine Learning y Tipos de error en Machine Learning, ¿los conoces?.

from sklearn.metrics import confusion_matrix

confusion_matrix(y_train, labels_euc)

	Normal	Anomalous
Normal	2563	64
Anomalous	201	1672

confusion_matrix(y_train, labels_kshape)

	Normal	Anomalous
Normal	2545	82
Anomalous	192	1681

confusion_matrix(y_test, km_euc.predict(X_test))

	Normal	Anomalous
Normal	283	9
Anomalous	12	196

confusion_matrix(y_test, ks.predict(X_test))

	Normal	Anomalous
Normal	283	11
Anomalous	12	196

El desempeño de ambos algoritmos es muy parecido. Dado que estamos trabajando con datos clínicos de pacientes, lo que deseamos es tener el menor error de tipo I posible. Es decir, debemos minimizar la posibilidad de que a una persona con una patología cardiaca no se le diagnostique correctamente, considerándolo una persona sana. Para ambas aproximaciones, esto sólo ocurre en 12 muestras en el test dataset.

Representando los clusters y sus centroides

Con las siguientes líneas de código podremos representar los distintos grupos de series temporales junto con los centroides que los definen. Estos centroides se consiguen con model.cluster_centers_.

def plot_groups(model):

    if model == 'kmeans-euclidean':
        model = km_euc
    elif model == 'kshape':
        model = ks
    else:
        sys.exit('Please, provide a valid model string.')
        
    labels = model.labels_

    plt.figure(figsize=(10, 10))
    for yi in range(2):
        plt.subplot(2, 1, 1 + yi)
        for xx in X_train[labels == yi]:
            plt.plot(xx.ravel(), "k-", alpha=.15)
        plt.plot(model.cluster_centers_[yi].ravel(), "r-")
        plt.xlim(0, size-1)
        
        if yi == 0:
            plt.title("Normal")
        else:
            plt.title("Anomalous")
        plt.grid()

    plt.tight_layout()
    plt.show()

plot_groups(model='kmeans-euclidean')

def compare_cluster_centers():

    plt.figure(figsize=(10, 10))
    for yi in range(2):
        plt.subplot(2, 1, 1 + yi)
        plt.plot(km_euc.cluster_centers_[yi].ravel(), "-", alpha=1., label='K-Means euclidean')
        #plt.plot(km_softdtw.cluster_centers_[yi].ravel(), "-", alpha=0.5, label='K-Means softdtw')
        plt.plot(ks.cluster_centers_[yi].ravel(), "-", alpha=.5, label='K-Shape')
        if yi == 0:
            plt.title("Normal")
        else:
            plt.title("Anomalous")
        
        plt.xlim(0, size-1)
        plt.grid()
    plt.legend()

    plt.tight_layout()
    plt.show()

Conclusiones

Los centroides de los clusters son bastante similares comparando entre los dos algoritmos usados, tal y como se preveía dadas las tasas de acierto y error mostradas por las matrices de confusión.

def compare_cluster_centers():

    plt.figure(figsize=(10, 10))
    for yi in range(2):
        plt.subplot(2, 1, 1 + yi)
        plt.plot(km_euc.cluster_centers_[yi].ravel(), "-", alpha=1., label='K-Means euclidean')
        plt.plot(ks.cluster_centers_[yi].ravel(), "-", alpha=.5, label='K-Shape')
        if yi == 0:
            plt.title("Normal")
        else:
            plt.title("Anomalous")
        
        plt.xlim(0, size-1)
        plt.grid()
    plt.legend()

    plt.tight_layout()
    plt.show()

compare_cluster_centers_classes()

def compare_cluster_centers_classes():

    plt.figure(figsize=(10, 6))
    plt.plot(km_euc.cluster_centers_[0].ravel(), "-", c='green', alpha=.75, label='Normal')
    plt.plot(km_euc.cluster_centers_[1].ravel(), "-", c='red', alpha=.75, label='Anomalous')
    plt.xlim(0, size-1)
    plt.title("Cluster Centers comparison per classes")
    plt.tight_layout()
    plt.legend()
    plt.grid()
    plt.show()

compare_cluster_centers_classes()

Resulta muy interesante ver cómo con unas pocas líneas de código en Python, se puede obtener un perfil promedio de latido que nos permite diferenciar entre un ritmo cardíaco normal y un individuo con una patología.

Otra opción alternativa a la clusterización de series temporales de este tipo podría abordarse a través de una reducción de dimensionalidad con PCA. Por ejemplo, se podría convertir el training dataset en un pandas dataframe de 4500 filas y 140 columnas y, sobre ese array, realizar un Principal Component Analysis. Con 2 componentes mantendríamos un 65% de varianza y con 3 componentes subiríamos hasta el 75%. Con esta aproximación y también de manera muy rápida, podríamos encontrar patrones comunes entre secuencias muy bien agrupados.

Os invitamos a probar a usar tslearn en problemas de clasificación, clustering y regresión de series temporales con éste y otros de los muchos datasets disponibles.

Para mantenerte al día con LUCA visita nuestra página web,  suscríbete a LUCA Data Speaks o síguenos en Twitter, LinkedIn y YouTube.