🐴

競馬データを用いたお試し解析(R)!

2024/12/18に公開

analytics

date

tech

はじめに

みなさま、データ解析をどのように行っていますか？
最近、Pythonを使って解析をしている方をよく見かけますが、Rもとても使いやすいツールです。特にRStudioを利用することで、さらに快適にデータ解析を進めることができます。

今回は、私が初めて競馬場で馬券を買って敗北⭐︎した菊花賞の競馬データ（https://race.netkeiba.com/top/race_list.html）
を題材に、Rを使った統計的な解析の一例をご紹介します。

この方法は競馬データに限らず、さまざまな実験データや研究データの解析にも応用できる部分があるので、ぜひ参考にしてみてください。(実験や研究データを用いて解析する場合はデータセットの作成以降が参考になると思います)

この記事を通じて、RとRStudioの使いやすさを感じていただければ幸いです！

今回作るもの

本記事では、以下の内容を進めていきます：

菊花賞の競馬データを用いたデータ解析
データスクレイピングと加工
統計的な解析

環境

R version: 4.3.1 (2023-06-16)
RStudio version: 2023.6.0.421
OS: macOS

全体の流れ

環境設定
データのスクレイピング
データセットの作成
統計的な解析

環境設定

まずはRとRStudioのインストールを済ませましょう。インストールが完了したら、RStudioを開き、以下の手順を進めてください。

RStudio左上の「File」から「New File」→「R Notebook」を選択します。
Option + Command + I を押して新しいRチャンクを挿入します。
(これ以降はコードを書くたびにRチャンクを挿入、実行の繰り返し)
今回の解析に必要なパッケージをインストールします。

# 必要なパッケージをインストール
install.packages("rvest")
install.packages("httr")
install.packages("dplyr")
install.packages("readr")
install.packages("RSelenium")
install.packages("data.table")
install.packages("writexl")
install.packages("ggplot2")
install.packages("ARTool")

パッケージを読み込みます。

library(rvest)# Webスクレイピング用
library(httr) # 文字化け対策
library(dplyr)# データ操作用
library(readr)
library(RSelenium)
library(data.table)
library(writexl)
library(ggplot2)# 図のプロット
library(ARTool)

データのスクレイピング

では早速データの取得から始めましょう。今回は2024年の菊花賞データを以下のコードで取得します。

url <- "https://db.netkeiba.com/race/202408050611/"  # 2024 菊花賞のURLを指定
html <- read_html(url, encoding = "euc-jp") # ページを取得しHTMLを解析

Sys.sleep(1) # アクセス負荷を軽減するための待機

tbl <- html_element(html, css = "#contents_liquid > table") # 指定したCSSセレクタから表を抽出
df <- html_table(tbl) # HTML表をデータフレームに変換

# データの保存
fwrite(df, "kikka_2024_results.csv", row.names = FALSE, bom = TRUE) # CSV形式で保存
write_xlsx(df, "kikka_2024_results.xlsx") # Excel形式で保存

ここでは、fwrite 関数を用いてCSVファイルに保存しています。この関数は data.table パッケージの一部で、大規模データの書き出しを高速に行うことができます。また、bom = TRUE を指定することで、ファイルをExcelで開いた際に文字化けを防止します。

さらに、2019年から2024年までのデータをまとめて取得してみましょう。

# 菊花賞 (2019-2024) のデータ取得
years <- 2019:2024 # 対象年
base_url <- "https://db.netkeiba.com/race/" # 基本URL
race_ids <- c("201908040711", "202008040611", "202109040611", "202209040711", "202308020711", "202408050611") # レースID
urls <- paste0(base_url, race_ids, "/") # 各年のURLを生成

all_data <- list() # データ格納用リスト

for (i in seq_along(urls)) {
  url <- urls[i]
  year <- years[i] # 対応する年
  
  html <- read_html(url, encoding = "euc-jp") # HTMLを読み込み
  Sys.sleep(1) # アクセス負荷軽減
  
  tbl <- html_element(html, css = "#contents_liquid > table") # 表データの抽出
  df <- html_table(tbl, fill = TRUE) # データフレームに変換
  
  df$Year <- year # 年の列を追加
  
  all_data[[i]] <- df # リストに追加
}

# データの結合
final_data <- bind_rows(all_data) # 全データを結合

# CSVファイルに保存
fwrite(final_data, "kikka_2019_2024_results.csv", row.names = FALSE, bom = TRUE)

# Excelファイルに保存
write_xlsx(final_data, "kikka_2019_2024_results.xlsx")

ここでは、html_table(tbl, fill = TRUE) を用いることで欠損値を含む場合でもデータフレームとして正しく読み込むことができます。また、ループ内で Sys.sleep(1) を使用してアクセス間隔を調整し、サーバーへの負荷を軽減しています。

データセットの作成

次に、先ほど作成したデータファイルを読み取ってデータセットを作成してみましょう。

path <- "./kikka_2019_2024_results.csv"
data <- read.csv(path, na = "")　# CSVファイルを読み込み
head(data)

ここでは、na = "" を指定することで、空欄を欠損値 (NA) として読み込みます。これにより、解析の際に欠損値を適切に処理できるようになります。

実行結果としてhead(data)によりデータセットの最初６つのデータが表示されます

統計的な解析

今回は、馬体重の増減とタイムの関係(なんとなく)や、特定の騎手（ルメール騎手）が本当に勝率が高いのか（ルメール騎手に賭ければ大丈夫とよく聞くから笑）を調べていきます。
データ解析を行う上で個人的に大事だと思うのは、まずデータの分布を確認することです。
そのため、解析の手順として以下のように進めていきます：

データを確認 (分布を見てどのような傾向があるのかな？とか正規性ありそうだなとかどのようなデータがあるのかなとか)
データの正規性を確認 (正規性があるかないかで使う統計モデルが変わってくるため)
適したモデルを用いて検定

データを確認

まずは、どのようなデータがあるのかを確認しましょう。

colnames(data) #コンソールで確認すると便利

このコードを実行すると、データフレームの列名が取得できます。出力は以下の通りです：

 [1] "着順"       "枠番"       "馬番"       "馬名"       "性齢"       "斤量"       "騎手"      
 [8] "タイム"     "着差"       "ﾀｲﾑ指数"    "通過"       "上り"       "単勝"       "人気"      
[15] "馬体重"     "調教ﾀｲﾑ"    "厩舎ｺﾒﾝﾄ"   "備考"       "調教師"     "馬主"       "賞金.万円."
[22] "Year"

次に、特定の列にどのようなデータが含まれているかを確認します。

table(data$性齢)

出力例：

牝3 牡3 
  2 105

この結果から、菊花賞の過去6年間で牝馬は2頭しか出場していないことがわかります。（牡馬が多い！）

次に、タイム列のデータ型を確認します。

str(data$タイム)

その結果、以下のようにchr型だと表示されます。

chr [1:107] "3:06.0" "3:06.0" "3:06.2" "3:06.3" "3:06.4" "3:06.6" "3:06.6" "3:06.6" "3:06.7" ...

タイム列が文字型（chr型）であるため、数値型に変換する必要があります。
以下のコードで新しいTime列を作成します。

data$Time <- sapply(data$タイム, function(x) {
  # 分と秒に分割
  time_parts <- strsplit(as.character(x), ":")[[1]]
  minutes <- as.numeric(time_parts[1])  # 分を抽出して数値型に変換
  seconds <- as.numeric(time_parts[2])  # 秒を抽出して数値型に変換
  # 合計秒数を計算
  total_seconds <- minutes * 60 + seconds
  return(total_seconds)
})

続いて、タイムの分布をヒストグラムで可視化します。

ggplot(data, aes(x = Time)) +
  geom_histogram(binwidth = 2, fill = "blue", color = "black", alpha = 0.7) +
  labs(
    title = "Distribution of Horse Race Times (Seconds)",
    x = "Time (seconds)",
    y = "Frequency"
  ) +
  theme_minimal()

この結果から、基本的に185秒より前のタイムを持つ馬が速い馬と言えそうですね！

実際に2024の菊花賞で1位のアーバンシックは184.1秒でした。

data[data$Year == 2024 & data$着順 == 1, c("馬名", "Time")]

             馬名  Time
90 アーバンシック 184.1

馬体重の増減とタイムの関係

まず、現在の馬体重列は"340(+2)"のような形式になっているため、増減部分と基準馬体重に分解し、増減がプラス・マイナス・変化なしのカテゴリに分類します。

# 増減部分を抽出してカテゴリ化
data$馬体重変化 <- gsub(".*\\(([-+0-9]+)\\)", "\\1", data$馬体重)  # 数値部分を抽出
data$馬体重変化 <- as.numeric(gsub("\\+|\\(|\\)", "", data$馬体重変化))  # "+"や括弧を除去して数値化

# カテゴリ列を作成
data$馬体重カテゴリ <- ifelse(data$馬体重変化 > 0, "プラス", 
                            ifelse(data$馬体重変化 < 0, "マイナス", "変化なし"))

次に、タイムの正規性をカテゴリごとに確認します。

by(data$Time, data$馬体重カテゴリ, function(x) shapiro.test(x))

by()関数は、データをグループ化して適用するのに便利な関数です。
この場合、馬体重カテゴリごとにShapiro-Wilk検定を実行しています。

出力例：

data$馬体重カテゴリ: プラス

	Shapiro-Wilk normality test

data:  x
W = 0.94438, p-value = 0.03123

------------------------------------------------------------------------ 
data$馬体重カテゴリ: マイナス

	Shapiro-Wilk normality test

data:  x
W = 0.68764, p-value = 3.635e-08

------------------------------------------------------------------------ 
data$馬体重カテゴリ: 変化なし

	Shapiro-Wilk normality test

data:  x
W = 0.91634, p-value = 0.08425

プラス・マイナスカテゴリはp値が0.05未満であり、正規分布ではないと判断されます。
変化なしカテゴリはp値が0.05以上で正規分布に従う可能性があります。

すべてのカテゴリで正規性が満たされないため、今回はノンパラメトリック手法のART ANOVAを使用して検定を行います。

# 馬体重カテゴリを因子型に変換
data$馬体重カテゴリ <- as.factor(data$馬体重カテゴリ)

そしてART ANOVAを実行します。

art_model <- art(Time ~ 馬体重カテゴリ, data = data)
anova(art_model)

出力例：

Analysis of Variance of Aligned Rank Transformed Data

Table Type: Anova Table (Type III tests) 
Model: No Repeated Measures (lm)
Response: art(Time)

                 Df Df.res F value  Pr(>F)  
1 馬体重カテゴリ  2    104  0.4489 0.63956  
---
Signif. codes:   0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p値が0.05を超えているため、馬体重の増減がタイムに与える影響は統計的に有意ではないと判断されます。

騎手とトップ3の関係

最後に、騎手の1人であるルメール騎手がトップ3に入りやすいかを検証します。
流れはこのようになります。

トップ3に入ったかどうかを示す新しい列を作成
ルメール騎手とその他の騎手を比較する列を作成
トップ3に入りやすいかを統計的に検証
ルメール騎手がトップ3に入った確率を計算

まず、トップ3に入ったかどうかを示す新しい列を作成します。
そこで着順が1～3位であれば「1」、それ以外は「0」とするバイナリ変数を作成します。

data$トップ3 <- ifelse(data$着順 <= 3, 1, 0)

先ほどと同じように騎手がルメール騎手であれば「1」、それ以外は「0」とするバイナリ変数を作成します。

data$ルメール <- ifelse(data$騎手 == "ルメール", 1, 0)

そして「ルメール騎手」と「その他の騎手」で、トップ3に入る確率に有意差があるかどうかを検証するためにカイ二乗検定を行います。

# クロス集計表を作成
table_top3 <- table(data$トップ3, data$ルメール)

# カイ二乗検定を実行
chi_test <- chisq.test(table_top3)

# 結果を表示
print(chi_test)

出力例：

Pearson's Chi-squared test with Yates' continuity correction

data:  table_top3
X-squared = 10.6, df = 1, p-value = 0.001131

p値が0.05未満であるため、ルメール騎手がトップ3に入りやすいことが統計的に示されました。

さらに、ルメール騎手がトップ3に入った確率を計算します：

prop_ルメール <- sum(data$トップ3[data$ルメール == 1]) / sum(data$ルメール)

# 結果を表示
cat("ルメール騎手のトップ3率:", prop_ルメール, "\n")

出力例：

ルメール騎手のトップ3率: 0.8

ルメール騎手は直近で出場した試合のうち80%の確率でトップ3に入ることが示されました！
(これぜ〜たいルメール騎手に賭けた方が良さそうじゃん！笑)

おわりに

最後までお読みいただきありがとうございました。
競馬データは分析の題材として非常にユニークであり、楽しみながら統計解析を学ぶことができます。
この記事が少しでも皆さんの解析の参考になれば幸いです。

ぜひ、自分のデータでもRを使った解析に挑戦してみてください！

Rの補足

(個人的に)よく使うショートカットキー

Option + Command + I: 新しいRチャンクを挿入
Option + Command + C: コメントアウト

(個人的に)よく使うコード

colnames(data) :　データフレームやマトリックスの列名を取得
table(data$列名) :　特定の列のデータにおける値の頻度をカウント
data_new <- data[c("列名1", "列名2"), ]:　特定の列だけを抽出して、新しいデータセットを作成

はじめに

今回作るもの

環境

全体の流れ

環境設定

データのスクレイピング

データセットの作成

統計的な解析

データを確認

馬体重の増減とタイムの関係

騎手とトップ3の関係

おわりに

Rの補足

(個人的に)よく使うショートカットキー

(個人的に)よく使うコード

Discussion