Ccmmutty logo
Commutty IT
0 pv57 min read

08 実践編: ディープラーニングを使ったモニタリングデータの時系列解析

https://picsum.photos/seed/e605ab16a7b147c58e99d432f2670c2f/1200/630
健康意識の高まりや運動人口の増加に伴って,活動量計などのウェアラブルデバイスが普及し始めています.センサーデバイスから心拍数などの情報を取得することで,リアルタイムに健康状態をモニタリングできる可能性があることから,近年ではヘルスケア領域での活用事例も増えてきています.2018年2月には,Cardiogram社とカリフォルニア大学が共同研究の成果を発表し,心拍数データに対してDeep Learningを適用することで,高精度に糖尿病予備群を予測可能であることを報告し,大きな注目を集めました.また,Apple Watch Series 4には心電図作成の機能が搭載されるなど,センサーデバイスも進歩を続け,より精緻な情報が取得できるようになってきています.こうした背景において,モニタリングデータを収集・解析し,健康管理に繋げていく取り組みは今後益々盛んになっていくものと考えられます.
本章では,心電図の信号波形データを対象として,不整脈を検出する問題に取り組みます.

環境構築

本章では, 下記のライブラリを利用します.
  • Cupy
  • Chainer
  • Scipy
  • Matplotlib
  • Seaborn
  • Pandas
  • WFDB
  • Scikit-learn
  • Imbalanced-learn
以下のセルを実行 (Shift + Enter) して,必要なものをインストールして下さい.
python
!apt -y -q install tree
!pip install wfdb==2.2.1 scikit-learn==0.20.1 imbalanced-learn==0.4.3
インストールが完了したら以下のセルを実行して,各ライブラリのインポート,及びバージョン確認を行って下さい.
python
import os
import random
import numpy as np
import chainer
import scipy
import pandas as pd
import matplotlib
import seaborn as sn
import wfdb
import sklearn
import imblearn

chainer.print_runtime_info()
print("Scipy: ", scipy.__version__)
print("Pandas: ", pd.__version__)
print("Matplotlib: ", matplotlib.__version__)
print("Seaborn: ", sn.__version__)
print("WFDB: ", wfdb.__version__)
print("Scikit-learn: ", sklearn.__version__)
print("Imbalanced-learn: ", imblearn.__version__)
Platform: Linux-4.14.65+-x86_64-with-Ubuntu-18.04-bionic Chainer: 5.1.0 NumPy: 1.14.6 CuPy: CuPy Version : 5.1.0 CUDA Root : /usr/local/cuda CUDA Build Version : 9020 CUDA Driver Version : 9020 CUDA Runtime Version : 9020 cuDNN Build Version : 7301 cuDNN Version : 7301 NCCL Build Version : 2307 iDeep: 2.0.0.post3 Scipy: 1.1.0 Pandas: 0.22.0 Matplotlib: 2.1.2 Seaborn: 0.7.1 WFDB: 2.2.1 Scikit-learn: 0.20.1 Imbalanced-learn: 0.4.3
また,本章の実行結果を再現可能なように,4章 (4.2.4.4) で紹介した乱数シードの固定を行います.
(以降の計算を行う上で必ず必要となる設定ではありません.)
python
def reset_seed(seed=42):
    random.seed(seed)
    np.random.seed(seed)
    if chainer.cuda.available:
        chainer.cuda.cupy.random.seed(seed)

reset_seed(42)

心電図(ECG)と不整脈診断について

**心電図(Electrocardiogram, ECG)**は,心筋を協調的,律動的に動かすために刺激伝導系で伝達される電気信号を体表から記録したものであり,心電図検査は不整脈や虚血性心疾患の診断において広く用いられる検査です[1, 2].
標準的な心電図は,手足からとる心電図(四肢誘導)として,双極誘導(),及び単極誘導(aVRaV_RaVLaV_LaVFaV_F)の6誘導,胸部からとる心電図(胸部誘導)として,V1V_1V2V_2V3V_3V4V_4V5V_5V6V_6の6誘導,計12誘導から成ります.このうち,特に不整脈のスクリーニングを行う際には,誘導とV1V_1誘導に注目して診断が行われるのが一般的とされています.
心臓が正常な状態では,ECGにおいては規則的な波形が観測され,これを**正常洞調律 (Normal sinus rhythm, NSR)**といいます.
具体的には,以下の3つの主要な波形で構成されており,
  1. P波:心房の脱分極(心房の興奮)
  2. QRS波:心室の脱分極(心室の興奮)
  3. T波:心室の再分極(心室興奮の収まり)
の順番で,下図のような波形が観測されます.
正常心電図の概略図
([1]より引用)
こうした規則的な波形に乱れが生じ,調律に異常があると判断された場合,不整脈などの疑いがあるため,診断が行われることになります.

使用するデータセット

ここでは,ECGの公開データとして有名なMIT-BIH Arrhythmia Database (mitdb)を使用します.
47名の患者から収集した48レコードが登録されており,各レコードファイルには約30分間の2誘導(IIIIV1V_1)のシグナルデータが格納されています.また,各R波のピーク位置に対してアノテーションが付与されています.(データとアノテーションの詳細についてはこちらを御覧ください.)
データベースはPhysioNetによって管理されており,ダウンロードや読み込み用のPythonパッケージが提供されているので,今回はそちらを利用してデータを入手します.
python
dataset_root = './dataset'
download_dir = os.path.join(dataset_root, 'download')
まずはmitdbデータベースをダウンロードしましょう. ※実行時にエラーが出た場合は,再度実行して下さい.
python
wfdb.dl_database('mitdb', dl_dir=download_dir)
Created local base download directory: ./dataset/download Downloading files... Finished downloading files
無事ダウンロードが完了すると, Finished downloading files というメッセージが表示されます.
ファイル一覧を確認してみましょう.
python
print(sorted(os.listdir(download_dir)))
['100.atr', '100.dat', '100.hea', '101.atr', '101.dat', '101.hea', '102.atr', '102.dat', '102.hea', '103.atr', '103.dat', '103.hea', '104.atr', '104.dat', '104.hea', '105.atr', '105.dat', '105.hea', '106.atr', '106.dat', '106.hea', '107.atr', '107.dat', '107.hea', '108.atr', '108.dat', '108.hea', '109.atr', '109.dat', '109.hea', '111.atr', '111.dat', '111.hea', '112.atr', '112.dat', '112.hea', '113.atr', '113.dat', '113.hea', '114.atr', '114.dat', '114.hea', '115.atr', '115.dat', '115.hea', '116.atr', '116.dat', '116.hea', '117.atr', '117.dat', '117.hea', '118.atr', '118.dat', '118.hea', '119.atr', '119.dat', '119.hea', '121.atr', '121.dat', '121.hea', '122.atr', '122.dat', '122.hea', '123.atr', '123.dat', '123.hea', '124.atr', '124.dat', '124.hea', '200.atr', '200.dat', '200.hea', '201.atr', '201.dat', '201.hea', '202.atr', '202.dat', '202.hea', '203.atr', '203.dat', '203.hea', '205.atr', '205.dat', '205.hea', '207.atr', '207.dat', '207.hea', '208.atr', '208.dat', '208.hea', '209.atr', '209.dat', '209.hea', '210.atr', '210.dat', '210.hea', '212.atr', '212.dat', '212.hea', '213.atr', '213.dat', '213.hea', '214.atr', '214.dat', '214.hea', '215.atr', '215.dat', '215.hea', '217.atr', '217.dat', '217.hea', '219.atr', '219.dat', '219.hea', '220.atr', '220.dat', '220.hea', '221.atr', '221.dat', '221.hea', '222.atr', '222.dat', '222.hea', '223.atr', '223.dat', '223.hea', '228.atr', '228.dat', '228.hea', '230.atr', '230.dat', '230.hea', '231.atr', '231.dat', '231.hea', '232.atr', '232.dat', '232.hea', '233.atr', '233.dat', '233.hea', '234.atr', '234.dat', '234.hea']
ファイル名の数字はレコードIDを表しています.各レコードには3種類のファイルがあり,
  • .dat : シグナル(バイナリ形式)
  • .atr : アノテーション(バイナリ形式)
  • .hea : ヘッダ(バイナリファイルの読み込みに必要)
となっています.

データ前処理

ダウンロードしたファイルを読み込み,機械学習モデルへの入力形式に変換するデータ前処理について説明します.
本節では,以下の手順で前処理を行います.
  1. レコードIDを事前に 学習用 / 評価用 に分割
    • 48レコードのうち,
      • ID =(102, 104, 107, 217)のシグナルはペースメーカーの拍動が含まれるため除外します.
      • ID = 114のシグナルは波形の一部が反転しているため,今回は除外します.
      • ID = (201, 202)は同一の患者から得られたデータのため,202を除外します.
    • 上記を除く計42レコードを,学習用とテスト用に分割します(分割方法は[3]を参考).
  2. シグナルファイル (.dat) の読み込み
    • 誘導シグナルとV1V_1誘導シグナルが格納されていますが,今回は誘導のみ利用します.
    • サンプリング周波数は360 Hz なので,1秒間に360回のペースで数値が記録されていることになります.
  3. アノテーションファイル (.atr) の読み込み
    • 各R波ピークの位置 (positions) と,そのラベル (symbols) を取得します.
  4. シグナルの正規化
    • 平均0,分散1になるように変換を行います.
  5. シグナルの分割 (segmentation)
    • 各R波ピークを中心として2秒間(前後1秒ずつ)の断片を切り出していきます.
  6. 分割シグナルへのラベル付与
    • 各R波ピークに付与されているラベルを,下表(※)に従って集約し,今回の解析では正常拍動 (Normal),及び心室異所性拍動 (VEB) に対応するラベルが付与されている分割シグナルのみ学習・評価に利用します.
※ Association for the Advancement of Medical Instrumentation (AAMI)が推奨している基準([3])で,5種類に大別して整理されています.
AAMIの分類基準
([4]より引用)
まずは以下のセルを実行して,データ前処理クラスを定義しましょう.
前処理クラス内では,以下のメンバ関数を定義しています.
  • __init__() (コンストラクタ) : 変数の初期化,学習用とテスト用への分割ルール,利用するラベルの集約ルール
  • _load_data() : シグナル,及びアノテーションの読み込み
  • _normalize_signal() : methodオプションに応じてシグナルをスケーリング
  • _segment_data() : 読み込んだシグナルとアノテーションを,一定幅(window_size)で切り出し
  • preprocess_dataset() : 学習データ,テストデータを作成
  • _preprocess_dataset_core() : preprocess_datataset()内で呼ばれるメインの処理.
python
class BaseECGDatasetPreprocessor(object):

    def __init__(
            self,
            dataset_root,
            window_size=720,  # 2 seconds
    ):
        self.dataset_root = dataset_root
        self.download_dir = os.path.join(self.dataset_root, 'download')
        self.window_size = window_size
        self.sample_rate = 360.
        # split list
        self.train_record_list = [
            '101', '106', '108', '109', '112', '115', '116', '118', '119', '122',
            '124', '201', '203', '205', '207', '208', '209', '215', '220', '223', '230'
        ]
        self.test_record_list = [
            '100', '103', '105', '111', '113', '117', '121', '123', '200', '210',
            '212', '213', '214', '219', '221', '222', '228', '231', '232', '233', '234'
        ]
        # annotation
        self.labels = ['N', 'V']
        self.valid_symbols = ['N', 'L', 'R', 'e', 'j', 'V', 'E']
        self.label_map = {
            'N': 'N', 'L': 'N', 'R': 'N', 'e': 'N', 'j': 'N',
            'V': 'V', 'E': 'V'
        }

    def _load_data(
            self,
            base_record,
            channel=0  # [0, 1]
    ):
        record_name = os.path.join(self.download_dir, str(base_record))
        # read dat file
        signals, fields = wfdb.rdsamp(record_name)
        assert fields['fs'] == self.sample_rate
        # read annotation file
        annotation = wfdb.rdann(record_name, 'atr')
        symbols = annotation.symbol
        positions = annotation.sample
        return signals[:, channel], symbols, positions

    def _normalize_signal(
            self,
            signal,
            method='std'
    ):
        if method == 'minmax':
            # Min-Max scaling
            min_val = np.min(signal)
            max_val = np.max(signal)
            return (signal - min_val) / (max_val - min_val)
        elif method == 'std':
            # Zero mean and unit variance
            signal = (signal - np.mean(signal)) / np.std(signal)
            return signal
        else:
            raise ValueError("Invalid method: {}".format(method))

    def _segment_data(
            self,
            signal,
            symbols,
            positions
    ):
        X = []
        y = []
        sig_len = len(signal)
        for i in range(len(symbols)):
            start = positions[i] - self.window_size // 2
            end = positions[i] + self.window_size // 2
            if symbols[i] in self.valid_symbols and start >= 0 and end <= sig_len:
                segment = signal[start:end]
                assert len(segment) == self.window_size, "Invalid length"
                X.append(segment)
                y.append(self.labels.index(self.label_map[symbols[i]]))
        return np.array(X), np.array(y)

    def preprocess_dataset(
            self,
            normalize=True
    ):
        # preprocess training dataset
        self._preprocess_dataset_core(self.train_record_list, "train", normalize)
        # preprocess test dataset
        self._preprocess_dataset_core(self.test_record_list, "test", normalize)

    def _preprocess_dataset_core(
            self,
            record_list,
            mode="train",
            normalize=True
    ):
        Xs, ys = [], []
        save_dir = os.path.join(self.dataset_root, 'preprocessed', mode)
        for i in range(len(record_list)):
            signal, symbols, positions = self._load_data(record_list[i])
            if normalize:
                signal = self._normalize_signal(signal)
            X, y = self._segment_data(signal, symbols, positions)
            Xs.append(X)
            ys.append(y)
        os.makedirs(save_dir, exist_ok=True)
        np.save(os.path.join(save_dir, "X.npy"), np.vstack(Xs))
        np.save(os.path.join(save_dir, "y.npy"), np.concatenate(ys))
データ保存先のrootディレクトリ(dataset_root)を指定し, preprocess_dataset() を実行することで,前処理後のデータがNumpy Array形式で所定の場所に保存されます.
python
BaseECGDatasetPreprocessor(dataset_root).preprocess_dataset()
実行後,以下のファイルが保存されていることを確認しましょう.
  • train/X.npy : 学習用シグナル
  • train/y.npy : 学習用ラベル
  • test/X.npy : 評価用シグナル
  • test/y.npy : 評価用ラベル
python
!tree ./dataset/preprocessed
./dataset/preprocessed ├── test │   ├── X.npy │   └── y.npy └── train ├── X.npy └── y.npy
2 directories, 4 files
次に,保存したファイルを読み込み,中身を確認してみましょう.
python
X_train = np.load(os.path.join(dataset_root, 'preprocessed', 'train', 'X.npy'))
y_train = np.load(os.path.join(dataset_root, 'preprocessed', 'train', 'y.npy'))
X_test = np.load(os.path.join(dataset_root, 'preprocessed', 'test', 'X.npy'))
y_test = np.load(os.path.join(dataset_root, 'preprocessed', 'test', 'y.npy'))
データセットのサンプル数はそれぞれ以下の通りです.
  • 学習用 : 47738個
  • 評価用 : 45349個
各シグナルデータは,2 (sec) * 360 (Hz) = 720次元ベクトルとして表現されています.
python
print("X_train.shape = ", X_train.shape, " \t y_train.shape = ", y_train.shape)
print("X_test.shape = ", X_test.shape, " \t y_test.shape = ", y_test.shape)
X_train.shape = (47738, 720) y_train.shape = (47738,) X_test.shape = (45349, 720) y_test.shape = (45349,)
各ラベルはインデックスで表現されており,
  • 0 : 正常拍動 (Normal)
  • 1 : 心室異所性拍動 (VEB)
となっています.
学習用データセットに含まれている各ラベル毎のサンプル数をカウントしてみましょう.
python
uniq_train, counts_train = np.unique(y_train, return_counts=True)
print("y_train count each labels: ", dict(zip(uniq_train, counts_train)))
y_train count each labels: {0: 43995, 1: 3743}
評価用データについても同様にラベル毎のサンプル数をカウントします.
python
uniq_test, counts_test = np.unique(y_test, return_counts=True)
print("y_test count each labels: ", dict(zip(uniq_test, counts_test)))
y_test count each labels: {0: 42149, 1: 3200}
学習用データ,評価用データ共に,VEBサンプルは全体の10%未満であり,大多数は正常拍動サンプルであることが分かります.
次に,正常拍動,及びVEBのシグナルデータを可視化してみましょう.
python
%matplotlib inline
import matplotlib.pyplot as plt
正常拍動の例を図示したものが以下になります.
P波 - QRS波 - T波が規則的に出現していることが確認できます.
python
idx_n = np.where(y_train == 0)[0]
plt.plot(X_train[idx_n[0]])
[<matplotlib.lines.Line2D at 0x7f04055b3320>]
png
一方でVEBの波形は規則性が乱れ,R波ピークの形状やピーク間距離も正常例とは異なる性質を示していることが分かります.
python
idx_s = np.where(y_train == 1)[0]
plt.plot(X_train[idx_s[0]])
[<matplotlib.lines.Line2D at 0x7f040354c4a8>]
png
本章の目的は,ECGシグナル特徴をうまく捉え,新たな波形サンプルに対しても高精度に正常/異常を予測するモデルを構築することです.
次節では,深層学習を利用したモデル構築について説明していきます.

深層学習を用いた時系列データ解析

学習

まず,前節で準備した前処理済みデータをChainerで読み込むためのデータセットクラスを定義します.
python
class ECGDataset(chainer.dataset.DatasetMixin):

    def __init__(
            self,
            path
    ):
        if os.path.isfile(os.path.join(path, 'X.npy')):
            self.X = np.load(os.path.join(path, 'X.npy'))
        else:
            raise FileNotFoundError("{}/X.npy not found.".format(path))
        if os.path.isfile(os.path.join(path, 'y.npy')):
            self.y = np.load(os.path.join(path, 'y.npy'))
        else:
            raise FileNotFoundError("{}/y.npy not found.".format(path))

    def __len__(self):
        return len(self.X)

    def get_example(self, i):
        return self.X[None, i].astype(np.float32), self.y[i]
続いて,学習(+予測)に利用するネットワーク構造を定義します.
今回は,画像認識タスクで有名な,CNNベースのResNet34と同様のネットワーク構造を利用します.[5]
ただし,入力シグナルは1次元配列であることから,ここでは画像解析等で一般的に利用される2D Convolutionではなく,前章の遺伝子解析と同様,1D Convolutionを利用します.
python
import chainer.functions as F
import chainer.links as L
from chainer import reporter
from chainer import Variable
      
    
class BaseBlock(chainer.Chain):

    def __init__(
            self,
            channels,
            stride=1,
            dilate=1
    ):
        self.stride = stride
        super(BaseBlock, self).__init__()
        with self.init_scope():
            self.c1 = L.ConvolutionND(1, None, channels, 3, stride, dilate, dilate=dilate)
            self.c2 = L.ConvolutionND(1, None, channels, 3, 1, dilate, dilate=dilate)
            if stride > 1:
                self.cd = L.ConvolutionND(1, None, channels, 1, stride, 0)
            self.b1 = L.BatchNormalization(channels)
            self.b2 = L.BatchNormalization(channels)

    def __call__(self, x):
        h = F.relu(self.b1(self.c1(x)))
        if self.stride > 1:
            res = self.cd(x)
        else:
            res = x
        h = res + self.b2(self.c2(h))
        return F.relu(h)


class ResBlock(chainer.Chain):

    def __init__(
            self,
            channels,
            n_block,
            dilate=1
    ):
        self.n_block = n_block
        super(ResBlock, self).__init__()
        with self.init_scope():
            self.b0 = BaseBlock(channels, 2, dilate)
            for i in range(1, n_block):
                bx = BaseBlock(channels, 1, dilate)
                setattr(self, 'b{}'.format(str(i)), bx)

    def __call__(self, x):
        h = self.b0(x)
        for i in range(1, self.n_block):
            h = getattr(self, 'b{}'.format(str(i)))(h)
        return h


class ResNet34(chainer.Chain):

    def __init__(self):
        super(ResNet34, self).__init__()
        with self.init_scope():
            self.conv1 = L.ConvolutionND(1, None, 64, 7, 2, 3)
            self.bn1 = L.BatchNormalization(64)
            self.resblock0 = ResBlock(64, 3)
            self.resblock1 = ResBlock(128, 4)
            self.resblock2 = ResBlock(256, 6)
            self.resblock3 = ResBlock(512, 3)
            self.fc = L.Linear(None, 2)

    def __call__(self, x):
        h = F.relu(self.bn1(self.conv1(x)))
        h = F.max_pooling_nd(h, 3, 2)
        for i in range(4):
            h = getattr(self, 'resblock{}'.format(str(i)))(h)
        h = F.average(h, axis=2)
        h = self.fc(h)
        return h
      

class Classifier(chainer.Chain):

    def __init__(
            self,
            predictor,
            lossfun=F.softmax_cross_entropy
    ):
        super(Classifier, self).__init__()
        with self.init_scope():
            self.predictor = predictor
            self.lossfun = lossfun

    def __call__(self, *args):
        assert len(args) >= 2
        x = args[:-1]
        t = args[-1]
        y = self.predictor(*x)

        # loss
        loss = self.lossfun(y, t)
        with chainer.no_backprop_mode():
            # other metrics
            accuracy = F.accuracy(y, t)
        # reporter
        reporter.report({'loss': loss}, self)
        reporter.report({'accuracy': accuracy}, self)

        return loss

    def predict(self, x):
        with chainer.function.no_backprop_mode(), chainer.using_config('train', False):
            x = Variable(self.xp.asarray(x, dtype=self.xp.float32))
            y = self.predictor(x)
            return y
学習を実行するための準備として,以下の関数を用意します.
  • create_train_dataset():学習用データセットを ECGDataset クラスに渡す
  • create_trainer():学習に必要な設定を行い,Trainerオブジェクトを作成
python
from chainer import optimizers
from chainer.optimizer import WeightDecay
from chainer.iterators import MultiprocessIterator
from chainer import training
from chainer.training import extensions
from chainer.training import triggers
from chainer.backends.cuda import get_device_from_id


def create_train_dataset(root_path):
    train_path = os.path.join(root_path, 'preprocessed', 'train')
    train_dataset = ECGDataset(train_path)

    return train_dataset


def create_trainer(
    batchsize, train_dataset, nb_epoch=1,
    device=0, lossfun=F.softmax_cross_entropy
):
    # setup model
    model = ResNet34()
    train_model = Classifier(model, lossfun=lossfun)

    # use Adam optimizer
    optimizer = optimizers.Adam(alpha=0.001)
    optimizer.setup(train_model)
    optimizer.add_hook(WeightDecay(0.0001))

    # setup iterator
    train_iter = MultiprocessIterator(train_dataset, batchsize)

    # define updater
    updater = training.StandardUpdater(train_iter, optimizer, device=device)

    # setup trainer
    stop_trigger = (nb_epoch, 'epoch')
    trainer = training.trainer.Trainer(updater, stop_trigger)
    logging_attributes = [
        'epoch', 'iteration',
        'main/loss', 'main/accuracy'        
    ]
    trainer.extend(
        extensions.LogReport(logging_attributes, trigger=(2000 // batchsize, 'iteration'))
    )
    trainer.extend(
        extensions.PrintReport(logging_attributes)
    )
    trainer.extend(
        extensions.ExponentialShift('alpha', 0.75, optimizer=optimizer),
        trigger=(4000 // batchsize, 'iteration')
    )

    return trainer
これで学習の準備が整ったので,関数を呼び出してtrainerを作成します.
python
train_dataset = create_train_dataset(dataset_root)
python
trainer = create_trainer(256, train_dataset, nb_epoch=1, device=0)
それでは学習を開始しましょう. (1分30秒程度で学習が完了します.)
python
%time trainer.run()
epoch iteration main/loss main/accuracy 0 7 1.1307 0.868304
0 14 0.273237 0.923549
0 21 0.0950281 0.96596
0 28 0.058416 0.982701
0 35 0.0658751 0.982143
0 42 0.0506687 0.985491
0 49 0.057125 0.986607
0 56 0.063658 0.986607
0 63 0.0600915 0.981027
0 70 0.0391555 0.988281
0 77 0.0325103 0.991629
0 84 0.0344455 0.987723
0 91 0.0281526 0.989955
0 98 0.0266191 0.991629
0 105 0.0318078 0.990513
0 112 0.0304052 0.991071
0 119 0.0293185 0.993304
0 126 0.0290823 0.989397
0 133 0.019204 0.996094
0 140 0.0177221 0.994978
0 147 0.0218593 0.990513
0 154 0.019589 0.994978
0 161 0.0257332 0.991629
0 168 0.0155559 0.99442
0 175 0.0161097 0.99442
0 182 0.0212924 0.992746
CPU times: user 1min 12s, sys: 14.7 s, total: 1min 27s Wall time: 1min 27s
学習が問題なく進めば,main/accuracyが0.99 (99%)付近まで到達していると思います.

評価

学習したモデルを評価用データに当てはめて識別性能を確認するため,以下の関数を用意します.
  • create_test_dataset() : 評価用データの読み込み
  • predict() : 推論を行い,結果の配列(正解ラベルと予測ラベル)を出力
  • print_confusion_matrix() : 予測結果から混同行列とよばれる表を出力
  • print_scores() : 予測結果から予測精度の評価指標を出力
python
from chainer import cuda
from sklearn.metrics import classification_report
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix


def create_test_dataset(root_path):
    test_path = os.path.join(root_path, 'preprocessed', 'test')
    test_dataset = ECGDataset(test_path)
    return test_dataset

  
def predict(trainer, test_dataset, batchsize, device=-1):
    model = trainer.updater.get_optimizer('main').target
    ys = []
    ts = []
    for i in range(len(test_dataset) // batchsize + 1):
        if i == len(test_dataset) // batchsize:
            X, t = zip(*test_dataset[i*batchsize: len(test_dataset)])
        else:
            X, t = zip(*test_dataset[i*batchsize:(i+1)*batchsize])
        X = cuda.to_gpu(np.array(X), device)
        y = model.predict(X)
        y = cuda.to_cpu(y.data.argmax(axis=1))
        ys.append(y)
        ts.append(np.array(t))
    return np.concatenate(ts), np.concatenate(ys)

  
def print_confusion_matrix(y_true, y_pred):
    labels = sorted(list(set(y_true)))
    target_names = ['Normal', 'VEB']
    cmx = confusion_matrix(y_true, y_pred, labels=labels)
    df_cmx = pd.DataFrame(cmx, index=target_names, columns=target_names)
    plt.figure(figsize = (5,3))
    sn.heatmap(df_cmx, annot=True, annot_kws={"size": 18}, fmt="d", cmap='Blues')
    plt.show()
    

def print_scores(y_true, y_pred):
    target_names = ['Normal', 'VEB']
    print(classification_report(y_true, y_pred, target_names=target_names))
    print("accuracy: ", accuracy_score(y_true, y_pred))
評価用データセットを用意し,
python
test_dataset = create_test_dataset(dataset_root)
評価用データに対して予測を行います. (17秒程度で予測が完了します)
python
%time y_true_test, y_pred_test = predict(trainer, test_dataset, 256, 0)
CPU times: user 15.6 s, sys: 2.73 s, total: 18.3 s Wall time: 18.3 s
それでは予測結果を確認していきましょう.
まずは, 混同行列 とよばれる,予測の分類結果をまとめた表を作成します.行方向(表側)を正解ラベル,列方向(表頭)を予測ラベルとして,各項目では以下の集計値を求めています.
  • 左上 : 実際に正常拍動であるサンプルが,正常拍動と予測された数
  • 右上 : 実際に正常拍動であるサンプルが,VEBと予測された数
  • 左下 : 実際にVEBであるサンプルが,正常と予測された数
  • 右下 : 実際にVEBであるサンプルが,VEBと予測された数
python
print_confusion_matrix(y_true_test, y_pred_test)
png
続いて,予測結果から計算される予測精度の評価指標スコアを表示してみましょう.
特に,以下のスコアに注目してみてください.
  • 適合率 (Precision) : それぞれの予測診断結果 (Normal or VEB) のうち,正しく診断できていた(正解も同じ診断結果であった)割合
  • 再現率 (Recall) : それぞれの正解診断結果 (Normal or VEB) のうち,正しく予測できていた(予測も同じ診断結果であった)割合
  • F1値 (F1-score) : 適合率と再現率の調和平均
  • 正解率 (Accuracy) : 全ての診断結果 (Normal and VEB) のうち,正しく予測できていた(予測も同じ診断結果であった)割合
また,クラスごとのスコアの下に,複数の平均スコアが表示されていますが,それぞれの意味は以下の通りです.
  • マイクロ平均 (micro avg) : 各クラスを区別せずに,混同行列全体からスコアを算出.計算結果はいずれも正解率と一致
  • マクロ平均 (macro avg) : クラスごとに算出されたスコアの単純平均
  • 重み付き平均 (weighted avg) : クラスごとに算出されたスコアをサンプル数の比率で重み付けした加重平均
python
print_scores(y_true_test, y_pred_test)
precision recall f1-score support
      Normal       0.99      0.94      0.96     42149
         VEB       0.52      0.90      0.66      3200

   micro avg       0.94      0.94      0.94     45349
   macro avg       0.76      0.92      0.81     45349
weighted avg       0.96      0.94      0.94     45349

accuracy:  0.9351694634942336
サンプル数が多い正常拍動に対する予測スコアは高い値を示す一方で,サンプル数の少ないVEBに対しては,スコアが低くなる傾向があります.今回のデータセットのように,サンプルが占めるクラスの割合が極端に偏っている不均衡データでは,こうした傾向がしばしば観測されることが知られています.
次節では,このようなクラス不均衡問題への対応をはじめとして,予測モデルを改善するための試行錯誤について幾つか紹介していきます.

精度向上に向けて

本節では,前節にて構築した学習器に対して,「データセット」「目的関数」「学習モデル」「前処理」といった様々な観点で工夫を行うことで,精度改善に寄与する方法を模索していきます.
機械学習を用いて解析を行う際には,どの工夫が精度改善に有効なのか予め分からない場合が多く,試行錯誤が必要となります.ただし,手当たり次第の方法を試すことは得策では無いので,対象とするデータセットの性質に基づいて,有効となり得る手段を検討していくことが重要となります.
まずはじめに,前節でも課題として挙がっていた,クラス不均衡の問題への対処法から検討していきましょう.

クラス不均衡データへの対応

前節でも触れたように,クラス不均衡データを用いて学習器を構築する際,大多数を占めるクラスに偏った予測結果となり,少数のクラスに対して精度が低くなってしまう場合があることが一般的に知られています.一方で,(今回のデータセットを含めて)現実世界のタスクにおいては,大多数の正常サンプルの中に含まれる少数の異常サンプルを精度良く検出することが重要であるというケースは少なくありません.こうした状況において,少数クラスの検出に注目してモデルを学習するための方策が幾つか存在します.
具体的には,
  1. サンプリング
    • 不均衡データセットからサンプリングを行い,クラス比率のバランスが取れたデータセットを作成.
      • Undersampling : 大多数の正常サンプルを削減.
      • Oversampling : 少数の異常サンプルを水増し.
  2. 損失関数の重み調整
    • 正常サンプルを異常と誤分類した際のペナルティを小さく,異常サンプルを正常と誤分類した際のペナルティを大きくする.
    • 例えば,サンプル数の存在比率の逆数を重みとして利用.
  3. 目的関数(損失関数)の変更
    • 異常サンプルに対する予測スコアを向上させるような目的関数を導入.
  4. 異常検知
    • 正常サンプルのデータ分布を仮定し,そこから十分に逸脱したサンプルを異常とみなす.
などの方法があります.本節では,「1.サンプリング」,「3.目的関数の変更」の例を紹介していきます.
サンプリング
UndersamplingOversamplingを組み合わせて,データセットの不均衡を解消することを考えます.
今回は以下のステップでサンプリングを行います.
  1. Undersamplingにより,正常拍動サンプルのみ1/4に削減 (VEBサンプルは全て残す)
    • ここでは,単純なランダムサンプリングを採用します.ランダム性があるため,分類にとって重要な(VEBサンプルとの識別境界付近にある)サンプルを削除してしまう可能性があります.
    • ランダムサンプリングの問題を緩和する手法も幾つか存在しますが,今回は使用しません.
  2. Oversamplingにより,Undersampling後の正常拍動サンプルと同数になるまでVEBサンプルを水増し
    • SMOTE (Synthetic Minority Over-sampling TEchnique) という手法を採用します.
    • ランダムにデータを水増しする最も単純な方法だと,過学習を引き起こしやすくなります.SMOTEでは,VEBサンプルと,その近傍VEBサンプルとの間のデータ点をランダムに生成してデータに追加していくことで,過学習の影響を緩和しています.
サンプリングを行うために, SampledECGDataset クラスを定義します.
また,そのクラスを読み込んで学習用データセットオブジェクトを作成する create_sampled_train_datset() 関数を用意します.
python
from imblearn.datasets import make_imbalance
from imblearn.over_sampling import SMOTE


class SampledECGDataset(ECGDataset):

    def __init__(
            self,
            path
    ):
        super(SampledECGDataset, self).__init__(path)
        _, counts = np.unique(self.y, return_counts=True)
        self.X, self.y = make_imbalance(
            self.X, self.y,
            sampling_strategy={0: counts[0]//4, 1: counts[1]}
        )
        smote = SMOTE(random_state=42)
        self.X, self.y = smote.fit_sample(self.X, self.y)

        
def create_sampled_train_dataset(root_path):
    train_path = os.path.join(root_path, 'preprocessed', 'train')
    train_dataset = SampledECGDataset(train_path)

    return train_dataset
python
train_dataset = create_sampled_train_dataset(dataset_root)
それでは先程と同様に,trainerを作成して学習を実行してみましょう.(1分程度で学習が完了します.)
python
trainer = create_trainer(256, train_dataset, nb_epoch=2, device=0)
%time trainer.run()
epoch iteration main/loss main/accuracy 0 7 1.55422 0.71317
0 14 0.282043 0.902344
0 21 0.161596 0.939174
0 28 0.109996 0.960379
0 35 0.0750687 0.976004
0 42 0.0759563 0.969866
0 49 0.0531007 0.979911
0 56 0.0477525 0.984375
0 63 0.0657228 0.981585
0 70 0.0559838 0.979911
0 77 0.0410638 0.985491
0 84 0.0311444 0.989397
1 91 0.0260796 0.993304
1 98 0.0306111 0.990513
1 105 0.0267346 0.987165
1 112 0.0161048 0.995536
1 119 0.0202833 0.991629
1 126 0.0112008 0.998326
1 133 0.0164346 0.99442
1 140 0.0189331 0.992746
1 147 0.0153746 0.99442
1 154 0.0173289 0.994978
1 161 0.01027 0.996094
1 168 0.01294 0.995536
CPU times: user 55.1 s, sys: 13.4 s, total: 1min 8s Wall time: 1min 8s
学習が完了したら,評価用データで予測を行い,精度を確認してみましょう.
python
%time y_true_test, y_pred_test = predict(trainer, test_dataset, 256, 0)

Discussion

コメントにはログインが必要です。