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()
					

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:

In [16]:

						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).

In [17]:

						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.

In [18]:

						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í:

In [19]:

						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í.

In [20]:

						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:

In [21]:

						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.

In [22]:

						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.

Přidat názor

18. 5. 2022 14:31

misch

Těším se na další díl! Jsem totální laik, ale tohle mi přijde jako velice praktická ukázka. S omezením testovací sady dat se toho nedá moc dělat, ale i přesto (nebo právě proto?) jsem zvědavý na pokračování a na info jak moc se dá aktuální detekce zlepšit.

Přidat názor

Zdrojová data

Naivní klasifikátor

RandomForest klasifikátor

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

Plně propojené vrstvy

Konvoluční vrstvy

Sdílet

Jiří Raška

Nejčtenější články autora

Poslední názory

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

Zdrojová data

Naivní klasifikátor

RandomForest klasifikátor

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

Lineární klasifikátor

Plně propojené vrstvy

Konvoluční vrstvy

Sdílet

Jiří Raška

Nejčtenější články autora

Poslední názory

Dále u nás najdete

Jak českým firmám pomáhá digitální transformace?

Co dělat, když vám zavřou dodavatele IT služeb?

Prokletí slavného spisovatele: duševní poruchy a sebevraždy

Vyhněte se chybám a prodlužte životnost USB flash disku

SSD už také mají obrovskou kapacitu – až 122 TB

Čeští středoškoláci očekávají vypuštění své družice

Šlamastika kolem pozastavené výplaty důchodu

Svařák bez alkoholu je skvělá varianta pro těhotné či abstinenty

Pokud potřebujete půjčit, nemusíte prodávat bitcoiny. Stačí je zastavit

Chrome OS se zřejmě promění v Android

Souvisí zákaz Huawei i s nákupem stíhaček F-35?

Osvědčené způsoby, jak snížit vysoký krevní tlak

V aplikaci Záchranka je rozcestník péče pro nevyléčitelně nemocné

Stow vytváří regálové systémy pro Alzu, Lidl i Rohlík

Vánoční slevy operátorů? Daleko za očekáváním

Od pokeru do dílny: Vyrábí v ní nábytek z ohýbaného betonu

Přílišné uklízení škodí plicím stejně jako cigarety

Vyznejte se v kuřatech. Takhle se liší selské od venkovského

Lumbální punkce sice vypadá děsivě, ale nebolí

World of Warcraft slaví 20 let a stále se hraje