Klasifikace EKG křivek - výlet do světa neuronových sítí

13. 5. 2022 0:00 Jiří Raška

Tak jsem si jednou řekl, že by nebylo špatné se procvičit v návrhu a využítí neuronových sítí. A nejlépe se to dělá, pokud si člověk vybere nějaký reálný problém. Takto řečeno je to hodně široký pojem. Co by to ale mělo být? Chtěl jsem si vyzkoušet nějaký klasický klasifikační problém, tedy něco v oblasti supervised learning. Navíc jsem se chtěl podívat na to, jak si neuronová síť poradí s extrakcí vlastností ze surových dat. No a protože se mnoho let pohybuji v oblasti informačních systémů pro zdravotnictví, tak jsem hledal něco v této oblasti. Nakonec volba padla na klasifikaci EKG křivek.

Zdrojová data

Prvním a velice podstatným problémem je, kde vzít nějaká relevantní data pro vybraný úkol. A nejenom tak ledajaká data, ale musí být klasifikována do tříd. Vydal jsem se tedy hledat na server Kaggel, kde je velké množství datových balíčků určených pro soutěžení a hraní si v oblasti strojového učení. Navíc se zde můžete podívat a pochytit hodně podnětů od lidí, kteří již před vámi tento problém řešili.

Vybral jsem si tedy datový soubor ECG Heartbeat Categorization Dataset, který obsahuje dvě sady dat z projektů the MIT-BIH Arrhythmia Dataset a The PTB Diagnostic ECG Database. V dalším svém bádání budu používat data pouze z první sady.

Jedná se o křivky nasbírané z EKG vyšetření pacientů, které byly pro účely trénování sítí segmentovány do samostatných křivek pro jednotlivé srdeční údery. Dále byly lékaři klasifikovány do pěti tříd. Bližší infomrace k datům jistě naleznete na výše připojeném odkazu.

Nejdříve se tedy podívám na to, co se z těch dat dá vyčíst v jejich surové podobě:

In [1]:

						import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style('darkgrid')

np.random.seed(0)
					

Data mám stažena do lokálního adresáře ve formě CSV souborů, takže načtení do DataFrame je docela přímočaré:

In [2]:

						df_data = pd.read_csv('./data/ecg/mitbih_train.csv', header=None)
df_test = pd.read_csv('./data/ecg/mitbih_test.csv', header=None)
					

Oba DateFrame obsahují v jednom řádku hodnoty pro jednu křivku. Posledním sloupcem je číslo třídy, ve které byla křivka klasifikována. Takže pokud chci oddělit vstupní data od očekávaného výsledku, musím poslední sloupec vytáhnout zvláště.

In [3]:

						X_data, y_data = df_data.iloc[:, :-1].to_numpy(), df_data.iloc[:, -1].to_numpy()
X_test, y_test = df_test.iloc[:, :-1].to_numpy(), df_test.iloc[:, -1].to_numpy()

print(f"DATA shape: {X_data.shape}, {y_data.shape}")
print(f"TEST shape: {X_test.shape}, {y_test.shape}")

labels = ['Normal beat','Supraventricular premature beat','Premature ventricular contraction','Fusion of ventricular and normal beat','Unclassifiable beat']
					

DATA shape: (87554, 187), (87554,)
TEST shape: (21892, 187), (21892,)

Pokud se podívám na data křivek ze základní sady, pak je vidět, že jsou všechny normalizovány do rozsahu <0, 1>. Dá se předpokládat, že tak jsou upraveny i data z testovací sady (ověřoval jsem si to).

In [4]:

						pd.DataFrame(X_data).describe().T
					

Out[4]:

	count	mean	std	min	25%	50%	75%	max
0	87554.0	0.890360	0.240909	0.0	0.921922	0.991342	1.000000	1.0
1	87554.0	0.758160	0.221813	0.0	0.682486	0.826013	0.910506	1.0
2	87554.0	0.423972	0.227305	0.0	0.250969	0.429472	0.578767	1.0
3	87554.0	0.219104	0.206878	0.0	0.048458	0.166000	0.341727	1.0
4	87554.0	0.201127	0.177058	0.0	0.082329	0.147878	0.258993	1.0
…	…	…	…	…	…	…	…	…
182	87554.0	0.003681	0.037193	0.0	0.000000	0.000000	0.000000	1.0
183	87554.0	0.003471	0.036255	0.0	0.000000	0.000000	0.000000	1.0
184	87554.0	0.003221	0.034789	0.0	0.000000	0.000000	0.000000	1.0
185	87554.0	0.002945	0.032865	0.0	0.000000	0.000000	0.000000	1.0
186	87554.0	0.002807	0.031924	0.0	0.000000	0.000000	0.000000	1.0

187 rows × 8 columns

Dále je dobré se podívat na zastoupení vzorků v jednotlivých třídách.

Takto to vypadá v základní sadě:

In [5]:

						pd.Series(np.bincount(y_data.astype(int)))
					

Out[5]:

0    72471
1     2223
2     5788
3      641
4     6431
dtype: int64

A takto v sadě testovací:

In [6]:

						pd.Series(np.bincount(y_test.astype(int)))
					

Out[6]:

0    18118
1      556
2     1448
3      162
4     1608
dtype: int64

Rozložení vzorků v základní i testovací sadě je obdobné, takže to je dobrá zpráva. Horší je to, že vzorků ve třídě 0 je výrazně více naž vzorků v třídách ostatních, a to by mohla být komplikace.

Nyní se ještě podívám na nějaké vzorky křivek rozdělených dle tříd, abych měl bližší představu o jejich průběhu:

In [7]:

						np.random.seed(0)
sample_used = 10

for group_id, group_name in enumerate(labels):
    sns.relplot(data=pd.DataFrame(X_data[y_data == group_id]).sample(sample_used).T, kind='line', dashes=False, height=4, aspect=4).set(title=group_name)
plt.show()
					

Naivní klasifikátor

Na první pohled to bude vypadat jako nesmysl, ale jakou úspěšnost klasifikace bych dosáhl, kdybych prostě střílel naslepo? Ono to není vůbec špatné si hned na začátku udělat nějakou představu o úspěšnosti, kterou potřebuji dosáhnout, aby celé modelování mělo nějaký smysl.

Náhodné zařazení do tříd asi moc úspěšné nebude, ale co zařazení všech křivek do třídy 0?

In [8]:

						from sklearn.metrics import accuracy_score,confusion_matrix

y_pred = np.random.randint(low=y_test.min(), high=y_test.max()+1, size=y_test.shape)
print(f"Accuracy Score for random values: {accuracy_score(y_test, y_pred):.1%}")
print()

y_pred = np.zeros(y_test.shape)
print(f"Accuracy Score for all zeroes:    {accuracy_score(y_test, y_pred):.1%}")
print(f"Confusion Matrix: \n{confusion_matrix(y_test, y_pred)}")

sns.heatmap(confusion_matrix(y_test, y_pred, normalize='true'), annot=True, fmt='.1%', cmap='Blues')
plt.show()
					

Accuracy Score for random values: 20.0%

Accuracy Score for all zeroes:    82.8%
Confusion Matrix:
[[18118     0     0     0     0]
 [  556     0     0     0     0]
 [ 1448     0     0     0     0]
 [  162     0     0     0     0]
 [ 1608     0     0     0     0]]

Z výsledku je vidět, že jakýkoliv můj model bude mít smysl, pokud jeho úspěšnost na testovacích datech bude větší jak 83%.

RandomForest klasifikátor

Toto není klasifikátor z oblasti neuronových sítí, ale jen tak pro zajímavost jsem si chtěl vyzkoušet, jak úspěšný bude algoritmu RandomForest.

Nebudu nikterak ladit jeho parametry, pouze do něj pustím zdrojová data a uvidím, jak si povede na těch testovacích:

In [9]:

						from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier()
rf.fit(X_data, y_data)

y_pred = rf.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
conf_matrix = confusion_matrix(y_test, y_pred)

print(f"Accuracy Score:  {accuracy:.1%}")
print(f"Confusion Matrix: \n{conf_matrix}")

sns.heatmap(confusion_matrix(y_test, y_pred, normalize='true'), annot=True, fmt='.1%', cmap='Blues')
plt.show()
					

Accuracy Score:  97.4%
Confusion Matrix:
[[18102     5     9     0     2]
 [  219   334     2     0     1]
 [  151     0  1278    14     5]
 [   49     0    12   101     0]
 [   90     0     3     0  1515]]

Na první pokus to vůbec není špatné. Úspěšnost 97% je výrazně lepší než výsledky naivního klasifikátoru. Podle očekávání mám největší problém s klasifikací třídy 1 a 3 (to jsou třídy s nejmenším zastoupením vzorků v trénovací sade). Navíc je vidět velké ovlivnění třídou 0, kdy chybně klasifikované vzorky jsou převážně špatně klasifikované právě do této třídy.

Úprava zdrojových dat pro trénování neuronové sítě

Na tomto místě je dobré si uvědomit, jak vzorky dat vznikly. EKG křivky každého pacienta byly rozsekány na jednotlivé segmenty podle srdečního úderu a zařazeny do dat jako samostatné vzorky. Z toho vyplývá, že v datech mám hodně křivek od jednoho pacienta, takže mezi nimi bude jistě úzký vztah. Z toho vyplývá, že mezi daty v trénovací a validační sadě bude nějaký vztah. Udělat s tím moc nemůžu, jen je potřeba si být toho vědom.

Doporučení pro úpravu dat před jejich použitím pro trénování neuronových sítí říkají, že je dobré provést jejich standardizaci. Zkusím se tohoto doporučení držet. Výsledkem by měly být upravené sady zdrojových a testovacích dat s průměrem nula a směrodatnou odchylkou jedna.

In [10]:

						from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaler.fit(X_data)
X_data = scaler.transform(X_data)
X_test = scaler.transform(X_test)
					

Pro další bádání si ještě zdrojová data rozdělím na sadu pro trénování a sadu pro validaci, a to v poměru 80:20.

In [24]:

						from sklearn.model_selection import train_test_split

X_train, X_valid, y_train, y_valid = train_test_split(X_data, y_data, train_size=0.8)

A ještě budu potřeboval importovat knihovny TensorFlow a Keras:

In [12]:

					import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers

tf.random.set_seed(0)
				

Lineární klasifikátor

Jako první pokus zkusím úplně jednoduchý model s jednou vrstvou, která by měla fungovat jako lineární klasifikátor.

In [13]:

						def build_model1():
    model = keras.Sequential([
        layers.Dense(5, activation='softmax'),
    ])
    model.compile(optimizer="rmsprop", loss="sparse_categorical_crossentropy", metrics=["accuracy"])
    return model
					

Jako algoritmus pro trénování používám RMSProp, což se obvykle jeví jako dobrá varianta. Optimalizovaná funkce je SparseCategoricalCrossentropy, což je funkce vhodná pro klasifikace tříd zakódovaných do jedné hodnoty. Navíc budu při trénování sledovat ještě metriku, jak úspěšný jsem při rozpoznání třídy.

Vyzkouším, jak bude model úspěšný pro různé velikosti dávky pro gradient descent:

In [14]:

						results = pd.DataFrame()
metric = 'val_accuracy'
for batch_size in (32, 64, 128, 256, 512):
    print(f"Batch size: {batch_size}")
    history = build_model1().fit(X_train, y_train, epochs=10, batch_size=batch_size, validation_data=(X_valid, y_valid), verbose=0)
    results[f'{metric}_{batch_size}'] = history.history[metric]

sns.relplot(data=results, kind='line', height=4, aspect=4)
plt.show()
					

Batch size: 32
Batch size: 64
Batch size: 128
Batch size: 256
Batch size: 512

Vypadá to, že velikost batch_size=128 bude docela dobrý kompromis mezi rychlostí a přijatelným výsledkem.

Dále se zkusím podívat na průběh optimalizované funkce pro více epoch. Můžu vůbec takový model přetrénovat?

In [25]:

						model1 = build_model1()
history = model1.fit(X_train, y_train, epochs=100, batch_size=128, validation_data=(X_valid, y_valid), verbose=0)
results = pd.DataFrame(history.history)

sns.relplot(data=results[['loss', 'val_loss']][10:], kind='line', height=4, aspect=4)
plt.show()
					