「卓球競技団体戦における"ダブルスの勝敗"が"団体戦の勝敗"に与える因果効果を定量化したい」ステップ②回帰モデルとバックドア基準を用いて因果効果を評価

始めに

卓球競技の団体戦において、「ダブルスが重要である」「ダブルスが団体戦の鍵をにぎる」と経験知としてよく言われます。
学生時代等に卓球に打ち込んだ方々は、顧問の先生やクラブチームのコーチから言われたことがある人も少なくないのではないでしょうか。
かくいう私も学生時代に卓球に夢中になった者の一人ですが、色んな方々から言われた覚えがあるので、卓球界においてある種の暗黙知として存在するような気がします。

そこで、統計的因果推論を用いて「卓球団体戦におけるダブルスの勝敗」と「団体戦の勝敗」の因果関係の定量分析を試みます。
もう少し統計的因果推論の用語を用いると、本記事の目的は、「卓球団体戦におけるダブルスの勝敗(=原因変数)」が「団体戦の勝敗(=結果変数)」に与える因果効果(介入効果)の定量化、になります。

以前Twitterで、高校時代の先輩が「『卓球の団体戦はダブルスを取った方が有利』って説、誰か定量的に分析してくれないかなー」みたいな事を言っていたような気もするので、まあ需要もゼロではないんじゃないでしょうか？笑

自身の統計的因果推論の練習を兼ねてますので、見識がある方からのツッコミは歓迎しております！
また私同様に、卓球が好きな方からのコメント・ツッコミも歓迎しております！
この記事はステップ①の続きになりますが、「分析パートだけ見ればいいや！」って方は前記事を読まなくても理解に支障ありません！

(統計的)因果推論って何だっけ?

実は、現実世界の観測データから因果関係を捉える事は中々難しいようです。
例えば、観測データから散布図を作成したとします。

しかし散布図だけをいくら眺めていても、どうやら「Xの値が大きいときに、Yの値も大きいという関係がありそうだ」という事が分かっても、「Xの値を大きくした時に、Yの値も大きくなるかどうか」は分かりません。
前者は相関関係(狭義の)、後者が因果関係に対応します。

(統計的)因果推論とは、「観測データからいかに"因果関係"に迫るか」をモチベーションにした方法論、のようです。

以下に、因果推論に関する用語を少し整理します。(飛ばして貰ってもOKです!)

• 相関関係(狭義の):
  • 2つの変数の間に、一方の変数の値が大きいときに多方の変数の値も大きい(小さい)という、直線的な関係がある場合.
• 因果関係：
  • 要因Xを変化させた時に、要因Yも変化する場合、「XとYの間に因果関係がある」という。
  • 要因Xを原因、要因Yをその結果、と呼ぶ.
  • 原因を示す変数＝原因変数、結果を表す変数＝**結果変数(アウトカム)**と呼ぶ.
    (今回の場合は、"ダブルスの勝敗"が原因変数、"団体戦の勝敗"が結果変数になります。)
  • また、原因変数を操作して結果変数を変化させる事を、介入(Intervention)、処置(Treatment)等という。
• 因果効果：
  • 人を対象とした場合には、「同じ人物が、介入を受けた場合の結果変数の値と、介入を受けなかった場合の結果変数の値の差」

交絡因子の存在

因果推論においてとにかく重要な事は、「相関関係は必ずしも因果関係を意味しない」事になります。
変数間に因果関係があれば相関関係も観測されますが、逆は必ずしも真ではありません。

なぜ相関関係は必ずしも因果関係を意味しないのか、その大きな要因は"交絡因子"の存在です。
"交絡因子"に関して、因果ダイアグラム(因果グラフ)と呼ばれる図を使って説明します。

因果関係のない変数XとYに、ある共通の要因Cが影響を与えている(因果関係を持つ)場合を考えます。
この構造において、要因Cが変化した場合、XとYも変化します。
そして仮に分析者がXとYの変化のみを観測していた場合、Xの変化によってYが変化した(もしくはその逆)と誤解してしまう可能性がありますよね...

このような構造を交絡(Confounding)、XとYに影響を与えている共通の要因Cを**交絡因子(Confounding Factor)**と呼びます。
そして、本来は関連のない変数(XとY)の間に、交絡によって生じる相関関係の事を、**疑似相関(Spurious Correlation)**と呼びます。

統計的因果推論は、「交絡因子にいかに対処するか」というのが重要であるようです。

回帰モデルを用いた統計的因果推論

統計的因果推論において最も基本的な手法の1つに、回帰モデルを用いた手法があります。
結果変数yを被説明変数、原因変数xを説明変数として、以下の回帰モデルを構成したとします。

y ～ Norm_y(\mu, \sigma^2) \\
\mu = \beta_0 +\beta_1 x

ここで、もし結果変数と原因変数の両者に影響を与える交絡因子がなければ、xの傾きのパラメータβ_1は結果変数と原因変数の関係の強さ=因果効果をそのまま表していると解釈できます。

しかし仮に両者に影響を与える交絡因子zが存在していた場合、上記の回帰モデルで得られる回帰係数β_1は、因果効果と一致しません。
交絡により疑似相関が発生している可能性があるからです。

この場合は、交絡因子zをモデルの説明変数に追加する事で対処します。

y ～ Norm_y(\mu, \sigma^2) \\
\mu = \beta_0 +\beta_1 x + \beta_2 z

もう一度思い出しておくと、交絡因子zが問題を起こしていたのは、zが変化した場合に、それに伴ってxとyも変化してしまう為に、yが変化したのは、(1)xが変化した為なのか、(2)zが変化してxとyが変化した為にxが変化してyも変化したように見えるだけなのか、の2つに区別が付かない点にあります。

従って、どうにかして交絡因子zの値を固定してしまえばこの問題は解決されます。
重回帰モデルでは、任意の変数の傾きのパラメータβ(偏回帰係数)は他の変数の値を固定した場合にその変数の値を1変化させた時のyの変化量を表している、と解釈できます。
つまり、上式のβ_1は、交絡因子zの値を固定してxの値を変化させた場合に、yがどれだけ変化するか=因果効果の推定値と見なせます。
ちなみに、交絡因子が複数ある場合はそれらを同時に説明変数としてモデルに投入してしまえば良いようです。

更にちなみに本記事では、結果変数(被説明変数)「団体戦の勝敗」が0か1かのBinary変数なので、ロジスティック回帰モデルを用います。ロジスティック回帰モデルは、被説明変数をBernouli分布に従う確率的事象と見なして、そのパラメータλが説明変数に依存して変化すると仮定します。
パラメータλは0~1の値である必要があるので、線形回帰の成分をシグモイド関数に入力して0~1に収まるようにしています。

y～Bernouli_y(\lambda) \\
\lambda = sigmoid(\beta_0 +\beta_1 x + \beta_2 z)

バックドア基準とは？＝どんな説明変数を加えるべき?

バックドア基準ってなんぞや？という事なんですが、意味合いとしては以下の事を判断する基準のようなものです。
「回帰分析において『X=>Yの偏回帰係数』が『X=>Yの因果効果』と一致する（ズレない）為には、どの変数を回帰モデルに加えるべきか/べきではないのか」
バックドア基準について本記事でがっつり話してしまうと、読む方が飽きてしまうかもしれないので、本記事ではこれ以上触れません！
ただ、「交絡因子を説明変数に加える事はバックドア基準を満たす」すなわち「交絡因子を説明変数に加える事で、『説明変数X=>被説明変数Yの偏回帰係数』が『原因変数X=>結果変数Yの因果効果』と一致する」、また「何でもかんでも説明変数に加えてしまうとバックドア基準を満たさなくなってしまい、『偏回帰係数』と『因果効果』が一致しなくなる」という事だけ覚えていれば、本記事の内容は理解できると思います。

いざ、「ダブルスの勝敗」が「団体戦の勝敗」に与える因果効果を捉える！

手法

繰り返しになりますが、本記事の目標は、「卓球団体戦におけるダブルスの勝敗(=原因変数)」が「団体戦の勝敗(=結果変数)」に与える因果効果(介入効果)の定量化、になります。
ただ、三点先取の団体戦において、各試合の勝利が団体戦の勝利にPositiveな効果を与える事はそりゃそうだ、自明といっても良いと思います。恐らくPositiveな因果効果が得られると想定されます。

従って今回は、「ダブルスの勝敗」と「シングルス1の勝敗」が、それぞれ「団体戦の勝敗」に与える因果効果を定量化します。
そして両者の因果効果の大きさを比較する事で、「団体戦においてダブルスが鍵なのか否か」を議論してみたいと思います。

交絡因子としては、「ホームチームの地力」と「アウェイチームの地力」を考えました。
当然、自チームに強い選手(ex. 世界ランク上位、中国選手, etc.)がたくさんいれば、各試合の勝率もあがるし、団体戦の勝率も上がると考えられます。
(ちなみにここ10年、団体戦において中国チームが他国チームに敗戦しているのを、男女ともに見た事がありません...!)
また、相手チームの方に強い選手がたくさんいればその逆もしかり、と言えると思います。
よって今回検討する因果ダイアグラムは、以下のようになります。

具体的には「ホームチームの地力」「アウェイチームの地力」の代理指標として、「ホームチームのSランク選手割合」「アウェイチームのSランク選手割合」を用います。

回帰分析の準備

まあ準備といっても、今回はステップ1の段階でダミー変数化もしていますし、交絡因子として用いる「Sランク選手割合」のカラムも構築済みなので、単にステップ1データをロードするのみです。

def load_matchData()->pd.DataFrame:
    '''試合データをpd.DataFrameで読み込み'''

    INPUT_DIR = r'保存場所\Tleague_scraping'
    datasets_list = []

    season_years = ["2019", "2020", "2021"]
    # 各シーズンのデータを読み込み
    for i, season_year in enumerate(season_years):
        datasets_list.append(pd.read_csv(os.path.join(INPUT_DIR, f'{season_year}_match_result.csv'), header=0, encoding='utf-8'))
        
        # シーズンカラムを付与する
        datasets_list[i]["season_year"] = int(season_year)

    # 縦に結合
    df_dataset = pd.concat(datasets_list)

    return df_dataset

df = load_matchData()
df.head()

読み込んだ結果は以下です。このデータセットを使って因果効果を分析していきます。

db	s1	s2	s3	vm	HomeTeam_name	AwayTeam_name	db_HomeWin	s1_HomeWin	s2_HomeWin	s3_HomeWin	vm_HomeWin	HomeTeam_win	HomeTeam_rate_AAA_and_S	AwayTeam_rate_AAA_and_S	HomeTeam_rate_S	AwayTeam_rate_S
12	30	31	23	10	木下マイスター東京	岡山リベッツ	0	1	1	0	1	1	0.5555555555555556	0.5	0.4444444444444444	0.08333333333333333
12	32	31	13	01	日本生命レッドエルフ	木下アビエル神奈川	0	1	1	0	0	0	0.6363636363636364	0.5	0.2727272727272727	0.2
20	32	23	32	00	T.T彩たま	琉球アスティーダ	1	1	0	1	0	1	0.45454545454545453	0.5	0.18181818181818182	0.1
02	31	30	23	10	岡山リベッツ	T.T彩たま	0	1	1	0	1	1	0.5	0.45454545454545453	0.08333333333333333	0.18181818181818182
21	31	13	30	00	木下アビエル神奈川	トップおとめピンポンズ名古屋	1	1	0	1	0	1	0.5	0.4166666666666667	0.2	0.16666666666666666
20	31	13	30	00	木下マイスター東京	琉球アスティーダ	1	1	0	1	0	1	0.5555555555555556	0.5	0.4444444444444444	0.1

各レコードは、1つの団体戦における各試合結果となっています。
各カラムの定義は、以下の様になっています。

• db, s1, s2, s3, vm (str): その団体戦におけるダブルス、シングルス1, シングルス2, シングルス3, ビクトリーマッチのスコア。ex)db="02"の場合、ゲームカウント0-2でアウェイチームの勝利。
• HomeTeam_name (str)：その団体戦におけるホームチーム名です。
• AwayTeam_name (str)：その団体戦におけるアウェイチーム名です。
• db_HomeWin, s1_HomeWin, s2_HomeWin, s3_HomeWin, vm_HomeWin : その団体戦内のダブルス、シングルス1, シングルス2, シングルス3, ビクトリーマッチにおける、ホームチームの勝利を示すダミー変数。
• HomeTeam_win：その団体戦における、ホームチームの勝利を示すダミー変数。
• HomeTeam_rate_AAA_and_S, AwayTeam_rate_AAA_and_S：ホームチーム、もしくはアウェイチームにおける、全選手に占めるAAAランクorSランク選手の割合。
• HomeTeam_rate_S, AwayTeam_rate_S：ホームチーム、もしくはアウェイチームにおける、全選手に占めるSランク選手の割合。（詳しくは次回の記事で説明しますが、今回の因果推論で"交絡因子"として使用します）

回帰モデルによる因果効果の推定結果

(交絡因子なしVer.)「ダブルスの勝敗」と「シングルス1の勝敗」の因果効果の比較

まず交絡因子考慮なしVer.で、「ダブルスの勝敗」を説明変数(+定数成分)、「団体戦の勝敗」を被説明変数としてロジスティック回帰モデルを構築してみます。

import statsmodels.formula.api as smf
import statsmodels.api as sm

result = smf.glm(formula="HomeTeam_win ~ db_HomeWin", data=df, family=sm.families.Binomial()).fit()
result.summary()

得られた回帰係数は以下です。
(**ちなみに回帰係数の推定値が"coef"カラムになります。この回帰モデルの推定値の表は、これ以降本記事で良く出てきますが、"coef"カラムだけ見ていれば十分です！見なくても十分です！**笑)
(一応、他のカラムは標準誤差、z統計量、z検定の結果のp値、等です。Interceptレコードは定数成分です。)

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.8650	0.188	-4.589	0.000	-1.234	-0.496
db_HomeWin	1.8635	0.273	6.822	0.000	1.328	2.399

続いて交絡因子考慮なしVer.で、「シングルス1の勝敗」を説明変数(+定数成分)、「団体戦の勝敗」を被説明変数としてロジスティック回帰モデルを構築してみます。

result = smf.glm(formula="HomeTeam_win ~ s1_HomeWin", data=df, family=sm.families.Binomial()).fit()
result.summary()

得られた回帰係数は以下です。

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.9217	0.194	-4.742	0.000	-1.303	-0.541
s1_HomeWin	1.8957	0.274	6.922	0.000	1.359	2.433

両回帰モデルで得られた回帰係数を比較すると、「ダブルスの勝敗」は1.86、「シングルス1の勝敗」は1.90となり、やや「シングルス1の勝敗」の方が「団体戦の勝敗」に与える因果効果が大きい結果となりました。
しかし前述した通り、これらは交絡因子を考慮しないVer.の推定結果です。ここで推定された因果効果(回帰係数)には恐らく、交絡因子によって発生した疑似相関の影響が含まれており、正確な因果効果ではない、と考えられます。

従って、因果効果をより正確に推定する為に、交絡因子を考慮したVer.の回帰モデルを構築します。

(交絡因子ありVer.)「ダブルスの勝敗」と「シングルス1の勝敗」の因果効果の比較

上記(交絡因子なしVer.)の結果を踏まえ、因果効果をより正確に推定する為に、交絡因子として「両チームの地力」を考慮したVer.の回帰モデルを構築します。
まず「ダブルスの勝敗」「ホームチームのSランク選手割合」「アウェイチームのSランク選手割合」を説明変数(+定数成分)、「団体戦の勝敗」を被説明変数としてロジスティック回帰モデルを構築してみます。

result = smf.glm(formula="HomeTeam_win ~ db_HomeWin + HomeTeam_rate_S + AwayTeam_rate_S", data=df, family=sm.families.Binomial()).fit()
result.summary()

得られた回帰係数は以下です。

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.9394	0.435	-2.160	0.031	-1.792	-0.087
db_HomeWin	1.9088	0.287	6.648	0.000	1.346	2.472
HomeTeam_rate_S	3.5559	1.054	3.374	0.001	1.490	5.621
AwayTeam_rate_S	-3.4833	1.060	-3.285	0.001	-5.561	-1.405

続いて「シングルス1の勝敗」「ホームチームのSランク選手割合」「アウェイチームのSランク選手割合」を説明変数(+定数成分)、「団体戦の勝敗」を被説明変数としてロジスティック回帰モデルを構築してみます。

result = smf.glm(formula="HomeTeam_win ~ s1_HomeWin + HomeTeam_rate_S + AwayTeam_rate_S", data=df, family=sm.families.Binomial()).fit()
result.summary()

得られた回帰係数は以下です。

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.9206	0.433	-2.127	0.033	-1.769	-0.072
s1_HomeWin	1.7802	0.281	6.342	0.000	1.230	2.330
HomeTeam_rate_S	2.7796	1.015	2.739	0.006	0.790	4.769
AwayTeam_rate_S	-2.6679	1.062	-2.513	0.012	-4.749	-0.587

結果として、両回帰モデルで得られた回帰係数を比較すると、「ダブルスの勝敗」は1.91、「シングルス1の勝敗」は1.78となり、「ダブルスの勝敗」の方が「団体戦の勝敗」に与える因果効果が大きい結果となりました。
(実際には、ロジスティック回帰における回帰係数は、そのままでは”因果効果の定義”とは一致しません。)
(通常の一般線形回帰と異なり、シグモイド関数を通して変換してるので...。)
(でも少なくとも因果効果の大小関係は回帰係数と一致するはずです。)
(他に疑似相関を与えうる交絡因子がなければですが...)

また、交絡因子として「両チームの地力」を説明変数に加えた事で、「シングルス1の勝敗」による因果効果の方が「ダブルスの勝敗」よりも大きく低下しました。
これは「両チームの地力(ホームチームのSランク選手割合, アウェイチームのSランク選手割合)」が、「シングルス1の勝敗」と「団体戦の勝敗」の間に発生させる疑似相関がより大きかった、と解釈できます。

言い換えれば「両チームの地力」は、「ダブルスの勝敗」と「団体戦の勝敗」の関係よりも、「シングルス1の勝敗」と「団体戦の勝敗」の関係に対して、より強く疑似相関を発生させる。
すなわち、「両チームの地力」が「団体戦の勝敗」に与える因果効果は一定とみなすと、「ダブルスの勝敗」の方が「シングルス1の勝敗」よりも「両チームの地力」の影響を受けづらい、と解釈できるのではないでしょうか？
(この「ダブルスの方が、強い選手にもワンチャンスある」ともいえる解釈は、卓球選手としても肌感覚・暗黙知にも合致するような気がします！)

ちなみに...「シングルス2の勝敗」「シングルス3の勝敗」による因果効果とも比較してみた！

同様に「シングルス2の勝敗」及び「シングルス3の勝敗」においても、交絡因子として「両チームの地力」を考慮したVer.の回帰モデルを構築し、因果効果の定量化・比較してみました。

まずは「シングルス2の勝敗」。

result = smf.glm(formula="HomeTeam_win ~ s2_HomeWin + HomeTeam_rate_S + AwayTeam_rate_S", data=df, family=sm.families.Binomial()).fit()
result.summary()

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.6240	0.410	-1.523	0.128	-1.427	0.179
s2_HomeWin	1.4209	0.272	5.222	0.000	0.888	1.954
HomeTeam_rate_S	2.8155	1.014	2.776	0.006	0.828	4.803
AwayTeam_rate_S	-2.9341	1.021	-2.873	0.004	-4.935	-0.933

続いて「シングルス3の勝敗」です。

result = smf.glm(formula="HomeTeam_win ~ s3_HomeWin + HomeTeam_rate_S + AwayTeam_rate_S", data=df, family=sm.families.Binomial()).fit()
result.summary()

\	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.4921	0.406	-1.211	0.226	-1.289	0.304
s3_HomeWin	1.1983	0.272	4.399	0.000	0.664	1.732
HomeTeam_rate_S	3.5263	0.992	3.556	0.000	1.583	5.470
AwayTeam_rate_S	-4.0126	1.035	-3.878	0.000	-6.041	-1.98

結果として、各回帰モデルで得られた回帰係数を比較すると、「ダブルスの勝敗」は1.91、「シングルス1の勝敗」は1.78、「シングルス2の勝敗」は1.42、「シングルス3の勝敗」は1.20となり、団体戦の前半の試合ほど因果効果が大きく、後半の試合ほど因果効果が小さい結果が得られました。

この結果から得られる解釈としては、

• 同じシングルスでも、団体戦の前半か後半かで因果効果の大きさが異なる。
• 従って、今回得られた「ダブルスの勝敗」と「シングルス1の勝敗」が「団体戦の勝敗」に与える因果効果の差は、「ダブルスかシングルスか」という試合形式の違いだけではなく、「団体戦の試合の順番」の違いが影響している可能性がある。
• この現象は、「ダブルスが団体戦の1試合目に設定されているデータ」のみを抽出した事で発生した疑似相関、いわゆる「選択バイアス(Selection Bias)」の問題が原因？

もちろんTリーグにおいては、全ての団体戦においてダブルスが1試合目に行われているので、Tリーグに限定していえば、「ダブルスの勝敗」が「団体戦の勝敗」に与える因果効果はある程度、定量評価できたように考えられます。経験知通り、少なくともTリーグにおいては、他の試合と比較して「ダブルスの勝敗」が「団体戦の勝敗」に与える因果効果が最も大きく、「ダブルスが団体戦の鍵をにぎる」と言えそうです。

まとめ

今回はステップ2として、スクレイピングして作成したデータセットを元に、回帰モデルを用いた統計的因果推論によって、「卓球団体戦におけるダブルスの勝敗」と「団体戦の勝敗」の因果関係の定量分析を試みました。
結果として、「ダブルスの勝敗」の方が「シングルス1の勝敗」と比較して「団体戦の勝敗」に与える因果効果が大きい結果となりました。
この結果は「両チームの地力」を交絡因子として考慮した事で明らかになりました。
また、この交絡因子の考慮の有無によって「シングルス1の勝敗」による因果効果(回帰係数)がより大きく変動したことから、「ダブルスの方が、両チームの地力の差による影響が小さい」、言い換えれば、「ダブルスの方が、強い選手に対してワンチャンスある」という、卓球界における肌感覚・経験知に合致する解釈も得られました。
一方で今回の結果には、「ダブルスかシングルスか」という試合形式の違いだけではなく、「団体戦の試合の順番」の違いが影響している可能性も示唆されました。

もちろん、あくまで「上記のアプローチを採用して分析した結果、こういう解釈・示唆が得られた」という事で、今回の因果推論の結果が絶対的に正確であるという事は全くありません。
今回採用した手法も長所短所がありますし、本記事で考慮しきれなかった交絡因子も存在している可能性もあります。
統計的因果推論にも様々な手法があり、色んな分析手法・色んなデータセット等を用いた分析で、同様の結果・解釈が得られれば、より「ダブルスが団体戦の鍵をにぎるか否か」が明らかになっていくのかなと思います：）

ただ個人的には今回の内容にはとても満足していて、交絡因子として「両チームの地力」の考慮の有無によって、「ダブルスの勝敗」と「シングルス1の勝敗」の因果効果の大小関係が逆転したので、そこは工夫してみて良かったな～と思っています笑
統計的因果推論の良い練習になったような気がします！

次回以降はまだ未定ですが、"回帰モデルによる因果推論"の発展として階層型の回帰モデルを用いて、「ダブルスの勝敗」の因果効果の「Tリーグの各チームにおける個体差」みたいなところを定量化してみたいなと思っています。
(このチームは特にダブルス大事！このチームはダブルスの重要性低い、みたいな事が定量評価できれば、より意思決定に繋がる分析になる気がします：）)
その後は、他の統計的因果推論の手法を試してみたいですね！(傾向スコアマッチングとか！機械学習を使った手法とか！)

2回目の投稿という事で、文章的にもMarkdownの構成的にも、稚拙な点が多々あるとは思いますが、最後までお読みいただきありがとうございました！

繰り返しになりますが自身の統計的因果推論やコーディングの練習を兼ねてますので、上記分野が好きな方々からツッコミ・コメントをいただけたらとても喜びます！
また私同様に、卓球が好きな方からのコメント・ツッコミもいただけたらとても嬉しいです！

参考文献

• 岩波データサイエンス vol.3