M-1グランプリの結果の納得度を適切な仮説検定で調べてみる
はじめに
つい最近『P値 ―その正しい理解と適用』という本を読みました。
P値 ―その正しい理解と適用― (統計スポットライト・シリーズ)
- 作者: 柳川堯,島谷健一郎,宮岡悦良
- 出版社/メーカー: 近代科学社
- 発売日: 2018/11/28
- メディア: 単行本
- この商品を含むブログを見る
本当はベイズについて学びたいと思っているのですが,周りに理解者は少ないし,知識欲が高い人でも「学習コストが〜〜」という意見があったりして,やはり学位を取るまでは頻度主義で行くしかないかなぁと考えていたところ,ちょうどよい機会だったので手に取ってみました。
読みやすく薄い本で,短時間でP値の理解を深めることができます。ただ後半に進むほど数式が増えて,理解が難しくなったので8章以降は読んでいません。統計ユーザーにはそれくらいでいいかなと思います。
この『P値』ですが 「p値について正しく理解して統計的仮説検定を適切に行いましょう」という啓蒙書として読むことができます。p値への誤解を解く3章は「p値は小さいほどよいというのはまちがい」という基礎知識を知っている人にとっても,読めばさらにp値への理解が深まると思います。
5章では適切な仮説検定を行う手順が解説されています。簡単にまとめると
- 意味ある差と見なせる基準(効果量)を決めましょう。
- 有意水準,検出力(検定力)を設定して,サンプルサイズを決めましょう。
- データを集めて分析し,p値を使って有意かどうか2分法(有意傾向とかいう謎ワードを使わない)で判断しましょう。
です。
サンプルサイズが大きいほどp値は小さくなるので,集められるならたくさんデータがあったほうがいいというわけでないわけです。
せっかく勉強したので実践しみようと思います。
本エントリーのねた
2018年のRアドベントカレンダー(5日目)のエントリーなので,もちろんRを使います。では,今回の分析ネタについて説明します。
先日(2018年12月2日)にM-1グランプリが開催されました。優勝したのは「霜降り明星」だったのですが,いっしょに視聴していた妻が「結果に納得いかない」みたいな態度だったので(妻は和牛ファン。メジャーになる前から応援していた),世の大勢がこの結果(審査員の判定)についてどう思っているのか調べてみることにしました。
分析手法とサンプルサイズの設計
M-1グランプリの結果について「納得しているか/していないか」について調査したいので,二項検定(もしくは母比率の検定)を使います。サンプルサイズを決めるために効果量を設定します。
効果量は「意味ある違い(インパクトの大きさ)」があると見なせる評価指標です。今回はそこそこ小さくてもOKと考えて,0.2にします。有意水準と検出力については,それぞれ一般的とされる5%と0.8に設定して両側検定を行います。
サンプルサイズを出すには,pwr パッケージを使います。pwr については,こちらのページでまとまっています。
http://monge.tec.fukuoka-u.ac.jp/r_analysis/effect_size_01.html
今回は二項検定なので,サンプルサイズを決めるために pwr.p.test() を利用します。決めたいサンプルサイズ n を NULL すればOK。
> library(pwr) > pwr.p.test(h = 0.2, n = NULL, sig.level = 0.05, power = 0.8, alternative ="two.sided") proportion power calculation for binomial distribution (arcsine transformation) h = 0.2 n = 196.2215 sig.level = 0.05 power = 0.8 alternative = two.sided
結果から,サンプルサイズ(上の結果のn)は200ほど必要だとわかりました。次に,データ集めに移ります。
データ集めと概要
サンプルはM-1グランプリの結果を知っている人から無作為に集めるのがいいのでしょう。しかし,現実的に個人でやるにはムリなので,アンケート調査の「アンとケイト」を利用しました。
https://research-ssl.ann-kate.jp
質問は『今年(2018年)のM-1グランプリの結果をご存じの方に質問します。優勝が「霜降り明星」であることについて,あなたは納得していますか? していませんか? どちらかお答えください』です。回答は「納得している」「納得していない」の2択で,ランダムに表示されるようにしました。
費用は10800円でした。本年のブログで5000円ほどアマゾンアフィリエイト収益があったので,費用の半分をそれにあてたと考えています。今年もいろんな人のブログや twitter の情報にたいへん助けられた(身近に量的研究している人がいないので孤独なんです)ので,個人的歳末還元祭としました。
データはBOXからダウンロードできます。ご自由にお使いください(アクセス後、右上の「ダウンロード」からダウンロード可能です)。
Box
https://app.box.com/s/891kjk5adtt442wsjbgk2jyyvm4ovl3r
得られたデータの概要は次のとおりです。
性別 | 人数 |
---|---|
女性 | 118 |
男性 | 82 |
世代 | 人数 |
---|---|
15歳未満 | 2 |
15歳~19歳 | 11 |
20歳~29歳 | 44 |
30歳~39歳 | 41 |
40歳~49歳 | 36 |
50歳~59歳 | 32 |
60歳以上 | 34 |
納得している? | 人数 |
---|---|
納得していない | 79 |
納得している | 121 |
「納得していない」が79名,「納得している」が112名でした。納得している人が多いですね。帰無仮説検定に進みます。
分析と結果
今回は二項検定なので,binom.test() を使います。特別なパッケージは必要ないです。世間の結果への納得度に偏りがないならば,1人ごとの回答が「納得している/していない」のどちらかになる確率(二項分布のパラメータ)は0.5になります。
帰無仮説は「納得度に偏りがない(パラメータ=0.5)」,対立仮説は「納得度に偏りがある(パラメータ≠0.5)」です。データより,200人中121人が「納得している」と答えたので,分析コードは次のとおりになります。x が「納得している」と答えた人の数,n がサンプルサイズ,p がパラメータです。
> binom.test(x = 121, n = 200, p = 0.5, alternative = "two.sided", conf.level = 0.95) Exact binomial test data: 121 and 200 number of successes = 121, number of trials = 200, p-value = 0.003635 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.5336036 0.6732350 sample estimates: probability of success 0.605
分析結果より p 値が 0.05 より小さく,有意水準5%で差があると言えました(95%CI [0.53, 0.67])。
おわりに
どうやら世間一般では優勝の結果に納得している人が有意に多いようです。このことを妻に伝えると「ふ〜ん。それは何も考えていないから」と辛辣な答えが返ってきました・・・。現場からは以上です。
補足
なお,最近の類書としては『伝えるための心理統計』と『心理学のためのサンプルサイズ設計入門』があります。
- 作者: 大久保街亜,岡田謙介
- 出版社/メーカー: 勁草書房
- 発売日: 2012/01/26
- メディア: 単行本
- 購入: 9人 クリック: 164回
- この商品を含むブログ (13件) を見る
- 作者: 村井潤一郎,橋本貴充
- 出版社/メーカー: 講談社
- 発売日: 2017/03/08
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
『伝えるための〜』を以前読んだときは難しく感じたので,『P値』→『伝えるための〜』という順番がよさそうに思います。さらに実践していくときに『サンプルサイズ設計入門』を手に取るという流れがいいかもしれません(まだ積ん読状態)。より専門的なものには『サンプルサイズの決め方』がありますが,私には手が余るので紹介しません。
最後に,ブログでも許容度が過ぎたり不正確すぎる記述・まちがいがあれば教えてください。
purrr の map を使って検定を繰り返す
ずっとKJ法とかポストイットとかを使ってアンケートデータと格闘しているのでまったく統計分析はしていませんでした。少し目処が立ったので、次の研究のための予備調査を兼ねて昔のデータを見直すことにしました。
こんなデータがあるとします。
Xが性別みたいなカテゴリカル変数で、A~Dまでは7件法のデータだとします。A~Dでそれぞれt検定を行いたいのですが
t.test(hoge$A ~ hoge$X)
t.test(hoge$B ~ hoge$X)
・・・
と繰り返すのはちょっとイヤで何か簡単にt検定を繰り返すことはできないかと方法を探しました(今回のように5回くらいならガマンしますが、今後10とか20になりそうだったから)。
繰り返しだから purrr を使えばなんとかなるだろうと検索しても、なかなかお目当てのものが出てこないのでここにメモします。たぶん、検索のやり方が下手だったのだろうと思います。
実行
str(hoge) # データを確認 Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 12 obs. of 5 variables: $ X: num 1 2 2 1 1 2 1 1 2 2 ... $ A: num 4 2 1 1 4 2 3 5 3 4 ... $ B: num 7 7 3 3 1 2 3 5 3 6 ... $ C: num 7 4 6 7 3 5 6 7 6 1 ... $ D: num 1 4 7 4 5 3 2 3 2 7 ... hoge$X <- as.factor(hoge$X) # Xを因子化
ここまでで下準備OK。流れとしては
1. データからXをいったん省いて分析対象A~Dに絞り、
2. map にA~Dを流し込み、
3. t検定を繰り返し実行し、
4. 結果を kekka に保存する。
library(tidyverse) # library(purrr)もOK kekka <- hoge %>% select(-X) %>% map(~t.test(. ~hoge$X)) # t.test( の次の . が流し込まれるデータ)
無事できました。これでデータが10とか20になっても大丈夫。もっと上手な(オシャレな)方法があるのだと思います。あれば教えてください。
参考(やっぱりあった)
iPad 上で RStudio を動かす(小ネタ)
RStudio Advent Calendar 2017 の20日目です。先月の Hijiyama.R で RStudio.Cloud について少し報告したのですが,よく考えれば,ネットがつながれば使えるわけです。実際,iPhone で RStudio.Cloud を使ってみた,という報告も Twitter で見かけたことがあります(出典を探せませんでした)。だったら,iPad ならどうだろう,ということでやってみました。
Hijiyama.R についてはこちらに詳しいです。
niszet.hatenablog.com
使用したのは,iPad Pro 12.9インチ(SIMフリーで楽天モバイルのSIMを差しています)。まず,RStudio.Cloud にアクセスして,ログイン。私は Google アカウントでいつも接続しています。
以前作ったプロジェクトを開いてみると……
できた!!
まあ,当たり前ですね。
しかし,困ったことがあります。いざスクリプトを書こうとすると,キーボードの表示によって,画面半分がおおわれてしまいます。そのため,実行時にコンソールが見えず,いちいちキーボードを隠せねばなりません。こんな感じ。
物理のキーボードがあればなんとかなるんじゃないかと思い,Smart Keyboard を買ってきました。今日もアップルストア心斎橋は混雑してましたね。セッティングしたらこんな感じです。
これで通常の使い方(?)に近づいた感があります。気をつけたい点は,スクリプト実行のショートカットは Cmd + return ではなく,option + return であることです。
なんでしょう,利便性は乏しいんですが,「Rってこんなんだよ」って披露するときに使えるんじゃないでしょうか。すごくコスパ悪いです。
画像が暗いのは,宿泊先のホテルで撮ったたためです。すみません。
極私的Rの使い方5選(2017年版)
今年「このスライドはRで作っています」と言ったら「ええ!?」とドン引きされる経験を私もしたのを思い出し、「今年Rでやったこと」を書くのをやめて、自分のR環境と使い方について書きます。対象とするのは、Rを使うのが月に10回以下の方です(それ以上使っている人はきっと達人レベルに違いないので・・・)。
周りにRを使ったりする人がいないので、かなり効率の悪い方法になっていると思います。もっといい方法があれば、ぜひ教えてください。なお、本エントリーは2017年 R AdventCalender 21日目のものです。
1. %>% をショートカットキーなしで入力する
Ctrl + Shift + M、無意識に出せますか? 私はできません。でも、パイプ演算子のようによく使うものはすぐに出せるようにしておきたい。
3つのキーを組み合わせたショートカットは難しい(覚えるほど頻繁にRを使うわけではない。ほぼ1ヶ月使わないこともある)ので、スニペットツールで解決しています。
私の場合は、;;(セミコロン2つ)で %>%
が入力されるように設定しています。普段は Mac なので、ツールはランチャーアプリの「Alfred」です。パワーパックにするとスニペットの機能が使えます。
【あとで動画を載せる】
2.BetterTouchTool で Rmarkdown を使いやすく
Rmarkdownのコードチャンク(下のようなもの)を出すときのショートカットも同じように楽な設定をしています。具体的には、指3本で下に動かすとチャンクが挿入されます。
【あとで動画を載せる】
Cmd + Option + I はなんとか覚えられそうですが、このショートカットを使うときは、いちいちキーボードに目を落とせねばならず(タッチタイプがそのレベル)、それがストレスだったのです。
もちろん、これもスニペットツールでやればいいのですが、なぜか私は BetterTouchTool のジェスチャーでチャンクが出るようにしています。どうしてそうしたのは思い出せません。
3.いちいち library()するのが面倒なときにWコロン
「このパッケージのこの関数だけを1回だけ使いたい」ということがありますよね。たぶん。そういうときに::(Wコロン)です。
どのパッケージかはっきりさせるため ::
+関数があるのでしょうが、私の場合はどちらかというと不純な動機で「ちょっと1度だけ、各変数の平均と標準偏差をざっと見たい」なんてときに psych::describe(hogehoge)
とかやっています。
4.パッケージ管理は GUI で
install.packages()
はほぼ使わなくなりました。基本、RStudio でパッケージ管理しています。
「Install」ボタンを押して、インストールしたいパッケージを頭から数文字入れるとサジェスチョンを出してくれます。だから、スペルミスが減りました。あと、だいたい月に1回くらい「Update」ボタンを押します。一気にアップデートされるときは気持ちいいです。
実はここ最近、Jupyter Notebook も使おうとしていたのですが、パッケージのインストール・管理の機能がないので、使うのをやめたくらい。
5.データを渡すときは .xlsx のファイルで
なんとか周りにRを使う人を増やそうと、今年は月1のゼミ前に1時間ほど有志で統計分析の勉強会を3月から11月まで毎月行いました。経営学に近い分野の社会人学生は文系が多いので、あまり量的分析を好みません。しかし、それではいけないので、みんなでがんばりました。
はじめは例題として分析するデータをCSVやテキストファイルで渡していました。しかし、Windows の人や Mac の人が混在しているので、それぞれ対応する文字コードにしたファイルを渡すことになり、まず「文字コード???」となってつまずきます。
だったら、「もうエクセルでいいじゃない、RStudio のGUIだとエクセルファイル読み込めるし(正確にはread_excel()
を使っているけど)、文字コードの問題も起きない」ということでエクセルで配布するようになりました。
それ以来、文字コードのトラブルは起きず、インポート前に自主的にエクセルファイルの中身を見て(だいたい Office が入っていて、無意識にダブルクリックするみたい)、分析前に確認する習慣も付きそうです。
以上、ご意見があるかと思いますが、極私的なRの使い方でした。もっといいものがあれがお教えください。よろしくお願いします。
RStudio Cloud による社会人学生への統計教育の実践
Hijiyama.R the Final で LT をしてきました。LT のタイトルは「RStudio Cloud による社会人学生への統計教育の実践」です。RStudio Cloud を使えば,みんな同じ環境でRに触れることができるよ,という内容でした。
月1ペースで社会人学生で統計分析の勉強会をしているのですが,次のような理由で分析する環境を整えるのが結構大変だと思っていました。
PC環境がバラバラで環境構築であきらめることも。
アカウント名が日本語になっていることがよくある。
RStudioやパッケージを全員がちゃんとインストールするのが難しい…。
一人ひとりフォローすればなんてことはないのですが,勉強会は月1で1時間程度で時間がないのです。どうしたもんかと悩んでいたところに救世主が!
使い方はとても簡単。右上の Login から「Log in with Google」で利用できるようになります。ログイン後,しばらく待っているとスタンバイ完了。New Project でプロジェクトを作って,いつものようにRStudioを使うだけ。
動作は重いです。データは右下のペインの「Upload」からアップします。R スクリプトを書いて保存して,URLを共有すれば,すぐに他の人も使えます。違うアカウントでログインすると,右下に作成したスクリプトやデータがありますし,共有前にパッケージをインストールしていたら,そのパッケージも使えます(しかし,rstan とはダメです)。
まとめ
- R を使ってこんなふうに分析できるよという紹介用レベルには使える。
- 重い。
以上です。
補足
kazutan 先生によれば,次のとおりだそうです。
んー、これalpha版で多分test期間中なはず。
— kazutan v3.4.2 (@kazutan) 2017年11月26日
補足2
今日はいつもと違って,Jupyter Notebook + RISE でプレゼンしました。
シンプソンのパラドックスと階層ベイズ
丸の内オアゾに入っている書店をぶらついているときに『データ分析をマスターする12のレッスン』という本を見つけました。最近出た学部3〜4年生向けの本のようです。パラパラと立ち読みしたところ,商学と政策の先生が書いた本で自分にちょうどよさそうなレベルだったのでお買い上げ。
データ分析をマスターする12のレッスン (有斐閣アルマBasic)
- 作者: 畑農鋭矢,水落正明
- 出版社/メーカー: 有斐閣
- 発売日: 2017/10/12
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
その中で,おもしろい例題がありました。各都道府県の合計特殊出生率と女性労働力率の関係です。グーグルによると,それぞれ次のとおりです。
合計特殊出生率(tfr):一人の女性が一生に産む子供の平均数
女性労働力率(flp): 15歳以上の女性の人口に占める、実際に働いている、もしくは求職中の女性の割合。
横軸を女性労働力率,縦軸を合計特殊出生率でプロットすると,次のようになります。
散布図を見る限り,負の相関があります。実際に相関係数は -0.19。しかし,このデータは1980年と2000年のデータがいっしょになっていて,年代を色分けし,回帰直線を引いてみると次のようになります。青い回帰直線は年代別にしなかった場合のもの。
全体と層別にしたときに結果が異なるシンプソンのパラドックスが起きています。層別にすると,女性労働力率が大きいほど合計特殊出生率が大きくなって「女性が働くほど子どもが生まれやすい」なんて説明がつきそうです。
『データ分析をマスターする12のレッスン』ではトレンド項というのを入れたりして,うまく解決していくのですが,年代でネストされていて,都道府県の個別事情(効果)も踏まえたマルチレベル分析でもやれるのではないかと思い,練習がてら,やってみることにしました。
ベイズ推定でやる! ちょうど犬4匹本の輪読会も開催するし
各都道府県の年代別で線を引いてみると,そんなに傾きは変わらないようなので,ランダム切片だけに絞ります。
このレベルの Stan コードはまだ自前で書けないので,brms パッケージ頼み 。
prior = c(set_prior("cauchy(0, 3)", class = "sd", group = "prefID"), set_prior("cauchy(0, 3)", class = "sd", group = "yearID")) library(brms) kekka.1 = brm(tfr ~ 1 + flp + (1 | prefID) + (1 | yearID), data = d, iter = 5000, prior = prior, control = list(adapt_delta = 0.9, max_treedepth = 12))
上段はランダム切片の事前分布です。この事前分布にしたのは,stan レファレンスマニュアル(日本語版)を参考にして設定しました。もっといい設定があれば教えてください。下段がベイズ推定をしているところ。control の部分はオプションです。実行すると次のようなエラーが出ます。アドバイスに従って,最大値まで上げても消えなかったので,まあまあの落としどころとして設定しました。
警告メッセージ: There were ●●●● divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
推定結果
次のとおりで,Rhat だけで判断すると収束しています(有効サンプルサイズは気にしなくていいかな?)。
> kekka.1 Family: gaussian(identity) Formula: tfr ~ 1 + flp + (1 | prefID) + (1 | yearID) Data: d (Number of observations: 94) Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; total post-warmup samples = 10000 ICs: LOO = NA; WAIC = NA; R2 = NA Group-Level Effects: ~prefID (Number of levels: 47) Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sd(Intercept) 0.13 0.02 0.10 0.17 1232 1.00 ~yearID (Number of levels: 2) Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sd(Intercept) 1.22 1.21 0.15 4.48 953 1.01 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 1.69 1.09 -0.77 4.10 634 1.01 flp -0.00 0.00 -0.01 0.01 1095 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma 0.05 0.01 0.04 0.06 1590 1.00 Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1). 警告メッセージ: There were 402 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
固定効果(Population-Level Effects)を確認すると,傾きはほぼゼロ。切片のランダム効果を見ると,都道府県別の効果(~prefID)は0.13。年代による効果(~yearID )は1.22で,年代による影響が多いとわかります。この点推定値は事後平均です。
図示するとこんな感じです。b_Intercept は固定効果の切片で,sd_yearID_Intercept は年代による切片のランダム効果です。b_Intercept の95%信用区間が0をまたいでいますが,事後平均値も事後最大値も0にはなっていないので,大きな問題はないかと思います。
こちらの図のほうが見やすいでしょうか(個人的には上のほうが好き)。太い線が50%信用区間,細い線が90%信用区間です。点は事後中央値。
以上,階層ベイズモデルによる結果を見ると「(このデータだけを使うと)女性労働力率があがると合計特殊出生率も上がるという関係はなく,都道府県による違いも小さく,時代の変化によって合計特殊出生率は下がっている」と言えそうです。
いい練習になるので,また『データ分析をマスターする12のレッスン』をネタにするかもしれません。
P.S.1
divergent transitions after warmup はどうしてもなくなっていないので,分析結果の妥当性がバッチリあるとは言えなさそうです。このエラーは階層モデルのときに起こりやすいみたい。推定するパラメータとデータの数が釣り合わないのかな・・・?
P.S.2
Mplus を使ってマルチレベルSEM&ベイズ推定もやってみようと思いましたが,やり方がわかりませんでした・・・。
今回,参考にした本
データ分析をマスターする12のレッスン (有斐閣アルマBasic)
- 作者: 畑農鋭矢,水落正明
- 出版社/メーカー: 有斐閣
- 発売日: 2017/10/12
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版
- 作者: John K. Kruschke,前田和寛,小杉考司
- 出版社/メーカー: 共立出版
- 発売日: 2017/07/22
- メディア: 大型本
- この商品を含むブログを見る
StanとRでベイズ統計モデリング (Wonderful R)
- 作者: 松浦健太郎,石田基広
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
- この商品を含むブログ (8件) を見る
- 作者: 清水裕士
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2014/10/01
- メディア: 単行本
- この商品を含むブログを見る
- 作者: Richard McElreath
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 2016/02/19
- メディア: ハードカバー
- この商品を含むブログを見る
Proximity Matching による因果効果の推定?
2017年7月30日のKazutan.Rで「Proximity Matching による因果効果の推定」というタイトルで LT をしてきました。スライドは↓ですが,R のスクリプトは少なめなので,その補足説明をここに残しておこうと思います。
https://hikarugoto.github.io/kazutanR20170730/
LT の動機
今回の発表のモチベーションは,たまたま次の論考(何かの予稿集かな?)を見かけ,傾向スコアに代わる何かがあるのだったら,調べてみようかなと思った次第です。ランダムフォレストと傾向スコアについては,@TJO さんや@nakamichi さんなどがすでにブログで書いていらっしゃいますが,この論考では Proximity というのを使ってマッチングを行うものなので,ちょっと方法が違うようです。
Zhao, P., Su, X., Ge, T., & Fan, J. (2016). Propensity score and proximity matching using random forest. Contemporary Clinical Trials, 47, 85-92.
ランダムフォレストから Proximity を出して,マッチングを行う例は,すでに “Causal inference with observational data in R: A walk-through from some recent research” という YouTube 動画と資料があったので今回の LT ではそれをなぞってみました。
使用したデータ
統計的因果推論といえば,岩波データサイエンス vol.3 ! そちらのサポートページで公開されているデータを使いました。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2016/06/10
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (4件) を見る
10000行で35変数ありますが,「CM視聴の有無でゲーム利用回数はどの程度異なるか」をアウトカムにして,23変数を使うことにしました。イメージ的にはこんな感じ(あんまりいい図ではないけど)。
上の図は DiagrammeR で書きました。
grViz(" digraph boxes_and_circles { # a 'graph' statement graph [overlap = true, fontsize = 10] # several 'node' statements node [shape = box, fontname = Helvetica] 原因変数(CM視聴の有無); ゲーム利用回数 node [shape = circle, width = 1] // sets as circles 共変量 # several 'edge' statements 共変量->{原因変数(CM視聴の有無) ゲーム利用回数} 原因変数(CM視聴の有無)->{ゲーム利用回数} } ")
使用したパッケージ
使ったパッケージは次のとおりです。
library(ggplot2) library(randomForest) library(dplyr) library(magrittr) library(broom) library(knitr) library(DiagrammeR) library(ggpubr)
分析の手順
まず,使用する変数を元のデータから取り出します。select で -(マイナス)を使い,不要な列を削除しました。
CI_data %<>% select(-gamedummy, -area_keihan, -job_dummy8, -fam_str_dummy5, -T:-M3, -gamesecond) CI_data %>% glimpse() Observations: 10,000 Variables: 23 $ cm_dummy <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... $ area_kanto <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ... $ area_tokai <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ... $ area_keihanshin <int> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, ... $ age <dbl> 44.5, 34.5, 24.5, 44.5, 34.5, 24.5, 19.0, 54.5, 44.5, 34.5, 34.5, 34.5, ... <以下,省略>
%<>%
はあまり見かけませんが,左からデータフレームを右に流していって,その後,同じ名前のオプジェクトに置き換える破壊的な操作をするパイプです。今回は軽い分析なので使ってみました。 2行目はチラ見関数です。str()
よりもモダンな感じ?
ランダムフォレストの実行
proximity を得るために,ランダムフォレストを実行します。
cm_dummy = CI_data$cm_dummy==1 # 論理ベクトル rf = randomForest(cm_dummy ~ ., data = CI_data, proximity = TRUE) # 1分くらいかかる
1行目は CM視聴した有無を論理値ベクトルで入手し,それを randamForest
で使います。なぜ,論理値ベクトルにする必要があるのか,私には不明です。proximity = T
のオプションを忘れないのが重要です。
Proximity でマッチングさせたデータを作る
参考資料では,以下のようにして,マッチングを行っています。
prox = rf$proximity - diag(nrow(rf$proximity)) # 何やってるの? prox.true = prox[cm_dummy, !cm_dummy] # 何やってるの? untre.samples = apply(prox.true, 1, which.max) # 行ごとに proximity が最大の列を抜き出し? tre = CI_data[cm_dummy,] # 処置群のみ取り出し untre = CI_data[!cm_dummy,] # 対照群のデータのみ取り出し untre.samples2 = unique(untre.samples) # 重複している行を削除。 untre2 = untre[untre.samples2,] # 処置群と似ている対照群を抜き出している? whole.sample.u = rbind(tre, untre2) # それぞれのデータを結合。
もっとも重要なところなのでしょうが,勉強不足でよくわからん・・・。目をつぶって先に進みます。
共変量のバランスは取れたのかチェックしてみる。
せっかくなので,CM見た群と見なかった群の共変量の分布が近づいたかどうか,TVwatch_day と pmoney について見てみます。マッチング前後で見比べるため,ggplot で作図したものを並べます。ggplot を並べる方法はいくつかありますが,今回は ggpubr の ggarrange()
を使います。
a1 = CI_data %>% ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) a2 = whole.sample.u %>% ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) b1 = CI_data %>% ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5) b2 = whole.sample.u %>% ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)
TVwatch_day のマッチング前の密度分布(a1),マッチング後(a2)として,pmoney も同様にマッチング前後(b1 と b2)を ggplot()
で作ります。並べるには次のようにggpub パッケージの ggarrange()
を使います。
ggarrange(a1, a2, b1, b2, ncol = 2, nrow = 2, labels = c("Before Matching","After Matching","Before Matching","After Matching"), common.legend = TRUE)
図を2行2列で並べて,それぞれマッチング前とマッチング後のラベルをつけ,共通の凡例を表示させます。
確認しても,これで共変量が調整されているかなんとも言えないです・・・。なお,この並べ方は次のページを参考にしました。
因果効果の推定結果を確認
参考資料には「あとは回帰分析すればいいぜ」と書いてあったので,やってみます。なんか怪しい。まぁ,RCT に近づいたからやっていいという考え方だとは思うのですが・・・。
lm(gamecount ~ ., data = whole.sample.u) %>% tidy() %>% head(3) %>% kable()
lm()
で得た結果を2行目でデータフレームにして,3行目ではじめの数行だけ表示させ,4行目で下の表を出すためのマークダウンを得ました。
term | estimate | std.error | statistic | p.value |
---|---|---|---|---|
(Intercept) | 46.467720 | 7.462751 | 6.226621 | 0.0000000 |
cm_dummy | 5.080259 | 1.499469 | 3.388039 | 0.0007088 |
area_kanto | -3.428637 | 2.041777 | -1.679242 | 0.0931606 |
上の表から,CMを見た場合の平均処置効果(ATE)は 5.08 であることがわかりました。岩波DS vol.3 でIPWによってATEを出したときの結果が 5.32 ですから(96ページ),近いと言えば近い*1。
ATE | |
---|---|
岩波DS vol.3 | 5.32 |
RF で IPW @nakamichi*2 | 5.32 |
Proximity でマッチング | 5.08 |
分析してみた感想
やっていることをよく理解できず,今回の結果を信じるには至っていません。ちなみに傾向スコアマッチングは使わないほうがいいという論考が出ています*3。今回の方法は傾向スコアの弱点をクリアしているものの*4,IPW法など他のやり方があるので,わざわざマッチングをしなくてもいいように思いました。
最後に
やはり,Kazutan.R に参加してよかったです。主催者 kazutan 先生やハンズオンの MrUnadadon氏,LTされた方々に感謝。私にはデータ分析をする人が身近にいないので,情報の入手先は twitter やこういう勉強会がすべてなので(書籍の情報も twitter で知ることが多い),助かります。
*1:ちなみに,スライドで2.32と書いたのは ATT でまちがいでした。スミマセン
*2:https://rpubs.com/nakamichi_/study-iwanami-ds3
*3:King and Nielsen (2016)"Why Propensity Scores Should Not Be Used for Matching"
*4:Zhao らによれば,傾向スコアの限界は「model misspecifications, categorical variables with more than two levels, difficulties in handling missing data, and nonlinear relationships」とのことです。