【競馬AI開発#18】探索的データ分析（EDA）って何をやるの？

2024年10月19日 17:03

この【競馬AI開発】シリーズでは、競馬予想AIを作ることを通して、機械学習・データサイエンスの勉強になるコンテンツの発信や、筆者が行った実験の共有などを行っていきます。

■今回やること
今回は、競馬予想AI作成に使用しているデータを可視化して分析することで、精度改善のための示唆を得ていきます。
このような取り組みを、機械学習の用語で「探索的データ分析（EDA）」といい、これまでやってきたような「とりあえず特徴量に入れて精度が上がればOK」というアプローチと異なり、データに解釈を加える取り組みになります。

Pythonにはデータ可視化に便利なSeabornというライブラリがあり、例えば次のようなグラフを1行のコードで簡単に書くことができます。

馬の性別（0:牡馬, 1:牝馬, 2:セン馬）ごとの勝率が、月によってどう変わるか

この記事では、このSeabornの基本的な使い方についても入門しつつ、探索的データ分析で押さえるべきポイント・データ可視化のパターンについて学ぶことができます。
今回の内容は、実務でもよくある「データを可視化して示唆を得たい」という場面においても広く使える内容だと思っています。

■「探索的データ分析」の役割
探索的データ分析には大きく分けて二つの役割があります。

1.データの概要を確認する分析
まず一つ目は、与えられたデータに対して、ちゃんと十分なレコード数が入っているのか（欠損値の確認）や、変なデータが存在しないか（外れ値の確認）をチェックする目的の分析です。「データの検品」みたいなものです。
数値型の変数であれば分布を、カテゴリ変数であればカテゴリ数や、各カテゴリが何行存在するのか、などを確認しておきましょう。
また、機械学習のためのデータセットであれば、目的変数との相関も軽く確認しておくと良いでしょう。

2.仮説を検証する分析
二つ目は、ドメイン知識から仮説を立てて、ピンポイントでそれを検証しにいく分析です。
データ分析は無数の切り口が存在するため、闇雲に分析を進めると「大量の分析結果を出したが、何も示唆が得られない」という状況に陥りがちです。
そうならないためには「仮説」が重要で、この仮説と「それを検証するのに適切な分析」がセットになって初めて意味をなします。
幸い、競馬ドメインの場合は調べると様々な「仮説」がすぐに見つかるので、この中から今回は以下のものを検証します。

枠順の有利/不利
- 「競馬は基本内枠有利だが、ダートコースは外枠有利」は本当か？
- 「直線レースは外枠有利」は本当か？
- 「最初のコーナーまでの距離にもよるので、そもそもコースによって有利/不利が違う」は本当か？
「適切な出走間隔があり、連続しすぎていても、休ませすぎても良くない」は本当か？
「夏は牝馬が強い」は本当か？

今回使用するソースコードや、データ取得のためのソースコードは「ソースコードのダウンロード」からダウンロード・解凍すればそのまま使うことができるので、適宜手元で参照・実行しながらお読みください。

1. データの概要を把握する分析

ここからの内容は、動画の3:56〜でも解説しています。

まずは必要なデータを読み込んでいきます。
features.csvが、加工済みのデータを全て紐付けた特徴量テーブルになっているので、このデータを使用していきます。

▼v3_0_6/src/dev.ipynb

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()  # プロットの見た目を整える

df = pd.read_csv("../data/02_features/features.csv", sep="\t")

この記事から初めて来た人のために解説すると、このテーブルは以下のような、レースごとにデータを縦に繋げたものになっています。

レース結果テーブル（https://db.netkeiba.com/race/202306010101/ ）を繋げたもの

データの分布を確認する

まずは、目的変数の分布を確認していきます。
今作成している機械学習モデルでは、「1着に入ったかどうか」の0or1ラベルを目的変数に使用しているので、その正例割合を確認しておきます。

df["target"] = (df["rank"] == 1).astype(int)
df["target"].value_counts()

target
0    264003
1     20765
Name: count, dtype: int64

df["target"].mean()

0.0729190077536802

正例割合7%の不均衡データであることを念頭においておきましょう。

このように、正例と負例のラベルが偏ったデータを不均衡データといいます。
このような場合、機械学習モデルは多数派クラスに適応しやすくなったり、単純な「正解率」でモデルの精度を測定できなくなったりするので、注意が必要です。

次に、レースid, 馬id, 騎手id, 調教師idについて、年毎のユニーク数をカウントし、各項目どのくらいのメンバーがいるのか確認しておきます。

df["year"] = pd.to_datetime(df["date"]).dt.year
df.groupby("year")[["race_id", "horse_id", "jockey_id", "trainer_id"]].nunique()

このような、id列などの「カテゴリ変数のユニーク・カテゴリ数」を「カーディナリティ」ともいうのですが、この「カーディナリティ」によって、機械学習モデルに入れる時の適切な処理（エンコーディング）方法が変わってくるので、チェックしておきます。
今作成しているモデルではダミー変数化をしていないので関係ないですが、もし例えばhorse_id列に対してダミー変数化をしてしまうと、列数が爆発してしまうので、特にダミー変数化を検討している場合には、必ずチェックしておきましょう。

■棒グラフ：barplot(pd.Series)
次に、race_idごとに出現数をカウントすることで、1レースあたりの出走数の分布を確認してみます。

df["race_id"].value_counts()

race_id
201805020404    18
202004020105    18
202010011111    18
201807030307    18
202010011103    18
                ..
202008010606     5
202009050409     5
202302020205     5
202207010608     4
201907030504     4
Name: count, Length: 20729, dtype: int64

これに対してさらに.value_counts()を行うことで、出走数の分布が確認できます。

df["race_id"].value_counts().value_counts()

count
16    6116
15    2657
14    2101
13    1691
12    1686
11    1505
10    1239
18    1238
9      956
8      688
17     431
7      281
6       99
5       39
4        2
Name: count, dtype: int64

ここで、seabornのbarplot()を使うと、これを棒グラフで可視化することができます。

sns.barplot(df["race_id"].value_counts().value_counts())

値の大小が可視化され、解釈しやすくなりました。
1レース16頭が出走するパターンが一番多いですが、少ないと4頭しか出走しないレースもあるようです。
4頭と16頭では、シンプルに出走数が多い16頭の方が勝ちにくい（＝確率値としては小さい値になる）ことが予想されますが、この「出走数」はまだ特徴量に入れていないので、入れて検証してみたいですね。

機械学習モデル開発におけるEDAの際には、このように、分析結果から「何か効きそうな特徴量を作れないか」を常に想像しながら進めていくと良いでしょう。

■カテゴリ列の出現数を可視化：countplot()
上のプロットは、sns.countplot()を使っても描くことができます。
その場合は、以下のようなDataFrameを作成しておきます。

agg_df = df["race_id"].value_counts().to_frame().reset_index()

これに対して、以下のように実行できます。

sns.countplot(data=agg_df, x="count")

このようにcountplot()は、x=で指定した列に対して「どの値が何個存在するのか」カウントして棒グラフで出力します。

この後、countplot()の他にもいくつかの可視化関数が出てきますが、seabornの関数は基本的に全て以下の形で使うことができます。

【seabornの基本的な型】
sns.hogehogeplot(data=df, x="{列名}", y="{列名}"})

次に、horse_id列についても同様に同一idの出現数を可視化してみます。

agg_df = df["horse_id"].value_counts().to_frame().reset_index()
agg_df.rename(columns={"count": "id_count"}, inplace=True)
sns.countplot(data=agg_df, x="id_count")

今回は、x軸の目盛りが重なってしまいました。
このような場合は、以下のようにplt.figure()を使うと出力サイズを調整することができます。

▼common/src/dev.ipynb

plt.figure(figsize=(20, 5))
sns.countplot(data=agg_df, x="id_count")

また、x軸が長くなってしまう場合は、以下のようにy=に集計対象列を指定することで、横棒グラフにすることもできます。

▼common/src/dev.ipynb

plt.figure(figsize=(5, 10))
sns.countplot(data=agg_df, y="id_count")

【補足】
ここで、「seabornというライブラリを使って可視化しているのに、なぜ別ライブラリであるmatplotlibのplt.figure()でサイズの調整を行うのか？」という疑問を持った方もいるかもしれません。
実は、seabornは内部でmatplotlibを使用しており、最終的な描画もmatplotlibで行っています。
そのため、seabornで作成したプロットは、matplotlibの各種関数を使って様々な調整をすることができますし、matplotlibと併用してプロットを書くこともできます。

■箱ひげ図：boxplot()
次に、jockey_idについても見てみます。

sns.countplot(data=agg_df, x="id_count")

横に広がってかなり極端な分布になっているため、先ほどのようなサイズ調整も難しそうです。

中身を見ると、320人の騎手のうち、人によっては5274回の出現がありますが、ほとんどの人は1〜2回程度の出現であることが分かります。

このような場合に特徴を把握する可視化としては、以下のような箱ひげ図というプロットもよく出てきます。

sns.boxplot(data=agg_df, x="id_count")

これはそれぞれ、

箱の両端：25パーセンタイル値と75パーセンタイル値
箱の中央の線：中央値
線（ひげ）の左端：25パーセンタイル値 - 1.5 x 箱の長さ
線（ひげ）の右端：75パーセンタイル値 + 1.5 x 箱の長さ
ひげからはみ出した点：外れ値

を表します。

https://mathlandscape.com/box-plot/ より引用

実際の統計量は.describe()で表示できるので、この値と図を比較してみると理解が深まると思います。

agg_df["id_count"].describe()

count     320.000000
mean      889.900000
std      1284.215532
min         1.000000
25%         5.000000
50%       213.500000
75%      1327.500000
max      5274.000000
Name: id_count, dtype: float64

■連続値に対する度数分布：histplot()
ここまでは「idの出現回数」などの「離散値（＝整数などのとびとびの値）」の度数（＝出現回数）を可視化する方法でした。
連続値の場合は、一定の間隔（＝ビン）にデータを区切って度数をカウントする「ヒストグラム」という可視化をします。
例えば、単勝オッズの分布をヒストグラムで確認してみます。

sns.histplot(data=df, x="tansho_odds")

一見分からないかもしれないですが、細かいビンに区切られて、ビンごとの出現回数がプロットされています。
bins=パラメータでビンの個数を指定できるので、試しに10を指定してみると、この様子がわかりやすいです。

sns.histplot(data=df, x="tansho_odds", bins=10)

粗めのヒストグラムになりました。

【発展】
sns.histplot()におけるデフォルトのビンの数は、Freedman-Diaconis ruleという手法で決められているようです。

また、離散値の列でなくても、数値型（量的変数）の列に対してはヒストグラムで度数を表現します。
例えば、馬体重はデータとしては離散値である整数で入っていますが、

df["weight"]

0         462
1         478
2         442
3         464
4         456
         ... 
284763    512
284764    498
284765    508
284766    472
284767    490
Name: weight, Length: 284768, dtype: int64

これに対してcountplot()を使うと、以下のようにx軸の表示が潰れてしまいます。

sns.countplot(df, x="weight")

countplot()はカテゴリ変数をカウントするプロットのため、x軸は全てのカテゴリを表示しようとします（「間の値」という概念がないことを意味します）。
馬体重のような列は、カテゴリ変数ではなく、数や量で計ることのできる量的変数であるため、ヒストグラムでプロットするのが適切です。

sns.histplot(df, x="weight")

ただ、入っているデータとしては離散的であるため、ビンの切り方によってガタガタなグラフになってしまいます。
そこで、このようなデータに対してはdiscrete=Trueに設定することで、データ点に存在する値でビンを設定してプロットしてくれます。

sns.histplot(df, x="weight", discrete=True)

目的変数との相関を確認する

■カテゴリ変数と目的変数の相関：barplot()
次に、いくつかの列を使って目的変数との相関を確認してみます。
相関を見る方法は様々ありますが、説明変数がカテゴリ変数の場合は、barplot()を使う方法がおすすめです。

例えば、以下のように書くと「枠番と目的変数（＝勝率）の相関」を簡単に可視化することができます。

▼common/src/dev.ipynb

sns.barplot(df, x="wakuban", y="target")

内部で「x（=wakuban）ごとのy（=target）の平均」が集計され、さらに誤差棒（＝この平均値がどのくらい信頼できるものなのか）も一発で描いてくれるという、優れモノです。

この誤差棒はデフォルトだと95%信頼区間を表しており、errorbar=パラメータで調整することができます。（公式document）

グラフを見てみると、「内枠有利」の傾向は、少なくともこのデータセットからは言えないということが分かります。（むしろ、外枠の方がやや有利に見える）
枠順については、また後ろの方で詳細分析します。

■量的変数と目的変数の相関：violinplot()
次に、説明変数が量的変数の場合には、バイオリンプロットを使って相関を確認する方法も有力です。
例えば、「馬体重と勝率の関係」を見たい場合、barplot()ではx軸にカテゴリ変数を指定しなければならないため、以下のようなプロットをすることになります。

sns.barplot(df, x="target", y="weight")

これでも平均値は比べられるのですが、バイオリンプロットを描くと分布も併せて確認することができます。

sns.violinplot(df, x="target", y="weight")

わずかですが、勝った馬の方が体重が重い方に分布していることが分かります。

■二重軸グラフ：twinx()
また、次のように直接分布を重ねてプロットしてしまう方法もあります。

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()  # 2重軸を作成

# 一つ目の軸（ax）に対するプロット
sns.histplot(
    df.query("target == 0"),
    x="weight",
    discrete=True,
    ax=ax1,
    label="target:0"
)
ax1.legend(loc="upper left")
# 二つ目の軸（ax）に対するプロット
sns.histplot(
    df.query("target == 1"),
    x="weight",
    discrete=True,
    ax=ax2,
    color="red",
    label="target:1"
)
ax2.legend()

ax1が1つ目の軸（=左側の軸）を表しており、ax2が2つ目の軸（=右側の軸）を表しているので、例えばax1.legend()は「一つ目の軸に対するデータの凡例（label="target:0"）をプロットする」コードを表しています。

他にも、このax1, ax2を使ってグラフのスタイル設定をすることができます。
例えば、軸のラベルは.set_ylabel()で設定できるので、以下のようにすると軸名を区別して分かりやすくできます。

# 2重軸
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

sns.histplot(
    df.query("target == 0"),
    x="weight",
    discrete=True,
    ax=ax1,
    label="target:0"
)
ax1.legend(loc="upper left")
ax1.set_ylabel("count(target:0)")

sns.histplot(
    df.query("target == 1"),
    x="weight",
    discrete=True,
    ax=ax2,
    color="red",
    label="target:1"
)
ax2.legend()
ax2.set_ylabel("count(target:1)")

ちなみに、各オブジェクト名・設定メソッドは、以下のような対応関係になっているので、何となく頭に入れておくと便利です。

https://www.skytowner.com/explore/getting_started_with_matplotlib より引用
今回であれば、一つのfigureの中に二つのax（axes）を2重軸（twinx)の形式でプロットしている

https://matplotlib.org/stable/users/explain/quick_start.html より引用

■連続値 x 連続値の可視化: scatterplot(), displot()
今回の目的変数は0or1の離散値なので登場しませんでしたが、一応、連続値 x 連続値で分布を可視化をしたい場合も見ていきます。
例えば、特徴量の一つに「rank_mean_10races_relative：直近10レースの着順平均をレース内で標準化したもの」があります。これと「単勝オッズ」の相関を確認してみましょう。
まず一番単純なのは、各点を全てプロットする散布図です。

sns.scatterplot(data=df, x="rank_mean_10races_relative", y="tansho_odds")

ただし、見て分かるように、これだと「どこにプロットが集中しているのか」を把握することが難しくなってしまいます。

そこで、displot()を使うと、密度も合わせてプロットすることができます。

sns.displot(data=df, x="rank_mean_10races_relative", y="tansho_odds")

追加で使用しているテーブルの紐付き率

また、特徴量作成などの際に別テーブルを使用している場合は「そのテーブルのデータがどのくらい紐づいているのか」も確認しておきます。
例えば、今のモデルでは馬の過去成績を集計した特徴量を作成していますが、そもそも「過去成績」が存在するレコードはどの程度の割合を占めるのでしょうか？
これは、"rank_3races"列が欠損値かどうかで判定できるので、集計してみます。

# 「馬の過去成績」の欠損率
df["rank_3races"].isnull().mean()

0.09870139903359929

約10%のレコードは、馬の過去成績が存在していないことが分かります。
今の特徴量重要度上位を見ると馬の過去成績特徴量がきていますが、この10%程度のデータはその特徴量を使えていないことになり、うまく学習ができていない可能性があるので、要確認です。
場合によっては、母集団から除外するような対処も考えられます。

次に、馬の過去成績が紐づいている場合、何レース分の成績が紐づいているのか、分布を確認します。
FeatrueCreator.create_baselog()を実行すると、「母集団に過去成績を紐付けたテーブル」が作成されるので、このテーブルで確認していきます。

fc = FeatureCreator()
fc.create_baselog()
fc.baselog

このFeatureCreatorクラスの中身・作り方については、以下の記事で詳しく解説しているので、適宜参照してください。

このテーブルは、母集団の（race_id, horse_id）に対して、そのレースより過去の「馬の過去成績」の行を全て紐付けたものになっているので、行数をカウントすると、「過去何走の成績が紐づいているか」を集計することができます。

n_races = fc.baselog[["race_id", "horse_id"]].value_counts().to_frame()

最大で、122レースの過去成績が紐づいている馬が存在するようです。
分布を確認してみましょう。

sns.histplot(n_races, x="count", discrete=True)

n_races["count"].describe()

count    258813.000000
mean         11.025057
std          10.315446
min           1.000000
25%           3.000000
50%           8.000000
75%          16.000000
max         122.000000
Name: count, dtype: float64

平均だと11レース、中央値だと8レースの過去成績が紐づいているようです。
この値は時系列データにおいて特徴量を作る際に、「過去何行のデータを集計するのが良いか」の目安の値になります。

この「過去何行のデータを集計するか」をwindow sizeともいいます。

2. 仮説を検証する分析

ここからの内容は、動画の37:28〜でも解説しています。

最初に書いたように、「仮説を立てて、それが本当なのかをデータから見極める」というアプローチも、データ分析では重要なものになってきます。

今回検証する仮説を再掲します。

枠順の有利/不利
- 「競馬は基本内枠有利だが、ダートコースは外枠有利」は本当か？
- 「直線レースは外枠有利」は本当か？
- 「最初のコーナーまでの距離にもよるので、そもそもコースによって有利/不利が違う」は本当か？
「適切な出走間隔があり、連続しすぎていても、休ませすぎても良くない」は本当か？
「夏は牝馬が強い」は本当か？

芝は内枠有利・ダートは外枠有利？

先ほどの分析で、データ全体で見ると、内枠・外枠の差はみられないことが分かりましたが、枠順に関しては「ダートの内枠はコースが荒れるので不利/芝はその影響を受けないので内枠有利」ということがよく言われています。

これが本当なのか確かめてみましょう。
今回も「勝率」を集計したいのでbarplot()を使いますが、その際にhue=を指定すると、以下のようにグループごとに分割してプロットすることができます。

sns.barplot(df, x="race_type", y="target", hue="wakuban")

見ると、確かにダートのレースに関しては外枠の方が有利な傾向がありそうです。
芝のレースはわずかに内枠有利な傾向がありそうですが、誤差の範囲にあり、そこまで有意な差にはなっていないように見えます。

このようにbarplot()では「その差分がどの程度信用できるのか」も含めて一目で把握できる点が、クロス集計にはない強みですね。

直線レースは外枠有利？

次に、直線レースの内枠は、外枠と距離が変わらないにも関わらず、コースが荒れているデメリットのみ食らうので、不利なことが仮説として考えられます。
around列で直線レースを判別できるので、これも可視化してみましょう。

sns.barplot(df, x="around", y="target", hue="wakuban")

これはかなり差が開いていますね。直線レースでは明らかに外枠有利な傾向が見られます。

出走間隔は何週空けるのがベストか？

「馬には適切な出走間隔がある」という話がよく言われていますが、本当でしょうか？これも検証してみましょう。
FeatureCreator.baselogのdate_horse列の最新日付からの経過日数を見れば良いので、まずは各レースの「直近過去成績の日付」を集計します。

latest_df = (
    fc.baselog
    .groupby(["race_id", "horse_id", "date"])["date_horse"]
    .max()
    .to_frame()
    .reset_index()
)

このdate（出走日）とdate_horse（直近過去成績日付）の差が出走間隔になるので、まずはこの分布を確認してみます。

latest_df["interval"] = (
    pd.to_datetime(latest_df["date"]) - pd.to_datetime(latest_df["date_horse"])
).dt.days
latest_df["interval"] = latest_df["interval"] // 7  # 週単位にする
sns.histplot(latest_df, x="interval", discrete=True)

次に、このintervalごとの勝率をプロットしてみます。
intervalが30週を超えるとかなりデータ数が少なそうなので、30未満に絞って可視化します。

df2 = latest_df.merge(
    df[["race_id", "horse_id", "target"]],
    on=["race_id", "horse_id"]
)
sns.barplot(df2.query("interval < 30"), x="interval", y="target")

見ると、確かに出走間隔が小さすぎても大きすぎても勝率が低くなってしまうことが分かります。
これは特徴量としても素直にそのまま入れるだけで効きそうですね。

夏は牝馬が強い？

「牝馬は春に生殖期を終えて、夏はホルモンバランスが安定するので、パフォーマンスが安定して勝ちやすい」ということが言われていますが、これも検証してみましょう。

これまで通りbarplot()を使っても良いですが、lineplot()を使うと時系列での変化も追うことができます。

df["month"] = pd.to_datetime(df["date"]).dt.month
sns.lineplot(df, x="month", y="target", hue="sex", marker=".")

確かに、牝馬は夏になると明確に勝率が上がり、8月には牡馬の勝率を上回っていることが分かります。
一方、少なくともこのデータセットでは、よく言われる「夏は牝馬」という言葉ほど「夏中牝馬が牡馬の勝率を上回っている」というわけではなく、「夏になると勝率がイーブンになる」くらいの変化なようです。（回収率が高いことを指した言葉なんですかね）

ここで、この集計は牡牝混合戦でないと意味がないので、混合戦に絞った結果も一応確認しておきます。

mix_races = df.groupby("race_id")[["sex"]].nunique().query("sex == 2").index
sns.lineplot(
    df.query("race_id in @mix_races"),
    x="month",
    y="target",
    hue="sex",
    marker="."
)

傾向としてはほぼ変わらないようです。

ソースコードのダウンロード

今回の部分に関係あるもののみ抜粋して載せているので、

#17までのソースコードを持っている方
- 元のソースコードのv3_0_6/src/にv3_eda/src/eda.ipynbのみコピーしてください
今回初めての方
- common/src/main.ipynbを実行し、必要なデータを収集してください
- v3_eda/src/main.ipynbを実行すると、v3_eda/data/02_features/features.csvが作成されます

また、定期購読をすると、今月更新される全ての記事に加えて#1〜#4の記事が980円で全て読めるので、単品購入より圧倒的にお得です。