Knowledge As Practice

JAIST(東京)で Transformative Service Research に取り組んでる社会人大学院生の研究・勉強メモ

経営学における階層モデルの使われ方

2017年10月28日 関大梅田キャンパス Kandai MeRise で Kobe.R annex(犬4匹本輪読会)第2回を開催しました! 参加者は3名でした。対象となったのは『ベイズ統計モデリング』通称・犬4匹本の第9章と10章。

 

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

 
 
難しいです。懇切丁寧なのは確かですが,難しいです・・・。両章ともテーマは階層モデリングでした。参加者の方から9章の内容が報告されました。私のほうは10章はさておき,「経営学における階層モデルの使われ方」というテーマで発表を行いました。

 
f:id:hikaru1122:20171028155157p:plain

 
内容は『組織科学』Vol.44 No.4 「マルチレベル分析への招待」をネタにしたものです。経営学では,まずマルチレベル分析は見ないです。その中で貴重な特集号となっています。犬4匹本は階層ベイズですが,この『組織科学』の特集はベイズ推定は使われていません。マルチレベル分析がこんなふうに使われてるよ,と経営学ではない人向けへの発表です。

 
f:id:hikaru1122:20171028154830p:plain

 
なお,この特集号は全文PDFで読めます。興味がある方は↓のリンクをたどってください。

http://www.aaos.or.jp/contents/committee/?page=184

Kobe.R annex 第3回は2017年12月を予定しています。ご案内は↓で行っています。

 
kobe-r-annex.connpass.com

シンプソンのパラドックスと階層ベイズ

丸の内オアゾに入っている書店をぶらついているときに『データ分析をマスターする12のレッスン』という本を見つけました。最近出た学部3〜4年生向けの本のようです。パラパラと立ち読みしたところ,商学と政策の先生が書いた本で自分にちょうどよさそうなレベルだったのでお買い上げ。

データ分析をマスターする12のレッスン (有斐閣アルマBasic)

データ分析をマスターする12のレッスン (有斐閣アルマBasic)

 
その中で,おもしろい例題がありました。各都道府県の合計特殊出生率と女性労働力率の関係です。グーグルによると,それぞれ次のとおりです。

 

合計特殊出生率(tfr):一人の女性が一生に産む子供の平均数
性労働力率(flp): 15歳以上の女性の人口に占める、実際に働いている、もしくは求職中の女性の割合。

 
横軸を女性労働力率,縦軸を合計特殊出生率でプロットすると,次のようになります。

f:id:hikaru1122:20171023220827p:plain

 
散布図を見る限り,負の相関があります。実際に相関係数は -0.19。しかし,このデータは1980年と2000年のデータがいっしょになっていて,年代を色分けし,回帰直線を引いてみると次のようになります。青い回帰直線は年代別にしなかった場合のもの。

f:id:hikaru1122:20171023220845p:plain

 
全体と層別にしたときに結果が異なるシンプソンのパラドックスが起きています。層別にすると,女性労働力率が大きいほど合計特殊出生率が大きくなって「女性が働くほど子どもが生まれやすい」なんて説明がつきそうです。

 
『データ分析をマスターする12のレッスン』ではトレンド項というのを入れたりして,うまく解決していくのですが,年代でネストされていて,都道府県の個別事情(効果)も踏まえたマルチレベル分析でもやれるのではないかと思い,練習がてら,やってみることにしました。

 

ベイズ推定でやる! ちょうど犬4匹本の輪読会も開催するし

kobe-r-annex.connpass.com

 
都道府県の年代別で線を引いてみると,そんなに傾きは変わらないようなので,ランダム切片だけに絞ります。

f:id:hikaru1122:20171023220950p:plain

 
このレベルの 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にはなっていないので,大きな問題はないかと思います。

f:id:hikaru1122:20171023221036p:plain

 
こちらの図のほうが見やすいでしょうか(個人的には上のほうが好き)。太い線が50%信用区間,細い線が90%信用区間です。点は事後中央値。

f:id:hikaru1122:20171023222642p:plain

 
以上,階層ベイズモデルによる結果を見ると「(このデータだけを使うと)女性労働力率があがると合計特殊出生率も上がるという関係はなく,都道府県による違いも小さく,時代の変化によって合計特殊出生率は下がっている」と言えそうです。

 
いい練習になるので,また『データ分析をマスターする12のレッスン』をネタにするかもしれません。

 
P.S.1
divergent transitions after warmup はどうしてもなくなっていないので,分析結果の妥当性がバッチリあるとは言えなさそうです。このエラーは階層モデルのときに起こりやすいみたい。推定するパラメータとデータの数が釣り合わないのかな・・・?

 
P.S.2
Mplus を使ってマルチレベルSEMベイズ推定もやってみようと思いましたが,やり方がわかりませんでした・・・。

 

今回,参考にした本

データ分析をマスターする12のレッスン (有斐閣アルマBasic)

データ分析をマスターする12のレッスン (有斐閣アルマBasic)

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

個人と集団のマルチレベル分析

個人と集団のマルチレベル分析

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

第1回 Kobe.R annex(犬4匹本 輪読会)開催しました!

2017年9月24日 大阪府立大学・なんばサテライトにて Kobe.R annex 第1回を開催しました。第1部『ベイズ統計モデリング―R,JAGS,Stanによるチュートリアル― 原著第2版』の輪読,第2部 分析ネタなんでもという構成でした。参加者は5名(+会場をご提供いただいた大阪府立大学の先生)。

f:id:hikaru1122:20170924191307j:plain

 
予定は第6章〜9章の予定だったのですが,8章までとしました。もちろん全員,犬4匹本を購入済みです。以下,輪読会の簡単メモです。

 

第7章

  • JAGS の情報が載っていてわかりやすいウェブサイトがある。

http://ecological-stats.com/tag/jags/

  • MCMCスライド集

hikaru1122.hatenadiary.jp

  • stan のモデルをコンパイルするとき次のエラーが出たら,Xcode を1度起動して使用許諾に同意する(Xcode のアップデートをした直後に出る恐れ)
追加情報:  警告メッセージ: 
 命令 ''/Library/Frameworks/R.framework/Resources/bin/R' CMD config CXX 2>/dev/null' の実行は状態 69 を持ちました
追加情報:  警告メッセージ: 
 命令 '/Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB file8dae4ce1404e.cpp 2> file8dae4ce1404e.cpp.err.txt' の実行は状態 1 を持ちました  

 

疑問点

  • Gelman-Rubin統計量(184頁)はRハットのことだろうか?

→どうやらそのようです。

  • 間引き(thinning)はすべき? すべきでない?

→犬4匹本は基本thinningしない的な立場(191頁)のようだが,thinning したほうがよいという文献もあるらしい。

  • thinning をした場合,サンプリング時間はどう変わるだろうか?

→あるモデルで「iteration 1000,thin 1」「iteration 1000,thin 2」「iteration 2000,thin 2」を stan で実行してみると,90秒・90秒・180秒となった。保存されたオブジェクトのデータは8MB, 4MB, 8MB。thinning の設定でサンプリング時間は変わらないが,保存されるオブジェクトの容量は変わる。

 

第8章

  • runjagsのほうが早くて便利!?

  • 統計学エンドユーザーのための JAGS 入門

hikaru1122.hatenadiary.jp

  • 参加者の方が8章の便利関数をパッケージ化されました!

 

LT枠(分析ネタなんでも)

  • 状態空間モデルについて

  • 共変量バランシング傾向スコアを使った統計的因果推論について

 
以上です。次回もいずれ開催したいと思います。

大きい iPad Pro で論文を快適に読む方法

自分の中でだいたい固めることができたので、大きい iPad Pro で快適に論文を読む方法を書いておきます。最近、質的研究ばかりやっていて、R とかあんまり触っていないため、こういうネタしかないんです。

 

結論から

 
まず先に結論を言うと、2017年9月20日現在、次の組み合わせが私にはピッタリです。

 
1. iPad Pro 12.9インチ
2. アップルペンシル
3. アップルペンシルのキャップ紛失防止ツール
4. 適度に滑りにくくするための保護シート
5. PDF Expert
6. スリーブケース

 

私の電子リーダー端末遍歴

 
タブレットスマホはいろいろ触ってきましたが、論文を読むに関しては、Retina ディスプレイの iPad(第3世代)→ソニーのデジタルペーパー(旧モデル・DPT-S1)→ 小さい iPad Pro(9.7インチ・2015年モデル)→大きい iPad Pro(12.9インチ・2017年モデル)とたどってきました。他にもKindleとかソニーのリーダーとかでもチャレンジしましたが、目には優しいものの動作が重くダメでした。

 

なぜ紙で論文を読まないのか?

 
私が紙で論文を読みたくないと思ったのは単純に「かさばるから」です。仕事柄出張が多く(昨年は224泊)、論文をたくさん持ち歩きたくありません。社会科学の論文は(少なくとも私の領域では)だいたい1本20ページくらいで、5本でA4・100ページ、10本で200ページ。出張で持ち歩くには思ったよりかさばるんです…。そのため電子化しかないと数年前からいろいろ工夫していました。

 
中には論文は「紙で読む派」という方もいらっしゃると思いますが、差し迫った問題がなければ、私はどちらでもいいと思います。お好きなほうを選べばいいです。

 

デジタルペーパー VS 小さいiPad Pro VS 大きい iPad Pro

 
電子化された論文を読む勝ちパターンがこの1年で固まりました。まず昨年、ソニーデジタルペーパー(DPT-S1)を大阪のソニーストアで実機を確認。さすがに定価で買うには抵抗があり、中古をヤフオクでゲット。それでもそれなりのお値段がしました。

 

デジタルペーパーのいいところ

○ 軽い、薄い。
○ なんか目に優しい感じがする。
○ 基本、読み書きに特化しているので、論文読みに集中できる。

 

デジタルペーパーのイマイチだったところ

△ 動作がもったり(新モデルは改善しているそうです)
△ 文字がギザギザで読みにくい(これも新モデルでは改善しているそうです)
△ 思ったよりペンが引っかかる。もっと滑らかなのが個人的には好き。
△ 英語論文を読んでいて、知らない単語が出てきたとき、別の端末で調べる必要がある。
△ 論文を読んでいて調べ物をしたい場合、ネットで快適に調べることができない。

 
数ヶ月ほど使って「ナンカチガウ」感が出てきて、だんだんと使わなくなっていきました。「これではあんまり論文が読めないなぁ」と困り、iPad Pro はどうだろうと思い立ちました。iPad Pro には大きいサイズ(12.9インチ)と小さいサイズ(当時9.7インチだった)があります。

 
大きい iPad Pro を仕事に使っている知人がいたので見せてもらって、ちょっと大きいかなと思ったので、小さい iPad Pro を使うことに決め、アップルの整備済製品で購入しました。お得です。

iPad整備済製品 - Apple(日本)

 

小さい iPad Pro のいいところ

 
○ 字がきれい。
○ 知らない英単語が出てきたとき、付属の英語辞書ですぐわかる。
○ wi-fi につなげて、快適に調べ物ができる。
○ アップルペンシルはなかなかいい感じ。

 
iPad の辞書機能は特に便利です。私は英語に不自由しているので、よく知らない単語が出てきます。これを他の電子辞書を使わずに、iPad 内蔵の辞書で調べることができるのでストレスがとても減り、QOLアップ。  

www.youtube.com

 
PDFを読むとき、PDF Expert を使うのですが、上のように知らない単語を選択して「定義」をタップすれば、内蔵辞書ですぐに意味を確認できます。超便利!

 
アップルペンシルの書き心地は動画では伝わりにくいですが、こんな感じ。

youtu.be

 

小さい iPad Pro のイマイチなところ

 
△ 使っていると、だんだん小さく感じてくる。
△ ケースに入れて使っていると、思ったより重く感じる。
△ アップルペンシルはいいけど、ちょっと滑りすぎる。

 
また半年ほど使って、また「ナンカチガウ」感が出てきたので、大きい iPad Pro を本格検討。新型の大きい iPad が出て、かなり迷いました。清水の舞台から飛び降りる気持ちで移行を決意。在庫があるなら全国どのアップルストアにも行く(出張ついでに)覚悟でした。

 
結果、銀座にアップルストアにあったので、すぐにゲット。そして現在に至ります。まだそんなに経っていないですけど、とても満足して使っています。

 

大きい iPad Pro のいいところ

 
○ 小さい iPad Pro のいいところは全部ある。
○ さらに画面が大きく、見やすい。

 

大きい iPad Pro のイマイチなところ

 
△ デカい。

 

大きい iPad Pro をより快適に使う方法

 
以下、より楽しく、快適に iPad Pro 12.9 インチを使って論文を読むために試した方法を書きます。

 

アップルペンシルのキャップカバー(必須)

 
まずアップルペンシルは必須です。しかし、アップルペンシルの最大の弱点はキャップ。小さい iPad Pro を使っているときにアップルペンシルはすでに活用していましたが、すぐにキャップをなくしました。

 
キャップがなくても普通に使えますが、なんかカッコ悪い。仕方がないので、大きい iPad Pro 導入時にキャップだけ再度購入しました。アップルストアで修理扱いです。修理の予約が必要だし、修理代は1000円くらい。

 
また紛失したら悲しいので、キャップカバーを導入。修理代より安いので、導入したほうがいいです。どのキャップカバーでも構いません。

www.youtube.com

 
ちょっと不格好なのですが、紛失のリスクを考えて、ここはガマンですかね。

 

ペーパーライク保護フィルム(必須)

 
ソニーデジタルペーパーと比べ、iPad Pro はツルッツルです。そのためアップルペンシルが滑りすぎます。多少の引っかかりが必要なため、ディスプレイが紙のような感じになる「ペーパーライク保護フィルム」を導入しました。

 
最近は安いのが出ているので、好きなのを選べばいいかと思います。ただ、気をつけてほしいことがあります。それは、せっかくの美麗なディスプレイが台無しになることです。

 
画面はちょっとざらついた見た目になりますし、指紋はむっちゃ付着します。きれいな画像・映像を見たい人にはガマンできないと思います。

 
しかし、この iPad Pro は論文を読むためにあるのです! 割り切ってエンタメは諦めましょう…。次第に目は慣れていきます。なんか違和感のある手触りも慣れます。ざらつきはありますが、画面はヌルヌル動くので、操作は問題ありません。

www.youtube.com

 

PDFリーダー(必須。お好きなのを)

 
先ほど書いたように、PDF Expert を使っています。他のPDFリーダーも試して、いちばんしっくりきたためです。かなり感覚的でちゃんとした理由はありません。ドロップボックスに読むべき論文のフォルダを作り、そのフォルダと PDF Expert を同期させます。

 
f:id:hikaru1122:20170920015459p:plain

 

格安SIM(お好きにどうぞ)

 
wi-fi でつながっていないときに(電車の中とか)で論文を読んでいて、ネットで調べ物をしたいときに便利です。ただ、つい余計なものを見てしまうときもあるので、導入するかどうかは好き好きで。私は楽天ポイント(ホテル手配は楽天トラベルだからけっこう貯まる)で支払える楽天SIMにしています。

 

www.youtube.com

 
アンテナが2本しか立っていなかったのでスピードが少し遅いですね。大きな問題はないです。

 

スリーブケース(お好きにどうぞ)

 
ケースは重いので、スリーブケースにしました。デスクにどかんと置いて使うか、手に持って使うだけなので、別に立てて使うことはないからです。アマゾンとかで好きなモノを選んでください。

 
私は純正品を使っていますが、キャップ付きのアップルペンシルだと収まりが悪く、ちょっとカッコわるいです。

f:id:hikaru1122:20170920012251j:plain

 
というわけで、大きい iPad Pro であなたも快適な論文読みライフを送りませんか? コスパには目をつむってになりますけど…。

第1回 Kobe.R 別館を開催します!

f:id:hikaru1122:20170824004200j:plain

来る2017年9月24日に Kobe.R の弟分的な存在の勉強会「Kobe.R annex(別館)」第1回を開催します。小さく、緩く(いい意味で)、楽しくやっていきたいと思います。2部制にして、第1部は『ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版』(犬4匹本)を輪読します。第2部では、初心者セッションおよびLT大会です。第1部から参加、もしくは第2部のみ参加をお選びいただけます。

 
参加者数は最大11名と小さいため、LT枠などが埋まらない場合は、参加者同士で時間の許す限りR関連の話題やRを使った研究ネタについて(アイデアレベルでも)と話し合おうと思っています。初心者セッションは “R for Datascience” または “Quantitative Social Science: An Introduction” を題材にしてもいいかな~とも思っています。

 
日曜の大阪・難波での開催です。すでに滋賀県から参加される方もいらっしゃいます。関西在住の方(もちろんその他の地域の方も)は奮ってご参加ください!

 
kobe-r-annex.connpass.com

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 ! そちらのサポートページで公開されているデータを使いました。

github.com

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

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

 
10000行で35変数ありますが,「CM視聴の有無でゲーム利用回数はどの程度異なるか」をアウトカムにして,23変数を使うことにしました。イメージ的にはこんな感じ(あんまりいい図ではないけど)。

f:id:hikaru1122:20170730181643p:plain  
上の図は 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列で並べて,それぞれマッチング前とマッチング後のラベルをつけ,共通の凡例を表示させます。

 
f:id:hikaru1122:20170730162121p:plain

 
確認しても,これで共変量が調整されているかなんとも言えないです・・・。なお,この並べ方は次のページを参考にしました。

http://www.sthda.com/english/wiki/ggplot2-easy-way-to-mix-multiple-graphs-on-the-same-page#.V5F9h53_d2E.twitter

 

因果効果の推定結果を確認

 
参考資料には「あとは回帰分析すればいいぜ」と書いてあったので,やってみます。なんか怪しい。まぁ,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。今回の方法は傾向スコアの弱点をクリアしているものの*4IPW法など他のやり方があるので,わざわざマッチングをしなくてもいいように思いました。

 

最後に

 
やはり,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」とのことです。

テキストマイニングの環境整備つまずきメモ

テキストマイニングをしようと決意して,その準備をしたところ,いくつかつまづいたところがあったので,そのメモを残しておきます。解決策はネット上に散らばっていて,同じようなトラブルにまたあったときのためです。なお,Mac です。

 

1. 環境構築の基礎知識をゲット

 
まず『Rによるやさしいテキストマイニング』(やさテキ)を読み,次に『Rによるテキストマイニング入門(第2版)』を読みました。2冊でざっとテキストマイニングの基礎知識と分析の流れをつかみます。どちらも読みやすいですが,まず『やさテキ』から『Rによる〜』の順がよいと思います。

 

Rによるやさしいテキストマイニング

Rによるやさしいテキストマイニング

Rによるテキストマイニング入門

Rによるテキストマイニング入門

 

2.MeCab のインストール時のつまずき

 
『Rによる〜』に書かれている方法で MeCab のインストールを行ったところ,うまく行きませんでした。ターミナル上で「mecab すもももももももものうち」と打っても結果が

 

??,????,,,,,?,???????,???????

 
と文字化けするのです。これは MeCab をもう1度次の方法でインストールして解決できました。

 

$ ./configure --enable-utf8-only --with-charset=utf8 --with-mecab-config=/usr/local/bin/mecab-config
$ make
$ make install

 
詳細はこちらにあります。

kumagonjp2.blog.fc2.com

 

3.RMeCab インストール時のつまずき

 
私の環境では『Rによる〜』に記載されている方法でインストールして RMecab を使うと RStudio ごと落ちるという現象が起きました(Rコンソールでもダメ)。tweet したところ,石田先生に助けていただきました。

 
これで無事解決。ありがとうございます。

 

4.mecab-ipadic-NEologd を使用するときのつまずき

 
テキストマイニングの対象はちょっと特別な言葉が入っているので,標準の辞書ではない NEologd を使うことにしました。この辞書のインストールは次のページに従って行いました。

qiita.com

 
しかし,RMeCab で NEologd を使うときは,ユーザ辞書が必要とのこと。次のページの「R」セクションの記述に従って,ユーザ辞書を作り直しました。

github.com

 
作った辞書(〜〜.csv.dic という名前のファイル)は,適当なところにコピーして使います。こんな感じ。

> RMeCabC("オールセラミック", 
      dic = "/Users/jibun/mecab-user-dict-seed.20170630.csv.dic")
[[1]]
              名詞 
"オールセラミック" 

 
標準の辞書だと,

> RMeCabC("オールセラミック")
[[1]]
  接頭詞 
"オール" 

[[2]]
        名詞 
"セラミック"

 
ちゃんと,オールセラミックを1つの言葉として扱ってくれてる。よかった。

 
まだまだテキストマイニングは準備すら難しい印象を持ちました(私の環境では)。でも『やさテキ』と『Rによるテキストマイニング第2版』とネット上の先達のおかげでなんとかなりました。KH Coder に流れたくなったけど,せっかくなんとか準備が整ったので R でがんばりたいと思います。

Rによるやさしいテキストマイニング

Rによるやさしいテキストマイニング

Rによるテキストマイニング入門

Rによるテキストマイニング入門

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 - 継承 4.0 国際 ライセンスの下に提供されています。