傾向スコアによる統計的因果推論をやろうとして悩むのはほどほどにしたい
観察研究に基づく因果推論は「よほど共変量についての先行研究が積み上がっていないかぎり、バイアスを取り除いたと言い切るのが大変に難しい」という事情があります。
なので私はweb業界というスピーディな業界にあって、観察研究的なアプローチが必要になったときに「因果関係についてどの程度言及すべきか?どれぐらい強い主張をしてもいいのか?」についてそれなりに悩んできました。
そのときに考えたこと、読んだ記事などについて記します。
良い共変量とは?
白米の健康上のリスク*1を啓蒙していらっしゃるので有名な(?)津川友介さんのブログの統計的因果推論に関する記事が大変面白いです。
疫学者の中にはPSモデルには交絡因子だけ含まれていれば良いと考えている人達がいます。(中略)しかし、ルービンによるとPSモデルには交絡因子とアウトカムの予測因子の両者が含まれている必要があります。
これと似た点は因果推論の赤い本*2でも言及されていて、この本のp123では次の記述があります。
「割り当てに強い関連がある変数」よりも、「従属変数に強い関連がある変数」を共変量として選ぶほうが因果効果の推定が偏りが少なく、かつ推定量の分散が小さくなる
言われてみればそりゃそうか、と思える指摘です。RCTをするにしても傾向スコア分析をするにしても、「共変量の分布を揃えてバイアスを取り除くこと」を目的にしているため、アウトカムに重大な影響を与える因子が存在するなら、共変量を考慮する上で見逃されないべきなんですね。
ビジネスの問題に関する共変量の選択が難しい件について
私はCtoCのweb系企業に勤めていますが、先行研究が積み上がっているわけでもなく、因果推論にふんだんに時間を割けるわけでもなく、現実的に「強く因果関係を主張できるほどの分析」は行えていません。
とはいってもスピーディに因果関係について意見が欲しい!と頼まれる機会はあるので、「単なる回帰分析よりは、やや因果に近い主張ができる」というレベルを目指すことにしています。ドメイン知識に照らして「それってこういう交絡因子が想定できるんじゃない?」という想定ツッコミを削っていくイメージです。
考慮されていなかった共変量を利用して解析を行なっていくことこそが、人文社会学や疫学などの分野で因果関係を徐々に明確化していくことに繋がると期待される
上記は因果推論の赤い本のp129の記述です。おそらくこの考え方に近いのでは、と思います。
Magic: The Gatheringの大会記録をスクレイピングした
有名な米国産カードゲームでMagic: The Gatheringというのがあります(通称MTG)。 Youtubeでよく対戦動画を観たりなどするのですが、ふと「データ分析してみたい」と思いました。
とはいえオープンデータな棋譜が存在するわけでもないので、手に入りやすいデータ資源ってなんだろうと考えると
- 大会TOP8のデッキレシピおよびプレイヤー名
- カードのテキスト
あたりが現実的です。
本記事ではステップ1として、大会top8の戦績が公開されているMTG Decks Databaseから大会TOP8のデッキをスクレイピングしてみたときの記録を残します。
※なお、本コードを書くことよりもMTG分析の先行事例を探したりしてネタを練ることの方に時間がかかりました。
コード
R言語にてrvestというパッケージを使ってスクレイピングしました。
スクレイピング先のページ構成を調べきれていないので、まずは「エラー起きたら無視して次のページを読みに行く」というスタンスで作りました。
# まずは自作関数の定義。 # 大会の戦績ページにアクセスして、TOP8のプレイヤー名・デッキ名 etc... をまとめたmatrixを得る。 list_best_decks <- function(event.id){ tryCatch({ event.url <- paste0("https://mtgtop8.com/event?e=", event.id) event <- read_html(event.url) event.exist <- event %>% html_node(xpath = "/html/body/div[4]/div/div") %>% html_text() if (event.exist == "No event could be found."){ return(c()) } event.decks <- event %>% html_nodes(xpath = "/html/body/div[4]/div/table/tr/td[1]/div/div") event.titles <- event %>% html_node(xpath = "/html/body/div[4]/div/div/table") %>% html_children() event.about <- event %>% html_node(xpath = "/html/body/div[4]/div/table/tr/td[1]/div/table/tr/td/text()[2]") %>% html_text() decks <- c() for (deck in event.decks){ deck.detail <- html_children(deck) decks <- rbind( decks, c( event.id, event.titles[1] %>% html_text(), event.about %>% str_extract("[0-9]{2}/[0-9]{2}/[0-9]{2}"), event.about %>% str_extract("^[0-9]+ players") %>% str_extract("[0-9]+"), deck.detail[1] %>% html_text(),# rank deck.detail[2] %>% html_text(),# name deck.detail[3] %>% html_text(),# player deck.detail[2] %>% html_nodes("a") %>% html_attr("href") %>% str_extract("d=[0-9]+") %>% str_remove("d="), # deck_id deck.detail[2] %>% html_nodes("a") %>% html_attr("href") %>% str_extract("f=[0-9A-Za-z]+") %>% str_remove("f=") # format ) ) } colnames(decks) <- c("id","title", "date", "n_players", "rank", "name", "player", "deck_id", "format") }, error = function(e) { message("ERROR!") message(e) decks <- c() }, silent = TRUE) return(decks) } decks <- c() for(i in (整数):(整数)){ # スクレイピング先のサイトはeパラメータに大会のidを指定してページを分けているので、ここでその範囲を指定する。 tryCatch({ decks <- rbind( decks, list_best_decks(i) ) }, error = function(e) { message("ERROR!") message(e) }, silent = TRUE) Sys.sleep(1) }
参考文献はRの宇宙本です。
上記コードを実行した結果
以上によって、まずは2015年10月〜2018年2月の期間で、3万件程度の大会上位デッキのリストを獲得しました。
概観1:レガシーのデータが一番多い
取得したデータの範囲では、レガシーの大会のデータが最も多かったです。
# A tibble: 13 x 2 format cnt <chr> <int> 1 BL 1 2 CHL 12 3 EDH 478 4 EDHM 243 5 EX 6 6 HIGH 60 7 LE 1406 # レガシー 8 LI 191 9 MO 1302 # モダン 10 PAU 337 11 PEA 333 12 ST 1097 # スタン 13 VI 609 # ヴィンテージ
概観2:その中でも最人気はグリデル
レガシーフォーマットの中でTOP8ランクイン件数が多い順にデッキを並べるとグリクシスデルバーが一番でした。
# A tibble: 902 x 2 name cnt <chr> <int> 1 Grixis Delver 569 2 Miracles 312 3 UW Miracles 307 4 Eldrazi Aggro 288 5 Elves 275 6 Death & Taxes 265 7 Storm 250 8 Lands 200 9 Burn 175 10 Shardless BUG 173
というわけで、戦績データを集めることが出来ました。これを元に何かプレイヤーの意思決定に資する分析成果を出せるよう料理していこうと思います。
需要があるかどうか分からないので、分析成果を読んでやるぞというプレイヤーの方がもしいらっしゃれば大変ありがたい、
matplotlibを使いこなして複数のword cloudをpdf1枚にまとめる
pythonの大変便利なライブラリにwordcloudがあります。
私はgensimでトピックモデルを扱ったときに特にこれにお世話になりました。
本記事では「matplotlibを使って複数のwordcloudを単一のpdfにまとめるコード」を紹介します。
import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages from wordcloud import WordCloud from gensim.models.ldamodel import LdaModel lda = LdaModel.load(file_path) # 訓練済みトピックモデルをロードしています。 wc = WordCloud(font_path='/usr/share/fonts/truetype/takao-gothic/TakaoGothic.ttf') pdf = PdfPages(pdf_path) fig, axes = plt.subplots(nrows=10, ncols=2, figsize=(10,30), dpi=dpi) for k, ax in enumerate(axes.flat): topic_word_freq = dict(lda.show_topic(k, 100)) wc_plt = wc.generate_from_frequencies(topic_word_freq) ax.imshow(wc_plt, interpolation='bilinear') ax.axis("off") ax.set_title("トピック{} ワードクラウド".format(k)) fig.tight_layout() pdf.savefig() pdf.close() fig.clf()
以上のコードでは、おおまかに言って以下の処理をしています。
plt.subplots()
で 10行2列のaxesを作成する。10行2列 = 20フィールドは、上記の例ではトピックの個数に一致します。- axesをイテレータとして、ここのフィールド(ax)に以下の処理をする。
wc.generate_from_frequencies()
によってワードクラウドの画像を作成する。ax.imshow()
によって、ax に画像を表示する。
- レイアウトを整えてpdfに出力。
axとはmatplotlibの基礎概念で「図表」のことをさします。axesとはその複数形です。
上記コードではplt.subplots()
で「白紙の図表20個」を作成した後で、axesをイテレーションしながら1つずつ白紙を埋めていくわけです。
私がこれまでに読んできたpython入門書(みんなのpythonやオライリーの機械学習本)では、この辺の「matplotlibの思想」的な部分が見えにくくなっていましたが、下記のチュートリアルをこなしておくとmatplotlibの全体像がスッキリと見えるのでおすすめです。ちなみに下記チュートリアルはmatplotlib公式からもリンクがあります。
BigQueryでクエリパラメータを連想配列的にSELECTする
apacheのサーバーログを「とりあえずBigQueryにETLしたみた」という状態で分析を進めるにあたり、クエリパラメータを効率的に取り扱うためのSELECT文を工夫した記録です。
なお、本文ではBigQueryで用いるSELECT文のバージョンはStandardSQLです。
サンプルデータ
今回は下記をサンプルデータとします。
ステップ1:パスとパラメータ部分をSPLITする
まずは"?"をデリミタとしてSPLITします。
SQLは以下の通りです。
SPLIT関数では配列が返されるので、OFFSET()によってインデックスを指定しています。
SELECT SPLIT(url, "?")[OFFSET(0)] AS path, SPLIT(url, "?")[OFFSET(1)] AS params FROM sample
ステップ2:paramsを"key=value"を単位として配列にする
ステップ1で作成したparamsに対して、更に"&"をデリミタとしてSPLITします。
SELECT path, SPLIT(params, "&") AS params FROM step1
ステップ3:文字列"key=value"をkeyカラムとvalueカラムに分離する
UNNEST(params)
をFROM句に指定することで、配列paramsに対するSELECT文を記述しています。
そしてSELECTにおいて「"="をデリミタとしてSPLITする」という処理を行うことでkeyとvalueに分離しています。
(for文のようにFROM句でイテレータを指定して、SELECT以下で処理を定義するイメージです。)
「"="をデリミタとしてSPLITする」という処理の結果、keyとvalueを列方向にまとめておくのにSTRUCTを利用しています。
SELECT path, ARRAY( SELECT STRUCT( SPLIT(param, "=")[OFFSET(0)] AS key, SPLIT(param, "=")[OFFSET(1)] AS value) FROM UNNEST(params) AS param ) AS params FROM step2
ステップ4:連想配列的に"user"の値を取り出す
次のようにサブクエリ内でWHERE句にkeyを指定してやることで、クエリパラメータの値を柔軟に取得することができます。
SELECT path, (SELECT value FROM UNNEST(params) WHERE key = "user") AS user FROM step3
gensimで「文書を全部メモリに載せる」をやめてジェネレータでやるようにした
10GB強のテキストファイルを自然言語処理するにあたり、「とりあえずcsvファイルをpandasで読み込む」というアプローチができないため、きちんとジェネレータを使って「ファイルを1行ずつ読む」という処理を実装した記録です。
なお、「csvファイルをpandasで」が好ましくないことはgensim公式がきちんと名言していました。1
コード
次の処理を行うサンプルコードを用意しました。
- テキストファイルに対してジェネレータを生やして、
- gensimの辞書を作成し、
- 辞書を元に生テキストをBag-of-Wordsへ変換するジェネレータを生やして、
- トピックモデル訓練する。
- 最後に、生テキストをトピックモデル のベクトルへ変換するジェネレータを生やす。
def tokenize_text(some_str): # 文字列を正規化したり、形態素解析したりする処理 docs = (tokenize_text(line) for line in open('documents.txt')) dic = corpora.Dictionary(docs) class Corpus2bow(object): # イテレータとして呼ばれるたびにBag-of-Wordsのジェネレータを渡す def __init__(self, file_path): self.file_path = file_path def __iter__(self): for line in open(self.file_path): yield dic.doc2bow(tokenize_text(line)) corpus = Corpus2bow("documents.txt") lda = LdaModel(corpus=corpus, num_topics=50, id2word=dic) #トピックモデル を学習 corpus_lda = lda[corpus] # トピックモデルのベクトルを得るジェネレータ
参考情報1:まずgensim公式のお手本をおさらい
公式の introduction と tutorial を読み直しました。
特に「メモリに優しく」という側面で次の点が重要でした。
- ジェネレータでテキストファイルを読み込めば、そのままgensimの辞書やモデル訓練に使える。
- 学習したモデルを使って文書を潜在意味空間上のベクトルへと写像する部分も、ジェネレータを使って処理できる。
つまり、生テキスト→トークン化→Bag-of-Words→トピックモデルのベクトル という変換を「生テキスト1行ずつイテレーティブに」行うことができます。
参考情報2:Effective Pythonの項目17
ジェネレータが一度末尾に到達すると、以後はStopIteration例外が起きるため、同一のジェネレータを使って再度イテレーションをすることができません。
よって、本書を参考にして「イテレータとして呼ばれるたびにジェネレータを返すクラス」を作りました。ちなみにgensimのチュートリアル中でも似た実装が紹介されていました。
-
gensim公式のintroductionに
Memory independence – there is no need for the whole training corpus to reside fully in RAM at any one time
とある↩
概念整理:ベイズ信頼区間とベイズ予測区間
Stanで統計モデリングをするにあたり、「信頼区間と予測区間を用語として意識して使い分けよう」と思ったので、下記書籍p14の辺りを復習。
StanとRでベイズ統計モデリング (Wonderful R)
- 作者: 松浦健太郎,石田基広
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
- この商品を含むブログ (10件) を見る
例えば、lambda = 3のポアソン分布に従う乱数を50個抽出すると、ヒストグラムは下図のごとくになります。
これに対してパラメータの最尤推定をすると、lambda は3前後となる道理です。
一方、ベイズ推定の枠組みではlambdaの値は事後分布 に従うと考えます。
ここで想定されるlambda(より一般にはパラメータ)について、各自がお好きなパーセンタイル点において(MCMCサンプルを)区切ったものを、ベイズ信頼区間といいます。
続いて、ベイズ予測区間について。
ベイズ予測区間とは、「パラメータの確率的な生成を加味した上で、目的変数yの予測の幅」を意味します。
最尤推定では「唯一のパラメータに基づき、yの予測の幅が決まる」と考えたところを、「パラメータが生成され、それに基づきyの予測の幅が決まる」のように考えるわけです。
先ほどのポアソン分布の例で言うと、yは次の確率分布に従うということです。
では、上記の確率分布を計算するにはどうすればよいのでしょうか?
ポアソン分布例では以下のいずれかの方策を取ります。
- lambdaについて周辺化志向で考える。MCMCではサンプリング結果を「すべてのlambdaが得られた」と考え、そのlambdaによる確率分布の和を考える。
- 実際に事後分布 にしたがって lambdaを生成して、ランダムなyの生成を繰り返す。(「lambdaを生成して」の部分にはMCMCサンプルからのランダムサンプリングを利用する)
とにかく「パラメータは分布する」をキーワードとしたデータ生成でもって確率分布を考えるわけですね。
大変すっきりしました。