もし一般企業のデータ分析者がはなおでんがんの「最強の勝負めし分析」企画をやったら
これまでの記事では因果推論の真面目な話を書いてきましたが、今回は趣向を変えてポップなパロディ企画をやります。
記事の概要
ひとことで言うと、タイトルの通り、はなおでんがんさん（以下「はなでんさん」と略します）の動画 【全200試合】藤井七段の勝負メシの勝率から、最強の勝負メシを決めようではないか。 の検証を、データ分析者の立場からより厳密に行ってみようというものです。
はなでんさんのファンとしては、一度動画をご覧になっていただいた上で本記事を読んでいただけると嬉しいですが、動画の背景知識なしでも理解できるように書く予定です。
対象読者、読後イメージ
以下の3タイプの人を対象にします:
- 分析結果だけ知りたい人: 分析の全体まとめだけ読んでいただけると、そういう感じなのか〜という理解に至るかと思います。
- 分析の中身を知りたい人: PythonとRで可視化、検定、回帰する例を載せているので、もしよければ参考にしていただければと思います。
- 因果推論ガチ勢: IPWのところあたりを楽しんでいただければと思います。もし間違っていたら教えていただきたいです*1。
目次は以下になっています:
記事の構成
本記事では、分析結果を共有するときのフォーマットとして、YWTを拡張したものを使います。
まず、分析のフェーズに合わせていくつか大きな段落を作ります。
次に、それぞれの段落の中で、レポートと実装を分離します。レポートにしか興味のない人が余計なものまで見るコストを減らすためです。
最後に、レポートの中身を以下のように構造化します:
- B (Background): 背景を書きます。資料の最初の段落だけに書くことが多いです。
- P (Purpose): 段落の目的をシンプルに書きます。
- H (Hypothesis): 事前に持っている仮説をふわっと書いておきます。
- C (Conclusion): 結論や考察を書きます。次の段落やその後の意思決定に値する情報をシンプルに書きます。
- Y (やったこと): やったことを書きます。
- W (わかったこと): わかったことを事実ベースで書きます。
- T (次やること): 次にやること、やり残したことを書きます。Cと少し被るので厳密にはわけられないですが、Cよりは緩く書きます。
レポート部分は常体、レポート以外は敬体で記事を書く予定です。
分析の全体まとめ
分析結果をまとめます。忙しい人はここだけを見ていただければ大丈夫です。
B
- 冒頭に述べたはなでんさんの動画では、藤井聡太選手の年度ごとのデータを手動で集計し、勝率とご飯の関連を調べることで、「勝負めしとして何が適切か」を調べようとしていた
- 分析結果として「冷えたご飯はよくない」「カレーうどんと天丼がいい」という主張がされていた
- この分析は公開データの分析という点で興味深い試みであるが、以下3点で課題があると考えられる:
- 実際、はなでんさんの動画のコメント欄にも交絡について指摘をする人が存在していたことからも、因果についての分析を行うニーズはあるということが期待される
P
- 手動集計を自動化することにより再現性や拡張性を保証し、因果推論の手法を用いて交絡調整を行うことで因果関係についての主張を行う
- ただし、あくまでこれは動画企画のパロディ・ビッグバネイト・茶番であり、はなでんさんの動画を批判したり否定したりする意図は一切ないということに留意していただきたい
- また、本分析結果はあくまで現時点での一例であり、今後の追加分析によって主張がアップデートされる可能性があることにも留意していただきたい
H
- はなでんさんの分析結果には妥当性がありそうだが、何らかの交絡があるのではないか？
- 本当に検定結果が有意になるのか？
C
- 「冷たいご飯」や「カレーうどん」と勝率の間に一定の傾向はあるが、有意な関連ではなかった（「冷たいご飯だからダメ」と考えずに、冷たいものが食べたくなったら冷たいものを食べてもまあ統計的には大丈夫でしょう！）
- むしろ「うどん」や「そば」の方が関連が大きい
- 交絡調整をすると、「うどん or そば」or 「それ以外」の2水準に処置をわけたときの処置効果は有意ではなくなる（交絡調整しないときは有意になるので注意！）
- 交絡調整をしても、「うどん」「そば」「それ以外」の3水準に処置をわけたとき、「そば」を食べることの効果は有意にゼロより大きい（勝ちたいときはそばを食べてみよう！）
Y
W
- Cと同じ
T
- 藤井聡太さん以外のデータを集めてより一般性のある主張をする
- トレンドや場所などの情報を考慮してより正確な分析をする
基礎集計と分析設計
基礎集計と分析設計について詳しく説明します。
レポート
P
H
- はなでんさんの動画で参考にしていたページからスクレイピングすれば、食事名と勝敗のデータを簡単に取得できるのではないか
C
Y
- スクレイピング
- 対局一覧を https://shogi-sanpo.com/fujiisouta/result/ から取得し、各対局の詳細データをさらにスクレイピングして "fujii_shogi_meal.csv" として保存
- ご飯一覧を https://shogi-sanpo.com/fujiisouta/gourmet/ から取得し、"fujii_shogi_meal2.csv" として保存
- 不整合調査
- 2つのcsvと実際の対局記録から、データのズレをチェック
- 前処理
- 対局が2日にまたがった2件はリストワイズで除去（スクレイピングや集計が面倒）
- 店の名前はいったんは無視（分析目的からずれるのと、データに揺れがあるので）
- 1回目の食事名をfirst_mealとしてtreatmentにする（データ量の観点と、意味的にも妥当そう）
- 「不明」というデータが1つだけあったので欠損に置き換え
- 変数の分布理解
- 各変数の欠損パターンを把握
- outcomeとして、勝敗数の推移を把握
- treatmentとして、ご飯名の出現頻度を把握
W
-
- ご飯一覧のページは途中からデータの定義が変わっており、ある時期以前は詳細データのリンクまでたどることができないため、対局相手の情報を取得できない
- 対局一覧のページではすべての詳細データのリンクをたどることができるので、対局相手の情報はここから取得する必要がある
不整合調査
対局一覧がおかしいところが2点
- 20170810: 昼食の店の名前が間違っている（https://www.nikkansports.com/general/nikkan/news/1869967.html）
- 20181123: 夕食が昼食になっている（https://shogi-sanpo.com/fujiisouta/kisen/eiou/20181123/）
ご飯一覧がおかしいところが1点
- 20171102: 昼食の店の名前が間違っている（https://www.sankei.com/west/photos/171102/wst1711020048-p1.html）
どちらがおかしいか不明なところが1点
- 20180125: 夕食の味噌煮込みうどんがご飯一覧にのみ存在する
それ以外は店名ベースで見ると表記揺れ以外問題がなかった
変数の分布理解
- だいたいどの年も勝率は80%overをキープしている
- 2016, 2020年度は試合数が少ない
- 1回目の食事情報が含まれているのは148/218件で、67%ほど
- 食事名については対局一覧の方が表記揺れが多い
T
- 分析に入る
- 藤井聡太さん以外のデータを集める
- 「1回目の食事情報」よりもう少し深い情報（例えば夕食とおやつの考慮）を集める
- データをもう少し綺麗にする
実装
Jupyter notebook 上でスクレイピングと基礎集計を行いました。Python3.6です。
スクレイピングのコードやデータは公開できないです（スクレイピングの処理が一番大変でした）が、基礎集計部分は気が向いたらuploadします。
主要な処理のみ記事の中で紹介します。
importはこちらです:
%matplotlib inline import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from IPython.display import display_html import researchpy
スクレイピングしてきた対局とご飯のデータはこのようになっています:
df = pd.read_csv("fujii_shogi_meal.csv")
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 218 entries, 0 to 217 Data columns (total 12 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 date 218 non-null object 1 title 218 non-null object 2 link 218 non-null object 3 description 218 non-null object 4 is_play_first 218 non-null bool 5 lunch_store 145 non-null object 6 lunch 146 non-null object 7 dinner_store 68 non-null object 8 dinner 70 non-null object 9 ct_win 215 non-null float64 10 ct_lose 215 non-null float64 11 win 218 non-null int64 dtypes: bool(1), float64(2), int64(1), object(8) memory usage: 19.1+ KB
218件なので多いとは言えないですが、一応集計が可能なレベルかとは思います。
変数の定義は以下です:
- date: 対局日付がyyyy/mm/ddで入っている - title: 対局名が入っている。竜王戦とか。 - link: 対局詳細のリンク。 - description: 対局の説明。 - is_play_first: 先行後攻のフラグ - lunch_store: 昼食をとった店の名前 - lunch: 昼食メニュー名 - dinner_store: 夕食をとった店の名前 - dinner: 夕食メニュー名 - ct_win: 相手の通算勝利数（対戦日まで） - ct_win: 相手の通算敗北数（対戦日まで） - win: 対局の勝敗（1が勝利、0が敗北）
データの不整合調査は泥臭い作業なので記事では省略します。
年度ごとの勝率の遷移はこのようになりました。めちゃくちゃ強いですね。年ごとのトレンドを考慮する必要性は（今のところ）それほどなさそうです。
df["date"] = pd.to_datetime(df["date"]) date_q = pd.PeriodIndex(df["date"], freq="Q") df["fiscal_year"] = date_q.map(lambda x: x.year - 1 if x.quarter == 1 else x.year) sns.pointplot(x="fiscal_year", y="win", data=df) plt.ylim((0, 1)) plt.title("勝率は80%以上をキープしている") plt.ylabel("勝率")
ちなみに、対局数は2016, 2020年が少ないので、解釈の際には注意してください:
df.groupby("fiscal_year")["win"].count().plot.bar() plt.title("2016, 2020年度は対局数が少ない") plt.ylabel("対局数")
対局一覧からご飯情報を集計してみると、以下のことがわかります:
- 食事情報が含まれているのは147/218回（うち1回だけ不明）で、67%ほど
- lunchなしdinnerありは1件しかない（ただしデータ不整合で確認したようにこの分類は必ずしも正確ではない）
- lunchがありのときにdinnerがあるかどうかは半分ずつくらいの割合
pd.crosstab(df.lunch.notnull(), df.dinner.notnull())
dinner False True lunch False 71 1 True 77 69
昼食夕食の揺れを排除するため、1回目の食事名に相当する変数を作ります:
df["first_meal"] = df["lunch"] df.loc[df["lunch"].isnull(), "first_meal"] = df["dinner"] df.loc[df.first_meal == "不明", "first_meal"] = np.nan
食事名の頻度を集計するとこのようになります:
plt.figure(figsize=(6, 10)) cnt = df.first_meal.value_counts() remove_cnt = sum(cnt == 1) cnt.loc[lambda x: x > 1][::-1].plot.barh() plt.xlabel("メニュー") plt.ylabel("注文回数") plt.title(f"1回目の食事の注文回数(1回しか注文されていない {remove_cnt} 種類のメニューは除外)")
しかし、実際に食事名で並び替えて一覧を見てみると、かなりの表記揺れ（大文字小文字、スペース有無 etc）が存在していることがわかりました。
そのため、食事名そのものではなく、表記揺れに気をつけつつざっと見て目立ったキーワードでデータを作り直す必要があると判断しました。
ご飯と勝率の分析
方針が立ったので、実際に関連を分析していきます。
やるだけなので、この段落はすごくシンプルです。
レポート
P
- キーワードをいくつか出して、勝率との関連を調べる
H
- 冷たいご飯は勝率にネガテイブな関連を示すのではないか
C
- 冷たいご飯ではなく、「うどん」や「そば」が有意な関連を示したが、因果的な解釈をできるかを追加で調べる必要がある
Y
- 冷たい食事フラグ、カレーうどんフラグを作って、勝率を検定
- 生データを見て頻繁に出てきたうどんとそばについての変数を作って、勝率を検定
W
- はなでん動画の仮説に反して、冷たいご飯やカレーうどんは、勝率との有意な関連はない
- 「うどん」と「そば」は勝率に対してポジティブな関連がないとは言えない
T
- キーワードの候補を増やす
- 交絡を調整しに行く
実装
先ほどと同様のJupyter notebook上で実装しました。
処置変数の候補を作ります:
df["cold_flag"] = df.first_meal.str.contains("冷") df["curry_udon_flag"] = (df.first_meal.str.contains("カレー")) & ((df.first_meal.str.contains("うどん")) ) df["hamburger_flag"] = df.first_meal.str.contains("ハンバーグ") df["udon_soba_else"] = df.first_meal.map( lambda x: x if x != x else "both" if "うどん" in x and "そば" in x else "udon" if "うどん" in x else "soba" if "そば" in x else "else") df["udon_or_soba"] = df.udon_soba_else.map(lambda x: x if x != x else x in ["udon", "soba"])
はなでんさんの動画で出ていた「冷たいご飯」の他に、目立ったキーワードをいくつか加えました。 カレーうどんとハンバーグは処置群が少なかったので分析しません。
さて、冷たいご飯について検定してみます:
crosstab, res = researchpy.crosstab(df.cold_flag, df.win, test= "chi-square", correction=True) print("分割表") display_html(crosstab) print("検定結果") display_html(res)
分割表 win win 0 1 All cold_flag False 16 103 119 True 6 21 27 All 22 124 146 検定結果 Chi-square test results 0 Pearson Chi-square ( 1.0) = 0.7276 1 p-value = 0.3937 2 Cramer's phi = 0.0706
ブログの都合上見にくくなっていますが、「冷たいご飯か否か」は勝率との関連は (5%) 有意でないという結果になりました。
分割表についての他の検定結果もすべて同様でした。
効果量も小さい（< 0.1）ので、この仮説は微妙そうだ、と言えそうです。
次に、「うどん」「そば」「その他」の3水準にわけた変数について調べてみます。
crosstab, res = researchpy.crosstab(df.udon_soba_else, df.win, test= "chi-square", correction=True) print("分割表") display_html(crosstab) print("検定結果") display_html(res)
分割表 win win 0 1 All udon_soba_else else 17 64 81 soba 0 22 22 udon 5 38 43 All 22 124 146 検定結果 Chi-square test results 0 Pearson Chi-square ( 2.0) = 6.5185 1 p-value = 0.0384 2 Cramer's V = 0.2113
この検定によって、udon_sobe_elseの水準は5%有意、つまり「関連がないとは言えない」ことがわかりました。
分割表についての他の検定結果もすべて同様でした。
効果量は中程度より少し弱め（0.21 ~ 0.25）ので、この仮説は比較的筋が良さそうだ、と言えそうです。
なお、「うどんかそば」と「それ以外」に水準をわけてみたところ、同様に「関連がないとは言えない」ことがわかりました。
次の段落ではこれらの処置についての検証を深ぼっていきます。
天気と相手の勝率による交絡調整
ようやく因果推論の話になります。天気と相手の勝率で交絡調整をしていきます。
レポート
B
- これまでの分析では、以下のような交絡が考慮されていない
- 「暑い日は冷たいものを注文しやすく、暑い日は負けやすいが、冷たいものを食べるかどうかと負けるかどうかには因果関係はない」という状況が真だったときに、「冷たいもの」と「勝敗」だけの関連を分析すると、「冷たいものを食べると負けやすい」といった誤った解釈に至る可能性がある
- このような、交絡を引き起こす「暑さ」という変数は交絡変数と呼ばれ、分析において調整を行う必要がある
P
- 交絡調整を行って分析を進め、因果的な解釈を行う。
H
- 天気、相手の勝率、体調、盤面の様子、気分 などが交絡変数になりそう（年度や場所の効果は2で微妙だったので行わない）
C
- 天気と相手の勝率の変数が交絡調整に使えると判断した
- 交絡調整をすると「うどん or そば」or 「それ以外」の処置効果は有意ではなくなる
- 交絡調整をしても「そば」を食べることの効果は有意にゼロより大きく、ポジティブと言えるであろう（勝ちたいときはそばを食べてみよう！）
Y
- 藤井聡太選手のインタビュー*4をもとに、食べるご飯をどのように決めているのかを調べて、交絡変数をリストアップした
- 重要そうな変数として「天気」と「相手の勝率」の2つを選択し、データを整形した
- 天気
- 多くの食事が昼間にとられるものなので、最高気温を使う
- 対局場所が簡単に見つからなかったので、東京のデータを気象庁から入手して代用した (http://www.data.jma.go.jp/gmd/risk/obsdl/index.php)
- 勝率
- 対戦相手の対戦日前までの通算勝率として表示されているデータから計算
- 1回しか戦績がないケースで勝率が100%になっている例は外れ値すぎるので、flat priorを使って補正（勝率の性質を考えると、真ん中らへんが下がっているJeffreys priorよりは妥当と判断）
- 天気
- 2つの変数でいくつかのtreatmentのmodelを作り、IPW によって交絡を調整した
W
- インタビュー記事を見ると、即決で食事を決めているとの発言がある
- 気温が高いと冷たいご飯を食べやすい
- 気温が低いとそばを頼みやすい
- 相手の勝率と冷たいご飯の関連はほとんどない
- 相手の勝率が高いときは「うどんやそば」以外を頼みやすい
- 気温の高さと藤井聡太さんの勝率には負の関連がある
- 相手の勝率と藤井聡太さんの勝率には負の関連がある
- 交絡調整の有無に関わらず、そばの処置効果だけが有意にゼロでなない
- 交絡調整をすると、「うどん or そば」の水準を一つにまとめた処置効果が有意でなくなる
- 交絡調整によってうどんやそばの回帰係数は小さくなった
T
- さらに交絡調整を行う。例えばご飯を食べる直前の盤面の有利不利のスコアなど
- 年度のトレンドをうまく反映して分析する
実装
前処理は先ほどと同じ Jupyter notebookで行い、回帰分析からは Rmarkdownで行いました。
Pythonでの処理
天気の変数は、気象庁から落としてきたデータから以下のように加工して結果を確認しました:
weather = pd.read_csv( "~/Downloads/temperature.csv", skiprows=[0,1,2,4,5], encoding='shift_jis') weather["date"] = pd.to_datetime(weather["年月日"]) weather["max_temperature"] = weather["最高気温(℃)"] plt.plot(weather.date, weather["max_temperature"])
:
特におかしな挙動はしていないですね。
相手の勝率は以下のように加工しました:
alpha = 1 beta = 1 df["ct_win_ratio"] = (df.ct_win + alpha) / (df.ct_win + df.ct_lose + alpha + beta)
joinして保存しました:
df_merge = df.merge(weather[["date", "max_temperature"]], on="date") df_merge[df_merge.cold_flag.notnull()].to_csv("~/Downloads/fujii_merge.csv", index=False) df_merge.to_csv("~/Downloads/fujii_merge_all.csv", index=False)
joinしたデータで相関係数を見ると、以下のことがわかります:
fig = plt.figure(figsize=(10, 7)) sns.heatmap( df_merge[["max_temperature", "ct_win_ratio", "win", "udon_or_soba", "cold_flag"]].astype(float).corr(), annot=True, fmt="1.2f" )
Rでの基礎集計
ここから先の分析は、Rを用いて行います。
あまりRに慣れていないので雰囲気で書いていきます。
まずは因果的な分析を行わないnaiveな集計と可視化を行います。
利用するライブラリを読み込みます:
library(ipw) library(tidyverse)
データを加工します:
df <- read.csv("Downloads/fujii_merge.csv") # 相手の勝率がない場合は0.5で埋める（1件しかない） df$ct_win_ratio[is.na(df$ct_win_ratio)] <- 0.5 # Rのipwでlogitモデルを使うときはtreatmentを0-1に変換しないとダメ df_cp <- df df_cp$cold_flag <- as.integer(df_cp$cold_flag) - 1
冷たいご飯のフラグについてnaiveな回帰をしてみると、「係数の推定値は負だが5%有意ではない」という先ほどと同様の傾向が得られました:
summary(lm(win ~ cold_flag, data=df_cp))
Residuals:
Min 1Q Median 3Q Max
-0.8656 0.1344 0.1344 0.1344 0.2222
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.86555 0.03287 26.332 <2e-16 ***
cold_flag -0.08777 0.07644 -1.148 0.253
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3586 on 144 degrees of freedom
Multiple R-squared: 0.009073, Adjusted R-squared: 0.002191
F-statistic: 1.318 on 1 and 144 DF, p-value: 0.2528
次に、「うどんとそばとその他」について、naiveに回帰をしてみると、「そばは有意に勝率を高めそう」という結果が得られました:
summary(lm(win ~ udon_soba_else, data=df_cp))
Residuals:
Min 1Q Median 3Q Max
-0.8837 0.0000 0.1163 0.2099 0.2099
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.79012 0.03926 20.127 <2e-16 ***
udon_soba_elsesoba 0.20988 0.08494 2.471 0.0147 *
udon_soba_elseudon 0.09360 0.06666 1.404 0.1625
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3533 on 143 degrees of freedom
Multiple R-squared: 0.04465, Adjusted R-squared: 0.03129
F-statistic: 3.341 on 2 and 143 DF, p-value: 0.03817
ただし、0 or 1のoutcomeを回帰した信頼区間が算出されているということに注意してください。
交絡を調べるために、covatiate と処置の関係性を可視化します。
まず、最高気温と冷たいご飯フラグの関係性を見てみると、気温が高いと冷たいご飯を頼みやすいことがわかります:
df_cp %>% mutate(cold_flag = as.factor(cold_flag)) %>% ggplot(aes(x=max_temperature, y = ..density..)) + geom_histogram(aes(fill=cold_flag), alpha=.2, position="identity", binwidth = 3) + geom_density(aes(colour=cold_flag), alpha=.2)
また、気温が低いとそばを頼みやすいこともわかります:
ggplot(df_cp, aes(x=max_temperature, y = ..density..)) + geom_histogram(aes(fill=udon_soba_else), alpha=.2, position="identity", binwidth = 3) + geom_density(aes(colour=udon_soba_else), alpha=.2)
さらに、相手の勝率が高いときは、「うどんやそば」以外を頼みやすい傾向が確認できます（考察: うどんやそばでは強敵に対峙するエネルギーが足りないのでしょうか？）:
ggplot(df_cp, aes(x=ct_win_ratio, y = ..density..)) +
geom_histogram(aes(fill=udon_soba_else), alpha=.2, position="identity", binwidth = .02) +
geom_density(aes(colour=udon_soba_else), alpha=.2)
次に、covariate と勝率の関係性について可視化します。
負けている試合の方が最高気温は高そうです:
df_all <- read.csv("Downloads/fujii_merge_all.csv") df_all %>% mutate(win = as.factor(win)) %>% ggplot(aes(x=max_temperature, y = ..density..)) + geom_histogram(aes(fill=win), alpha=.2, position="identity", binwidth = 3) + geom_density(aes(colour=win), alpha=.2)
また、負けている試合の方が相手の勝率は高そうです:
df_all %>% mutate(win = as.factor(win)) %>% ggplot(aes(x=ct_win_ratio, y = ..density..)) + geom_histogram(aes(fill=win), alpha=.2, position="identity", binwidth = .03) + geom_density(aes(colour=win), alpha=.2)
よって、最高気温と相手の勝率で交絡を調整する効果はありそうだということがわかります。
Rでの交絡調整
IPWを用いて交絡を調整します。Stabilized weightsを使います。
まずは冷たいご飯フラグを処置として処置のモデルを計算します:
model_ipw <- ipwpoint( exposure = cold_flag, family = "binomial", link = "logit", numerator = ~ 1, denominator = ~ max_temperature + ct_win_ratio, data = df_cp) df_cp$ipw_cold <- model_ipw$ipw.weights
weightsの分布です。いびつな分布にはなっていないのでよしとします:
ipwplot(weights = model_ipw$ipw.weights, logscale = FALSE, main = "Stabilized weights", xlim = c(0, 8))
IPWによって最高気温の分布が釣り合っていることはこちらで確認しました:
df_cp %>% mutate(cold_flag = as.factor(cold_flag), ipw = ipw_cold / sum(ipw_cold)) %>% ggplot(aes(x=max_temperature, y = ..density.., weight=ipw)) + geom_histogram(aes(fill=cold_flag), alpha=.2, position="identity", binwidth = 3) + geom_density(aes(colour=cold_flag), alpha=.2)
IPWで処置効果を推定した結果、やはり冷たいご飯フラグの係数は有意にはなりませんでした（95%信頼区間が0を含んでいます）:
set.seed(71) orig_boot <- function(N, df, treatment, weights){ ate <- numeric(N) for (i in 1:N){ ind <- sample(1:nrow(df), replace = TRUE) ate[i] <- ipw_hajek(df[ind,], treatment, weights[ind]) } c(mean(ate), sort(ate)[c(0.025*N, 0.975*N)]) } orig_boot(1500, df_cp, "cold_flag", model_ipw$ipw.weights)
# mean, 0.025%, 0.0975% [1] -0.06024249 -0.24196345 0.08885549
最後に、「うどんとそばとその他」の変数を処置として計算します。
model_ipw_udon_soba_else <- ipwpoint( exposure = udon_soba_else, family = "multinomial", link = "logit", numerator = ~ 1, denominator = ~ max_temperature + ct_win_ratio, data = df_cp) df_cp$ipw_udon_soba_else <- model_ipw_udon_soba_else$ipw.weights
先ほどと同じように結果を見ていきます。
weightの分布に大きな問題はなさそうです:
ipwplot(weights = model_ipw_udon_soba_else$ipw.weights, logscale = FALSE, main = "Stabilized weights", xlim = c(0, 8))
IPWによって最高気温の分布が釣り合っています:
df_cp %>% mutate(ipw = ipw_udon_soba_else / sum(ipw_udon_soba_else)) %>% ggplot(aes(x=max_temperature, y = ..density.., weight=ipw)) + geom_histogram(aes(fill=udon_soba_else), alpha=.2, position="identity", binwidth = 3) + geom_density(aes(colour=udon_soba_else), alpha=.2)
相手の勝率の分布も釣り合っています:
df_cp %>% mutate(ipw = ipw_udon_soba_else / sum(ipw_udon_soba_else)) %>% ggplot(aes(x=ct_win_ratio, y = ..density.., weight=ipw)) + geom_histogram(aes(fill=udon_soba_else), alpha=.2, position="identity", binwidth = .03) + geom_density(aes(colour=udon_soba_else), alpha=.2)
IPWで推定結果を見ると、やはりそばが強いです！有意！（うどんは有意でない）
orig_boot_3 <- function(N, df, weights){ ate1 <- numeric(N) ate2 <- numeric(N) for (i in 1:N){ ind <- sample(1:nrow(df), replace = TRUE) m <- lm(win~udon_soba_else, data=df_cp[ind,], weights = weights) ate1[i] <- m$coefficients[["udon_soba_elseudon"]] ate2[i] <- m$coefficients[["udon_soba_elsesoba"]] } c(mean(ate1), sort(ate1)[c(0.025*N, 0.975*N)], mean(ate2), sort(ate2)[c(0.025*N, 0.975*N)]) } orig_boot_3(1500, df_cp, model_ipw_udon_soba_else$ipw.weights)
# mean, 0.025%, 0.0975% (udon) | mean, 0.025%, 0.0975% (soba) [1] 0.09588222 -0.06569996 0.23156328 0.21019828 0.12049088 0.31626702
一方、udon_or_sobaという処置についてIPWで推定した結果、処置効果は有意ではないという結果になりました:
orig_boot(1500, df_cp, "udon_or_soba", model_ipw_udon_or_soba$ipw.weights)
# mean, 0.025%, 0.0975% [1] 0.05500092 -0.07112215 0.17096173
ついでに
うどんとそばとその他の変数が重要そうなので、もう少し調べてみます。
年度ごとの注文数を眺めてみると、1年目はうどんかそばを必ず食べていたという驚きの結果になりました。
また、今年度はまだ一度もそばを食べていません。
初めてそばを食べたときの結果が楽しみです。
df_cp %>% mutate(fiscal_year = as.factor(fiscal_year)) %>% ggplot(aes(fiscal_year)) + geom_bar(aes(fill=udon_soba_else))
勝率の遷移まで調べてみるとこのようになります。
序盤に「トレンドは考慮しない」と言いましたが、今年度はelseの勝率が高まってきているので、何かしら変化があったのかもしれません:
df_cp %>% mutate(fiscal_year = as.factor(fiscal_year)) %>% group_by(fiscal_year, udon_soba_else) %>% summarize(mean = mean(win_int), sd = sd(win_int)) %>% ggplot(aes(x=fiscal_year, y=mean, colour=udon_soba_else, group=udon_soba_else), alpha=.5) + geom_line() + geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd, width = 0.3))
トレンド情報、対局場所の情報、盤面の情報などまだまだ考慮できていないものがあるので、本分析結果は暫定のものとして捉えていただけると幸いです。
分析に誤りがありましたら、そっと教えていただけると幸いです。
謝辞
以下の方々に感謝します。皆さんの益々のご活躍をお祈り申し上げます。
おわりに
そばを食べたくなってきました。
*1:モデルの評価や残差分析は今後の課題とします
*2:興味のある処置とoutcomeの共通の原因によって発生するバイアス
*3:情報解析目的であること、特別な規約が書かれていなかったことから、負荷をかけすぎないアドホックな処理であれば問題ないと判断しました
*4:https://www.chunichi.co.jp/amp/article/92847
*5:こちらを参考にしました: https://qiita.com/saltcooky/items/451b96b3f346cb93a0b4
