Klasifikace EKG křivek - výlet do světa neuronových sítí
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ě:
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é:
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ě.
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).
pd.DataFrame(X_data).describe().T
| 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ě:
pd.Series(np.bincount(y_data.astype(int)))
0 72471 1 2223 2 5788 3 641 4 6431 dtype: int64
A takto v sadě testovací:
pd.Series(np.bincount(y_test.astype(int)))
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:
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?
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:
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.
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.
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:
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.
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:
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?
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()
Není to tedy ten klasický průběh se zhoršováním výsledků optimalizované funkce na validační sadě, ale jeho stagnace kolem 60 epoch je vidět.
Zkusím se tedy podívat na výsledky tohoto modelu pro celou datovou sadu v porovnání s testovaci:
model1 = build_model1()
model1.fit(X_data, y_data, epochs=60, batch_size=128, verbose=0)
y_pred = model1.predict(X_test).argmax(axis=1)
print(f"Accuracy Score: {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: 91.1% Confusion Matrix: [[17682 33 370 4 29] [ 285 240 31 0 0] [ 829 14 555 23 27] [ 108 0 8 46 0] [ 169 1 20 0 1418]]
Výsledek tedy zatím není nic moc, ale zlepšení proti naivnímu modelu to je.
Plně propojené vrstvy
Zkusím tedy do modelu doplnit dvě plně propojené vrstvy s jednotkami aktivovanými funkcí relu (rectified linear unit).
def build_model2(units=32):
model = keras.Sequential([
layers.Dense(units, activation='relu'),
layers.Dense(units, activation='relu'),
layers.Dense(5, activation='softmax'),
])
model.compile(optimizer="rmsprop", loss="sparse_categorical_crossentropy", metrics=["accuracy"])
return model
model2 = build_model2()
model2(X_train)
model2.summary()
Model: "sequential_7" _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= dense_7 (Dense) (70043, 32) 6016 dense_8 (Dense) (70043, 32) 1056 dense_9 (Dense) (70043, 5) 165 ================================================================= Total params: 7,237 Trainable params: 7,237 Non-trainable params: 0 _________________________________________________________________
Ve funkci mám jako parametr počet jednotek v jedné vrstvě, protože bych chtěl vyzkoušet, jaký vliv na výsledek bude mít jejich počet v modelu.
results = pd.DataFrame()
for units in (8, 16, 32, 64, 128, 256):
print(f"UNITS: {units}")
history = build_model2(units).fit(X_train, y_train, epochs=20, batch_size=128, validation_data=(X_valid, y_valid), verbose=0)
results[f'val_loss_{units}'] = history.history['val_loss']
sns.relplot(data=results, kind='line', height=4, aspect=4)
plt.show()
UNITS: 8 UNITS: 16 UNITS: 32 UNITS: 64 UNITS: 128 UNITS: 256
V grafu jsou hodnoty optimalizované funkce pro validační sadu dat. Jako nejpřínosnější se jeví počet jednotek 64. Více jednotek již nemá viditelný vliv na výsledek a zbytečně komlikuje model.
Zkusím tedy model aplikovat na celou zdrojovou sadu dat s ověřením výsledku na sadě testovací:
model2 = build_model2(64)
model2.fit(X_data, y_data, epochs=60, batch_size=128, verbose=0)
y_pred = model2.predict(X_test).argmax(axis=1)
print(f"Accuracy Score: {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: 97.5% Confusion Matrix: [[17872 140 48 35 23] [ 118 427 8 2 1] [ 52 8 1356 28 4] [ 18 0 20 123 1] [ 18 6 9 2 1573]]
Konvoluční vrstvy
Jedním z důvodů, proč jsem si pro bádání vybral klasifikaci EKG křivek, byla práce se surovými daty, ze kterých je potřeba vydolovat nějaké společné vlastnosti. K tomu by mně měly dopomoci konvoluční vrstvy, které předřadím před ty plně propojené.
Do modelu tedy přidám další tři vrstvy, které budou postupně redukovat dimenzi surových dat a zároveň vytvářet množinu vlastností.
def build_model3():
model = keras.Sequential([
layers.Conv1D(16, 3, activation='relu', input_shape=(X_data.shape[-1], 1)),
layers.MaxPooling1D(2),
layers.Conv1D(32, 3, activation='relu'),
layers.MaxPooling1D(2),
layers.Conv1D(64, 3, activation='relu'),
layers.Dropout(0.5),
layers.Flatten(),
layers.Dense(64, activation='relu'),
layers.Dense(64, activation='relu'),
layers.Dense(5, activation='softmax'),
])
model.compile(optimizer="rmsprop", loss="sparse_categorical_crossentropy", metrics=["accuracy"])
return model
model3 = build_model3()
model3.summary()
Model: "sequential_15" _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= conv1d (Conv1D) (None, 185, 16) 64 max_pooling1d (MaxPooling1D (None, 92, 16) 0 ) conv1d_1 (Conv1D) (None, 90, 32) 1568 max_pooling1d_1 (MaxPooling (None, 45, 32) 0 1D) conv1d_2 (Conv1D) (None, 43, 64) 6208 dropout (Dropout) (None, 43, 64) 0 flatten (Flatten) (None, 2752) 0 dense_31 (Dense) (None, 64) 176192 dense_32 (Dense) (None, 64) 4160 dense_33 (Dense) (None, 5) 325 ================================================================= Total params: 188,517 Trainable params: 188,517 Non-trainable params: 0 _________________________________________________________________
Na souhrnu modelu je vidět, jak jednotlivé konvoluční vrstvy redukují dimenzi vzorku a vytváří množinu vlastností.
Podívám se, jak bude model fungovat na trénovacích a validačních datech:
history = build_model3().fit(X_train, y_train, epochs=40, 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)
sns.relplot(data=results[['accuracy', 'val_accuracy']][10:], kind='line', height=4, aspect=4)
plt.show()
Vypadá to, že pro větší množství epoch začíná model degradovat.
Proto jsem pro další trénování použil dvě omezení ve formě callback. Za prvé se jedná o dřívější ukončení trénování v situaci, kdy se sledovaná velična (v tomto případě ztrátová funkce) opakovaně zhoršuje. Dále jsem doplnil ukládání modelu do souboru v případě, že jeho výsledek je lepší než všech předchozích. Pro predikci se dále použije ten nejlepší model, který se podařilo v průběhu trénování najít.
MODEL_PATH = 'models/ECG-model3.keras'
callbacks = [
keras.callbacks.EarlyStopping(monitor='accuracy', patience=5),
keras.callbacks.ModelCheckpoint(filepath=MODEL_PATH, save_best_only=True, monitor='loss')
]
model3 = build_model3()
model3.fit(X_data, y_data, epochs=100, batch_size=128, callbacks=callbacks, verbose=0)
y_pred = keras.models.load_model(MODEL_PATH).predict(X_test).argmax(axis=1)
print(f"Accuracy Score: {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: 98.4% Confusion Matrix: [[18075 15 10 3 15] [ 146 402 6 1 1] [ 69 4 1350 21 4] [ 29 0 10 123 0] [ 21 1 2 0 1584]]
Postupně jsem se dopracoval k zatím nejlepšímu výsledku, i když problém s rozpoznáním třídy 1 a 3 stále přetrvává.
A to je dnes vše.
