バルセロナのゴール期待値をLightGBMで予測する

前回のデータ可視化編に続き、モデル構築編になります。

僕がゆったりと記事を書いている間に新たに3シーズン分のデータが追加されてしまったので、改めて昨18/19シーズンのシュートのxGを予測していきます。

前々回、前回の記事を読んでおくと内容をスムーズに理解できると思うので、もしよろしければ以下からどうぞ。

今回解く問題をもう一度確認しておきます。

StatsBomb Open Dataを使用し、FCバルセロナの18/19シーズンにおける全シュートのゴール期待値（xG）を予測する。
※ゴール期待値（Expected Goals, xG）…シュートが得点につながる確率を0～1の範囲で表したもの
training set：04/05~17/18シーズンのバルセロナのシュートデータ
test set：18/19シーズンのバルセロナのシュートデータ
使用モデル：LightGBM

LightGBMとは

今回使用するLightGBMについて簡単に紹介しておきます。

LightGBMとは、決定木アルゴリズムをベースとした勾配ブースティング（GBDT）の機械学習フレームワークです。

速さと精度の高さが売りであり、データ分析コンペのKaggleでも高い使用率を誇っています。

“xG”をどう解くか

上述のとおり、ゴール期待値（xG）とはシュートが得点につながる確率を示すものです。

回帰の手法を使ってxGを算出するか、それともゴール or ノーゴールを予測する二値分類として捉えクラス事後確率をxGとするかは議論が別れるところです。

今回は後者の方法を採用しますが、その場合、事後確率が実際の発生確率とズレるモデルがあることに注意しなければなりません（参考）。

個人的には、クラス事後確率がcalibrateされるロジスティック回帰を使うのが良いと考えていましたが、GBDTでもlog lossを最適化することで事後確率がcalibrateされそうなので、今回はLightGBMを使うことにしました。

つまり、LightGBMを二値分類器として使い、ゴールラベルに対する事後確率をxGとします。

確率のcalibrationについては後述します。

モデルの学習

最初に、学習に使用する特徴量を指定します。

カテゴリカル変数に対してはLabel encodingを行うため、カテゴリカル変数と連続変数を別々に指定しておきます。

# select features
cat_features =  ['play_pattern', 'under_pressure', 'body_part', 'type', 'technique', 'first_time', 'one_on_one', 'open_goal'] # カテゴリカル変数
num_features = ['distance', 'angle'] # 連続変数

# select features

cat_features = ['play_pattern', 'under_pressure', 'body_part', 'type', 'technique', 'first_time', 'one_on_one', 'open_goal'] # カテゴリカル変数

num_features = ['distance', 'angle'] # 連続変数

試行錯誤した結果、以下の特徴量を学習に使うことにしました。特徴量の解説を以下に記しておきます。

play_pattern：シュートに繋がったプレーのパターン
under_pressure：shot-takerが相手選手にプレッシャーを受けていたかどうか
body_part：シュートに使った部位
type：シュートの種類
technique：シュートに使われたテクニック（ボレー、ヒールなど）
first_time：ダイレクトシュートか否か
one_on_one：GKとの1対1か否か
open_goal：シュート時にゴールが空いていたか否か
distance…シュート位置からゴール中心部までの距離
angle…シュート位置→右ポストのベクトルと、シュート位置→左ポストのベクトルが成す角度

次に、シュートデータを特徴量と教師ラベルに分割します。

教師ラベルには、当該シュートがゴールだった否かを表すgoal列を指定します。

# split into features and label
X = shots.loc[:, cat_features + num_features]
y = shots['goal']

# split into features and label

X = shots.loc[:, cat_features + num_features]

y = shots['goal']

カテゴリカル変数には、category_encodersというライブラリを使用しLabel encodingを施しておきます。

from category_encoders import OrdinalEncoder

# encodes categorical features as ordered feature
oe = OrdinalEncoder(cols=cat_features)
X = oe.fit_transform(X)

from category_encoders import OrdinalEncoder

# encodes categorical features as ordered feature

oe = OrdinalEncoder(cols=cat_features)

X = oe.fit_transform(X)

実際に学習器に渡すデータは以下のようになります。10個の特徴量をLightGBMに食わせ、特徴量とシュート結果の関係性を学習してもらいます。

データをtraining setとtest setに分け、学習を行っていきます。

# train-test splitting
num_test = len(shots[shots['season']=='2018/2019']) # number of 18/19 season data
# train: 04/05 - 16/17 season
X_train = X.iloc[:-num_test]
y_train = y.iloc[:-num_test]
# test: 18/19 season
X_test = X.iloc[-num_test:]
y_test = y.iloc[-num_test:]

# train-test splitting

num_test = len(shots[shots['season']=='2018/2019']) # number of 18/19 season data

# train: 04/05 - 16/17 season

X_train = X.iloc[:-num_test]

y_train = y.iloc[:-num_test]

# test: 18/19 season

X_test = X.iloc[-num_test:]

y_test = y.iloc[-num_test:]

グリッドサーチ→交差検証の流れで学習させていきます。

まずはグリッドサーチ。scikit-learnのGridSearchCVを使用し、最適なハイパーパラメータの組み合わせを見つけます。

import lightgbm as lgbm
from sklearn.model_selection import GridSearchCV

# define LightGBM model
lightgbm = lgbm.LGBMClassifier(objective='binary', random_state=0)

# candidates of parameter values
params = {
         'num_leaves': [10, 31, 50, 70],
         'max_depth': [3, 5, 7, 9],
         'min_child_weight': [1, 2, 3, 4, 5],
         'reg_alpha': [0, 1, 10, 100],
         'reg_lambda': [0, 1, 10, 100], 
         }

# grid searching
gscv = GridSearchCV(estimator=lightgbm, 
                   param_grid=params, 
                   scoring='neg_log_loss',
                   cv=5, 
                   n_jobs=-1,
                   verbose=2,
                   )
gscv.fit(X_train, y_train)

# print results
print('best parameters', gscv.best_params_)
print('best log loss:', gscv.best_score_)

import lightgbm as lgbm

from sklearn.model_selection import GridSearchCV

# define LightGBM model

lightgbm = lgbm.LGBMClassifier(objective='binary', random_state=0)

# candidates of parameter values

params = {

'num_leaves': [10, 31, 50, 70],

'max_depth': [3, 5, 7, 9],

'min_child_weight': [1, 2, 3, 4, 5],

'reg_alpha': [0, 1, 10, 100],

'reg_lambda': [0, 1, 10, 100],

}

# grid searching

gscv = GridSearchCV(estimator=lightgbm,

param_grid=params,

scoring='neg_log_loss',

cv=5,

n_jobs=-1,

verbose=2,

)

gscv.fit(X_train, y_train)

# print results

print('best parameters', gscv.best_params_)

print('best log loss:', gscv.best_score_)

私の環境では1分程度で終了しました。もし実行時間が長くなりすぎるなら、パラメータの組み合わせを少なくしてから試してみてください。

次に、LightGBMのplot_importanceメソッドを使用し、特徴量の重要度を可視化してみます。

当然ですが、ゴールへの距離と角度の重要度が高いですね。

試合の実況を聞いていても、シュート後に「角度がないところから～」「遠いところから～」というフレーズはよく聞きますし、直感にも合致するので納得です。

続いて交差検証を行い、検証データに対する精度を確かめていきましょう。

グリッドサーチで得た最適なパラメータを使いますが、過学習を防ぐため検証データに対してearly stoppingをかけます。

また、log lossを最小化しているので事後確率のcalibrationは十分に行われているとは思いますが、さらなる精度向上を目指しcalibrationを行っておきます。

scikit-learnにはCalibratedClassifierCVクラスが用意されているため、これを使って事後確率の較正を行います。

“isotonic”と”sigmoid”の2つの方法があり、前者はノンパラメトリックな手法、後者はPlatt scalingによって較正をします。

今回は2つとも試してみて、検証データに対するlog lossの値が最も小さいモデルを選ぶことにします。

from sklearn.model_selection import StratifiedKFold
from sklearn.calibration import CalibratedClassifierCV
from sklearn.metrics import log_loss

# cross validation

# initizalize lists contain the results
val_scores = np.zeros((3, 5))
val_diffs = np.zeros((3, 5))
test_scores = np.zeros((3, 5))
test_diffs = np.zeros((3, 5))
y_test_preds = np.zeros((3, 5, len(X_test))) 

# define Stratified KFold
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)

for fold, (train_idx, val_idx) in enumerate(skf.split(X_train, y_train)):
   X_train2 = X_train.iloc[train_idx]
   X_val = X_train.iloc[val_idx]
   y_train2 = y_train.iloc[train_idx]
   y_val = y_train.iloc[val_idx]

   print('Fold {}'.format(fold + 1))
   
   # LightGBM
   lightgbm = lgbm.LGBMClassifier(
                                  objective='binary', 
                                  n_estimators=10000,
                                  random_state=0,
                                  **gscv.best_params_,
                                  )
   lightgbm.fit(X_train2, y_train2,
                     eval_set=[(X_val, y_val)],
                     early_stopping_rounds=20,
                     verbose=-1,
                     ) 
   
   y_val_pred = lightgbm.predict_proba(X_val)[:, 1]
   y_test_pred = lightgbm.predict_proba(X_test)[:, 1]
   val_scores[0, fold] = log_loss(y_val, y_val_pred)
   test_scores[0, fold] = log_loss(y_test, y_test_pred)
   val_diffs[0, fold] = np.sum(y_val_pred) - np.sum(y_val)
   test_diffs[0, fold] = np.sum(y_test_pred) - np.sum(y_test)
   y_test_preds[0, fold] = y_test_pred

   # isotonic calibration
   isotonic = CalibratedClassifierCV(lightgbm, cv='prefit', method='isotonic')
   isotonic.fit(X_val, y_val)
   
   y_val_pred = isotonic.predict_proba(X_val)[:, 1]
   y_test_pred = isotonic.predict_proba(X_test)[:, 1]
   val_scores[1, fold] = log_loss(y_val, y_val_pred)
   test_scores[1, fold] = log_loss(y_test, y_test_pred)
   val_diffs[1, fold] = np.sum(y_val_pred) - np.sum(y_val)
   test_diffs[1, fold] = np.sum(y_test_pred) - np.sum(y_test)
   y_test_preds[1, fold] = y_test_pred

   # sigmoid calibration
   sigmoid = CalibratedClassifierCV(lightgbm, cv='prefit', method='sigmoid')
   sigmoid.fit(X_val, y_val)
   
   y_val_pred = sigmoid.predict_proba(X_val)[:, 1]
   y_test_pred = sigmoid.predict_proba(X_test)[:, 1]
   val_scores[2, fold] = log_loss(y_val, y_val_pred)
   test_scores[2, fold] = log_loss(y_test, y_test_pred)
   val_diffs[2, fold] = np.sum(y_val_pred) - np.sum(y_val)
   test_diffs[2, fold] = np.sum(y_test_pred) - np.sum(y_test)
   y_test_preds[2, fold] = y_test_pred
   print()

# compute the average values between folds
val_scores = np.mean(val_scores, axis=1)
val_diffs = np.mean(val_diffs, axis=1)
test_scores = np.mean(test_scores, axis=1)
test_diffs = np.mean(test_diffs, axis=1)
y_test_preds = np.mean(y_test_preds, axis=1)

# model names
model_names = ['LightGBM',
              'LightGBM + Isotonic',
              'LightGBM + Sigmoid',
              ]

# print the validation resuls 
for i, model_name in enumerate(model_names): 
   print(model_name)
   print('\tlog loss: {:.3f}'.format(val_scores[i]))
   print('\tgoal diff: {:.3f}\n'.format(val_diffs[i]))

from sklearn.model_selection import StratifiedKFold

from sklearn.calibration import CalibratedClassifierCV

from sklearn.metrics import log_loss

# cross validation

# initizalize lists contain the results

val_scores = np.zeros((3, 5))

val_diffs = np.zeros((3, 5))

test_scores = np.zeros((3, 5))

test_diffs = np.zeros((3, 5))

y_test_preds = np.zeros((3, 5, len(X_test)))

# define Stratified KFold

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)

for fold, (train_idx, val_idx) in enumerate(skf.split(X_train, y_train)):

X_train2 = X_train.iloc[train_idx]

X_val = X_train.iloc[val_idx]

y_train2 = y_train.iloc[train_idx]

y_val = y_train.iloc[val_idx]

print('Fold {}'.format(fold + 1))

# LightGBM

lightgbm = lgbm.LGBMClassifier(

objective='binary',

n_estimators=10000,

random_state=0,

**gscv.best_params_,

)

lightgbm.fit(X_train2, y_train2,

eval_set=[(X_val, y_val)],

early_stopping_rounds=20,

verbose=-1,

)

y_val_pred = lightgbm.predict_proba(X_val)[:, 1]

y_test_pred = lightgbm.predict_proba(X_test)[:, 1]

val_scores[0, fold] = log_loss(y_val, y_val_pred)

test_scores[0, fold] = log_loss(y_test, y_test_pred)

val_diffs[0, fold] = np.sum(y_val_pred) - np.sum(y_val)

test_diffs[0, fold] = np.sum(y_test_pred) - np.sum(y_test)

y_test_preds[0, fold] = y_test_pred

# isotonic calibration

isotonic = CalibratedClassifierCV(lightgbm, cv='prefit', method='isotonic')

isotonic.fit(X_val, y_val)

y_val_pred = isotonic.predict_proba(X_val)[:, 1]

y_test_pred = isotonic.predict_proba(X_test)[:, 1]

val_scores[1, fold] = log_loss(y_val, y_val_pred)

test_scores[1, fold] = log_loss(y_test, y_test_pred)

val_diffs[1, fold] = np.sum(y_val_pred) - np.sum(y_val)

test_diffs[1, fold] = np.sum(y_test_pred) - np.sum(y_test)

y_test_preds[1, fold] = y_test_pred

# sigmoid calibration

sigmoid = CalibratedClassifierCV(lightgbm, cv='prefit', method='sigmoid')

sigmoid.fit(X_val, y_val)

y_val_pred = sigmoid.predict_proba(X_val)[:, 1]

y_test_pred = sigmoid.predict_proba(X_test)[:, 1]

val_scores[2, fold] = log_loss(y_val, y_val_pred)

test_scores[2, fold] = log_loss(y_test, y_test_pred)

val_diffs[2, fold] = np.sum(y_val_pred) - np.sum(y_val)

test_diffs[2, fold] = np.sum(y_test_pred) - np.sum(y_test)

y_test_preds[2, fold] = y_test_pred

print()

# compute the average values between folds

val_scores = np.mean(val_scores, axis=1)

val_diffs = np.mean(val_diffs, axis=1)

test_scores = np.mean(test_scores, axis=1)

test_diffs = np.mean(test_diffs, axis=1)

y_test_preds = np.mean(y_test_preds, axis=1)

# model names

model_names = ['LightGBM',

'LightGBM + Isotonic',

'LightGBM + Sigmoid',

]

# print the validation resuls

for i, model_name in enumerate(model_names):

print(model_name)

print('\tlog loss: {:.3f}'.format(val_scores[i]))

print('\tgoal diff: {:.3f}\n'.format(val_diffs[i]))

結果を見ると最もlog lossが低いのはLightGBM + Isotonicのモデルだったので、このモデルをテストデータの予測に使うことにしましょう。

交差検証の各foldにおいてテストデータに対する予測値とスコアを格納しておいたので、平均をとることで最終的な予測値およびスコアとしました。

# specify the model for use 
model_num = 1 # LightGBM + Isotonic
test_score = test_scores[model_num]
test_diff = test_diffs[model_num]
xg = y_test_preds[model_num]

# specify the model for use

model_num = 1 # LightGBM + Isotonic

test_score = test_scores[model_num]

test_diff = test_diffs[model_num]

xg = y_test_preds[model_num]

それでは実際に結果を見ていきましょう。

xGのキャリブレーションプロットとヒストグラムを描画します。

from sklearn.calibration import calibration_curve

# test

# print the test results
print('log loss: {:.3f}'.format(test_score))
print('goal diff: {:.3f}'.format(test_diff))

# plot calibration curve and histgram
fig = plt.figure(figsize=(10, 10))
ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
ax2 = plt.subplot2grid((3, 1), (2, 0))

ax1.plot([0, 1], [0, 1], 'k:', label='Perfectly calibrated')

fraction_of_positives, mean_predicted_value = calibration_curve(y_test, xg, n_bins=10)

ax1.plot(mean_predicted_value, fraction_of_positives, 's-', label='{} ({:.3f})'.format(model_names[model_num], test_score))

ax1.set_ylabel("Conversion rate")
ax1.set_ylim([-0.05, 1.05])
ax1.set_xticks(np.linspace(0, 1, 11))
ax1.set_yticks(np.linspace(0, 1, 11))
ax1.legend(loc="upper left")
ax1.set_title('Calibration plots')

n, bins, patches = ax2.hist(xg, range=(0, 1), bins=10, label=model_names[model_num], histtype='step', lw=2)
ax2.set_xlabel("Expected goals (xG)")
ax2.set_ylabel("Count")
ax2.set_xticks(np.linspace(0, 1, 11))
ax2.legend(loc="upper center", ncol=2)

plt.tight_layout()
plt.show()

from sklearn.calibration import calibration_curve

# test

# print the test results

print('log loss: {:.3f}'.format(test_score))

print('goal diff: {:.3f}'.format(test_diff))

# plot calibration curve and histgram

fig = plt.figure(figsize=(10, 10))

ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)

ax2 = plt.subplot2grid((3, 1), (2, 0))

ax1.plot([0, 1], [0, 1], 'k:', label='Perfectly calibrated')

fraction_of_positives, mean_predicted_value = calibration_curve(y_test, xg, n_bins=10)

ax1.plot(mean_predicted_value, fraction_of_positives, 's-', label='{} ({:.3f})'.format(model_names[model_num], test_score))

ax1.set_ylabel("Conversion rate")

ax1.set_ylim([-0.05, 1.05])

ax1.set_xticks(np.linspace(0, 1, 11))

ax1.set_yticks(np.linspace(0, 1, 11))

ax1.legend(loc="upper left")

ax1.set_title('Calibration plots')

n, bins, patches = ax2.hist(xg, range=(0, 1), bins=10, label=model_names[model_num], histtype='step', lw=2)

ax2.set_xlabel("Expected goals (xG)")

ax2.set_ylabel("Count")

ax2.set_xticks(np.linspace(0, 1, 11))

ax2.legend(loc="upper center", ncol=2)

plt.tight_layout()

plt.show()

log lossの値は0.406、実際のゴール数とxGの合計の差分が1.64点と、まずまずの結果です。

差分が少ない場合でもシュート1本ごとの予測ゴール確率が実際の確率とズレている場合があるので、キャリブレーションプロットで確認します。

キャリブレーションプロットは、予測確率と実際の確率がどれほど一致しているかを示します。点線の上に各プロットが乗っているのが理想的な状態です。

グラフによると、ゴール確率が中程度のシュートに関しては確率のズレがありますが、ゴール確率が低いシュートに関してはよくcalibrateされているようです。

ただ問題の性質上テストデータのサンプルサイズが小さくなるので、なんとも言えない感も…。

pandasの機能を利用し、要約統計量もチェックしてみましょう。

# print summary statistics of xG
pd.Series(xg).describe()

1 2	# print summary statistics of xG pd.Series(xg).describe()

平均値は約0.15となっており、バルセロナのシュートは15%の確率でゴールに繋がることがわかります。

シュートに対するゴールの割合は大体10%と言われているので、バルセロナ選手のシュート能力の非凡さが伺えますね。

最小値が約0.009、最大値が0.98と分布の幅は広いものの、下から数えて75%のxGが0.20以下に収まっていることから、大多数のシュートのゴール確率が低いことがわかります。

最後に、予測したxGをシュートデータに結合しましょう。

# add the 'xg' column to 18/19 season data
shots_1819 = shots.iloc[-num_test:].copy()
shots_1819['xg'] = xg

# add the 'xg' column to 18/19 season data

shots_1819 = shots.iloc[-num_test:].copy()

shots_1819['xg'] = xg

18/19シーズンにおけるバルセロナのxG

前々回の記事で既に15/16シーズンのシュートについて考察を行いましたが、新しく追加された18/19シーズンのデータについても考えてみます。

今回作成した予測モデルにより弾き出されるxGの値は、「過去のメッシ出場試合における、バルサの平均的な選手がシュートを放った場合のゴール確率」と捉えることができるため、xGと実際の結果を比較することにより、選手の評価をすることができます。

まずは、選手ごとの平均xGと決定率の比較をしましょう。

平均xGでトップの値を示したのはスアレスで、質の高い決定機でシュートを放っていたことがわかります。

ただ決定率は低くなっており、もっと点を取れていてもおかしくありません。

結果編で取り上げた15/16シーズンでは高い決定率を示していたので、流石のスアレスも徐々に衰えてきているのかもしれません。

目立つのはメッシとウスマン・デンベレです。平均xGより高い決定率を示しています。

メッシについては言うまでもありませんが、問題児とされてきたデンベレが高い能力を有していることがわかります。

実際に映像を確認してみましたが、角度のない位置からのシュートやミドルを決めており、シュート能力はさすがです。

個人的な評価になりますが、デンベレはボールを持った状態のアクションが異常なほど速く、そこから生まれるシュートの質も高いので、腐らせるのは勿体なく感じています。

次にアシスト期待値（Expected Asissts, xA）と実際のアシスト数を比較してみましょう。

キーパス直後に発生したシュートのxGの値を、アシスト期待値としてキーパスに付与しています。

メッシやば…何なのこの人…。

13アシストはリーグ首位の結果でしたが、期待値でも同等の値を示しており、妥当な結果であったことが分かります。

最後にxGが最も低かったゴールを載せて、本企画を締めさせていただきます！

3位. スアレスの劇的同点弾

試合：18/4/2 ビジャレアル vs バルセロナ
xG : 0.036
再生時間：1分15秒～

La Liga公式コンテンツが含まれているため、ブログ内で視聴できません。YouTubeに飛んでご覧ください。

2位. ラキティッチのボレーミドル

試合：18/10/20 バルセロナ vs セビージャ
xG : 0.034
再生時間：1分25秒～

1位. ラキティッチのジャンピングボレー

試合：18/9/2 バルセロナ vs ウエスカ
xG : 0.028
再生時間：1分8秒～

コードはGitHubに置いています。

データがメッシ出場試合に限られてしまっていること、モデルのチューニング、確率のcalibrationなど、課題は多くありますので、みなさんもxG算出にチャレンジしてみてください。