2017-05-20

R による整列ランク変換を使った分散分析

整列ランク変換（Aligned Rank Transform：ART）という言葉を聞いたことありますか？　私はなかったです。ART（えーあーるーてぃー。「あーと」ではない）について、Tokyo.R 第61回で LT をしてきました。補足も兼ねて ART についてまとめます。

きっかけ

　
先日、自分の専門分野の先行研究を調べていて論文（Heidenreich, et al. 2015）を読んでいたら次のような記述に出会いました。

“As our data were not normally distributed (D(243)=2.193, p<0.01), we employed the aligned rank transform (ART) procedure for hypothesis testing. For non-normal distributed data, the ART procedure is more robust and powerful than the traditional analysis of variance (ANOVA)”

　
超訳すると「データが正規分布に従わなかったので、aligned rank transform (ART) を使って検定やった。非正規分布のデータなら、ARTを使ったほうが普通の分散分析よりロバストでパワフルだ」。

　
ART について、ネット上で検索してもほとんど日本語の情報がありません。かろうじて「整列ランク変換」という訳語が見つかるくらいです。そこで少し調べてみました。

ART法はタイプIエラーをコントロールできるエクセレントな方法だ。研究者がやりがちな古典的分散分析と比べてパワーがある（Mansouri, et al. 2004）

ART はどんな分布のデータにも使えるわけじゃないけど、ロバストだし一般的なノンパラメトリック検定より望ましい性質を持っている（Higgins, et al. 1990）

ART 検定はパワフルでロバストで使いやすい（Mansouri 1998）

交互作用を調べるとき、いくつかの仮定が成り立たない場合、パラメトリックな分析は必ずしもよいわけではない。代替の策として、調整ランク変換検定（adjusted rank transform test）がある。それはパラメトリックな手法とノンパラメトリックな手法の中間みたいなものだ（Leys & Schumann 2010）

　
つまり、分散分析したい、特に交互作用に注目したい時にデータが正規分布に従わない場合、従来の分散分析は（厳密には）使えないので、整列ランク変換を施してから分散分析すればいいよ、ということだと思います。

　
R パッケージがないか探したところ、ARTool というものがちゃんとありました。CRAN に登録されています。使い方もカンタン。以下、ART を使った分散分析の方法です。例題として、性別と世代の不倫の許容度に対する違い・交互作用を見たいとします。応答変数は「結婚していても場合によっては不倫を許せる？（１：許せない～４：許せる）」です。4件法でスミマセン。

０．データの分布の確認

ggplot(rensyu, aes(furin)) + 
geom_histogram(fill = "red", alpha = 0.7, binwidth = 0.5)

f:id:hikaru1122:20170520195014p:plain

　
あきらかに正規分布には従わないっぽい。もちろん、正規性の検定はすべきですが、ここでは飛ばします。

１．ARTool を使えるようにする。

install.packages("ARTool") # 1回のみでOK。
library(ARTool)

２．データを整列ランク変換する。

henkan = art(furin ~ nendai*sex, data = rensyu)
summary(henkan)

　
art() で変換しますが、式はよくあるタイプです。式は「art(furin ~ nendai + sex + nendai:sex, data = rensyu)」を短縮した書き方にしています。どちらでもOK。lm() のときと同じです。

　
変換後、summary(henkan) の出力結果に注意します。すべての数値が０になっていることが望ましいです。今回は残念ながらそうならなかったので、警告が出ています。しかし、練習なので先に進めます。

## Warning in summary.art(henkan): F values of ANOVAs on aligned responses not
## of interest are not all ~0. ART may not be appropriate.

## Aligned Rank Transform of Factorial Model
## 
## Call:
## art(formula = furin ~ nendai * sex, data = rensyu)
## 
## Column sums of aligned responses (should all be ~0):
##     nendai        sex nendai:sex 
##          0          0          0 
## 
## F values of ANOVAs on aligned responses not of interest (should all be ~0):
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0000000 0.0000000 0.0144000 0.0007531 0.0854000

３．分散分析を実行する。

anova(henkan)

## Analysis of Variance of Aligned Rank Transformed Data
## 
## Table Type: Anova Table (Type III tests) 
## Model: No Repeated Measures (lm)
## Response: art(furin)
## 
##              Df Df.res  F value     Pr(>F)       Sum Sq Sum Sq.res
## 1 nendai      4   1261   3.7856  0.0045704  **  1935456  161178877
## 2 sex         1   1261 100.7097 < 2.22e-16 *** 10944212  137033949
## 3 nendai:sex  4   1261   4.6317  0.0010275  **  2369093  161247964
## ---
## Signif. codes:   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　
「anova(変換したデータ)」でOK。実際には、anova.art() が動いて、デフォルトでタイプ３の分散分析をやってくれます。世代と性別の主効果、交互作用は有意になりました（なお、効果量 $\eta^{2}$ は世代がほぼ0、性別が0.02、交互作用が0.01でほとんどインパクトはありません）。いちおう、交互作用のグラフをプロットします。

library(ggplot2)
ggplot(rensyu) +
aes(x = nendai, color = sex, group = sex, y = furin) +
stat_summary(fun.y = mean, geom = "point", size = 2) +
stat_summary(fun.y = mean, geom = "line", size = 1) +
xlab("世代") + ylab("不倫の許容度") +
scale_colour_discrete(name ="性別")

f:id:hikaru1122:20170520224157p:plain

伝統的な分散分析と結果を比べてみる。

aaa = lm(furin ~ nendai*sex, data = rensyu)
library(car)
Anova(aaa, type="III")

## Anova Table (Type III tests)
## 
## Response: furin
##             Sum Sq   Df  F value  Pr(>F)    
## (Intercept) 224.03    1 535.4585 < 2e-16 ***
## nendai        2.39    4   1.4269 0.22278    
## sex           2.22    1   5.3061 0.02141 *  
## nendai:sex    3.02    4   1.8051 0.12541    
## Residuals   527.59 1261                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　
同じくタイプ３でやってみましたが、こちらでは世代の主効果と交互作用が有意になりませんでした。いい感じに対比できましたが、データが怪しいので、結果はあくまで参考用として受け取ってください。

まとめ

整列ランク変換を利用した分散分析では交互作用が有意になった。
ただし、今回のデータはちゃんと整列ランク変換できてはいないのでなんとも言えない。
今後、データが正規分布に従わない場合、使ったほうがよい手法かもしれない。

さいごに

　
多重比較もできます。繰り返しのデータも分析可能のようです。詳しくはパッケージ作者の Github ページを参考にしてください。今後、必要あれば、使っていきたい分析法です。

参考文献

ARTool パッケージ
https://github.com/mjskay/ARTool
Feys, J. (2016). Nonparametric Tests for the Interaction in Two-way Factorial Designs Using R. The R Journal, 8(1), 367-378.
Heidenreich, S., Wittkowski, K., Handrich, M., & Falk, T. (2015). The dark side of customer co-creation: exploring the consequences of failed co-created services. Journal of the Academy of Marketing Science, 43(3), 279-296.
Higgins, J. J., Blair, R. C., & Tashtoush, S. (1990). The aligned rank transform procedure.
Leys, C., & Schumann, S. (2010). A nonparametric method to analyze interactions: The adjusted rank transform test. Journal of Experimental Social Psychology, 46(4), 684-688.
Mansouri, H. (1998). Multifactor analysis of variance based on the aligned rank transform technique. Computational statistics & data analysis, 29(2), 177-189.
Mansouri, H., Paige, R. L., & Surles, J. G. (2004). Aligned rank transform techniques for analysis of variance and multiple comparisons. Communications in Statistics-Theory and Methods, 33(9), 2217-2232.
Wobbrock, J.O., Findlater, L., Gergle, D. and Higgins, J.J. (2011). The Aligned Rank Transform for nonparametric factorial analyses using only ANOVA procedures. Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ‘11). Vancouver, British Columbia (May 7-12, 2011). New York: ACM Press, 143-146.

目を通してないけど、重要そうな文献

Salter, K. C., & Fawcett, R. F. (1993). The ART test of interaction: a robust and powerful rank test of interaction in factorial models. Communications in Statistics-Simulation and Computation, 22(1), 137-153.
Sawilowsky, S. S. (1990). Nonparametric tests of interaction in experimental design. Review of Educational Research, 60(1), 91-126.

2017-03-30

ggdendro と PAC 分析を使ったサービス学会の発表

ノンプログラマーのためのR入門研究に役立つツール統計分析参加報告価値共創サービス研究

2017年3月28日（火）、広島県情報プラザというところで開催されたサービス学会第５回国内大会に参加してきました。仕事以外で広島に滞在したのは小学校の修学旅行ぶりです。

http://ja.serviceology.org/events/domestic2017.html

　
今回は研究をスタートし始めたばかりのネタで、調査・分析が進んでおらず、いろいろと参加者と相談したいという理由で、ポスター発表にしました。

f:id:hikaru1122:20170330135756j:plain

　
PAC 分析というものを用いたのですが、この分析法ではデンドログラムを使います。せっかくなので、ggdendro を使ってみようと直前までポスターをいじっていました。PAC分析の手法についてはこちらを参照してください（ネット上にもたくさん情報はあります）。

PAC分析実施法入門―「個」を科学する新技法への招待

作者: 内藤哲雄
出版社/メーカー: ナカニシヤ出版
発売日: 2003/04
メディア: 単行本
クリック: 10回
この商品を含むブログ (3件) を見る

　
R で普通にクラスター分析するとこんな感じ。 f:id:hikaru1122:20170330140907j:plain

　
ポスターを作るときはこれをパワーポイントに貼り付けて、回転させてトリミングして整えていきますが、ちょっと色合いがさみしい感じ。私のポスターは茶系で地味なので、ちょっと色がほしかったのです（ggplot のあのグレーが好き）。

　
というわけで、ggdendro を使うことにしました。何もしないと ggdendro は背景色が白です。

f:id:hikaru1122:20170330141429j:plain

　
ここから ggplot2 に持って行き、縦軸と横軸を消していきます。例によって、kazutan 先生の資料を参考にしました。

https://rpubs.com/kazutan/ggdendro_test

　
できあがりはこちら。ggplot っぽくなりました。 f:id:hikaru1122:20170330141625j:plain

　
1行目で ggplot に渡すデータを作り、5行目でデンドログラムを横倒しにして、6～8行目で縦軸と横軸を削除しています。あとはこの図をパワーポイントに持っていき、情報を追加したら完成。

　
f:id:hikaru1122:20170330142947j:plain

いつも直前にホテルでかんづめして作ってます。

　
ちなみに、A０サイズで7000円ほどかかりました。高いけど、直前まで出力を受け付けてくれたキンコーズ広島本通店さん、ありがとう。

2016-12-23

旅人と RStudio

ノンプログラマーのためのR入門統計分析研究に役立つツール

このエントリーは RStudio Advent Calendar 2016 の23日目のエントリーです。ちなみ，本日は Osaka.stan 第2回に参加してきます。

　
さて，私は仕事柄，出張が多く，2016年は12月末までに224泊する予定です。だいたい1日ごとに移動しています。たぶん，日本でいちばんビジネスホテルで RStudio を使っている人間だと思います（なんの自慢にもならないですけどね！）。

　
せっかく出張が多いので，全国の名所でRStudioを使っている様子の写真を撮ってアップしようと企画したのですが，なんだか恥ずかしいのでやめました。

　
というわけで，本エントリーでは R や RStudio を付き合い始めてだいたい2年。日頃，とても勉強させていただいているRコミュニティーへの感謝も込めて，自分なりの学びや使い方を振り返ってみたいと思います。対象は例によって，統計分析やRを独学で学ばないといけなくて，ある程度時間よりもお金を少しかけられる人（つまり社会人）向けです。

パイプ演算子の出し方

RStudio のショートカットで %>% を出すには Mac の場合，シフト＋コマンド＋M です。でも，どーしてもこのショートカットが苦手なんです。いちいちキーボードを見ないと打てないし，ホームポジションから離れてしまう……。

　
で，現在は，Alfred Powerpac のスニペット機能を使っています。textexpander でもいいのでしょうが，月額課金よりは買い切り型でいろいろできる Alfred Powerpac を選びました。35英ポンド。

　
↓こんな感じでセミコロン2つで %>% が出るように設定しています。 f:id:hikaru1122:20161223023707p:plain

f:id:hikaru1122:20161223023730g:plain

　
他にも，よく読み込むパッケージを登録していたりします*1。

R と RStudio に慣れるには

身近に詳しい人がいればいいのでしょうが，私の場合はいません。コツコツ学べる才能があればいいのでしょうが，私にはそんなにありません。やらざるを得ない環境というか，追い詰められないとできないタイプです。ある程度コストをかけると「やらねば！」という気持ちになるので，そうしています。私がけっこう勉強になったなぁと思うものを2つ挙げます。

　
1．DataCamp

無料でいくつか受講できます。例えば「Introduction to R」。最初に動画で簡単な説明があって，RStudio っぽい画面で操作をしていきます。英語ですが，だいたい雰囲気はつかめます。私は昨年，3ヶ月ほど課金していくつか受講しました。

www.datacamp.com

　
1ヶ月30ドルくらい。得るものと比べたら安いくらいです。3ヶ月やれば，だいたい「これできるようになりたいな〜」が受講しきれます。難しいのはやらないし。でも最近，またちゃんとやり直そうと思っています。

　
2．『RStudioではじめるRプログラミング入門』

簡単な例からスタートして，少しずつ難しくなっていきます。私が Rの本で最後まで演習をしながら読み通せたのはこの本でした。「プログラミングとかよくわかんない」という私のようなものにはちょうどよく，勉強になりました。ただ，まったく未経験者にはいきなりこの本は難しいと思うので『はじめてのR』をざっと目を通してからのほうが無難です。

RStudioではじめるRプログラミング入門

作者: Garrett Grolemund,大橋真也,長尾高弘
出版社/メーカー: オライリージャパン
発売日: 2015/03/25
メディア: 大型本
この商品を含むブログを見る

　
なお，いちいち本なんか読んでられないという方には近い内容のものがオンライン講座になって発売されています。こちらです↓

shop.oreilly.com

　
およそ160ドルですが，私の場合，たまたま何かのセールのときに99ドルでした。

　
英語はちょっと厳しいという場合は，Udemy がいいです。これもセールだったので登録しました（見てない）。いまもセールで1800円のようです。お得ですね。

www.udemy.com

　
上記のように，課題的なものがあるとか，ステップバイステップのものが私にはあっているようでした。

もし RStudio に飽きてしまったら…

最近はいろいろと選択肢があります。R のアドカレ15日目のように Atom内で動かせる Hydrogen もありますし，Jupyter notebook とか，nteract とか，Sublimetext 内で動かしたりとか。

hikaru1122.hatenadiary.jp

　
最近はもっぱら Beaker Notebook にしています。Jupyter と違って，初めから等幅フォントで分析結果が見やすかったし，なんだか雰囲気がいい。でも，ちょっと起動が遅い。今日，Osaka.stan で Beaker Notebook を使っているグレーのスーツを着ている人がいたら，それはたぶん私です。

beakernotebook.com

それでも RStudio を手放せない理由

実は個人的に RStudio でいちばん気に入っている機能はパッケージ管理です。新しく使うパッケージのインストールも更新も楽チン。しばらく放置していて，たくさん更新があると，一気にアップデートができるので気持ちがいい。かゆい所に手が届くというか，そういうのがRStudio のいいところだと私は思います。

f:id:hikaru1122:20161223024053p:plain

f:id:hikaru1122:20161223024110p:plain

　
以上，RStudio Advent Calendar 23日目でした。

*1:.Rprofile に書いておけ，という意見もあるでしょうが，毎回これらのパッケージを使うわけではないので，必要な時に読み込んでいます。

2016-12-04

rstanで分析をするのに必要なR力（りょく）

Stan と R でベイズ統計モデリングベイズデータ分析統計分析ノンプログラマーのためのR入門

Stan Advent Calendar 2016 の 4日目です。非マニア枠です。

　
本日のエントリーは Stan を普通に使える人、情報科学・データ科学な人を対象にしていません。少し R に慣れてきて、ベイズ統計モデリングにも興味を持ち始めた人向け。ずばりターゲットは次のようなレベルの人です。

いざ R を使って分析するときには使い方を忘れていて検索しながら使う。
dplyr の使い方に慣れておらず、めんどくさくなって結果エクセルを使ってしまう。

　
つまり私。そんな中途半端なレベルな人（私）向けに『Stan と R でベイズ統計モデリング』（以下，アヒル本とします*1）を読んでいるときにつまずくであろうところを書いておきます。アヒル本を読む・理解するには、適度に R と ggplot2 を使えるスキルが必要だからです。

　
さて，R を使った統計分析の本をいくつか読んで、簡単な分析をできるようになったという人が Stan を使ったベイズ統計モデリングをするときにつまずくのは次のようなところだと思います。

結果が出るまで待ち時間が長い。
リスト形式でデータを渡すのを戸惑う。
結果が入ったオブジェクトを表示させたらズラーッと出てきてしまい、見たい結果を見るために必死に上へスクロールする。
タイプミスする。
結果からほしいデータを取り出せない。

1. 待ち時間が長い

私の場合，簡単なモデルでも30秒くらいかかります。 f:id:hikaru1122:20161203224300p:plain

　
lm()と違って，パッと結果が帰ってこないので，学習のリズムが狂いますよね。これは仕方ない。オススメは BGM をかけることです。ガンダム UC で流れていた MOBILE SUIT です。無事ちゃんと結果が出るか，ドキドキ感が増します。

2．リスト形式にとまどう

真の R の初級者はあまりリスト形式に触れることはありません。だいたいデータフレーム形式で間に合うので。lm(y ~ x, data=d)な感じで。しかし，Stan を使う場合，

> data = list(N=row(d), A=d$A, Score=d$Score/200, M=d$M)

　
というように分析に使うデータを別に作らなければなりません。「リスト形式ってなんだよ、めんどくさい」となるんですよ*2。

　
リスト形式とは、Todo リストみたいなものだと思ってください*3。Todo リストには、いろんな種類、内容を書けます。R のリスト形式も同じ。いろいろなデータをヒトマトメにできるのです。上の list だと N には1つの数値しか入っていませんが、A や Score にはたくさんの数値が入っています。そういうのをまとめて1つのものにするのがリスト。

　
Stan を使うには専用のデータを作ってあげる必要があって、それはリスト形式なんだと納得しましょう。慣れるしかありません。

3．結果の表示が多過ぎ！

（RStudioの場合）「サンプリングもうまくいった。よし，結果見てみよう」

> fit

　
f:id:hikaru1122:20161203174841g:plain
どばっと結果が出ます。そしてパラメータを確認したいので上へ必死にスクロール。

　
こういうのを10回くらいやっていて「ひょっとしたらいい方法あるんじゃないか」と思いました。ありました。たとえば、b1・b2 というパラメータを見たい場合は次のようにします*4。

> print(fit, par="b")

Inference for Stan model: model5-6b.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

     mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
b[1] 3.57       0 0.10 3.39 3.51 3.57 3.63  3.76   976 1.01
b[2] 0.26       0 0.04 0.19 0.24 0.26 0.29  0.34  1820 1.00
b[3] 0.30       0 0.15 0.01 0.20 0.30 0.40  0.59   967 1.01

Samples were drawn using NUTS(diag_e) at Sat Dec  3 10:54:04 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

　
事後対数尤度も見たいというなら、

> print(fit, par=c("b","lp__"))

　
便利ですね。Rhat を小数点とか細かく見てみたいなという場合は例えばprint(fit, par="b", digits=3)とすればOKです。詳しい人は「当たり前」と思ってしまうレベルですが、真の初級者にはこういう tips を知っておくと，途中でめんどくさくなってやめるリスクが減ります。

4．タイプミス

私の場合、1つの Stan ファイルを実行するのに4～5回くらいエラーが出ます。本当にささいなミスで。英単語のつづりが間違っていたとか、閉じカッコが足りなかったとか、引用符がなかったとか*5。これの解決法はわかりません。誰か教えてほしいです。print を pritn とやっちゃうくらいなので。

5．ほしいデータを取り出せない。

アヒル本にあるようにrstan::extract(fit)で結果を取り出せます。取り出せますが，そのデータをいじるのには慣れが必要です。「fit の中身ってどうなっているんだろうと」と思ってstr()を使ってもよくわからない。なんか表示結果に圧倒されてしまう･･･。

> ms = rstan::extract(fit)
> str(ms)
List of 4
 $ b     : num [1:4000, 1:3] 3.48 3.57 3.56 3.52 3.62 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ iterations: NULL
  .. ..$           : NULL
 $ lambda: num [1:4000, 1:50] 3.64 3.69 3.67 3.66 3.71 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ iterations: NULL
  .. ..$           : NULL
 $ m_pred: num [1:4000, 1:50] 34 40 54 32 43 32 32 40 36 40 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ iterations: NULL
  .. ..$           : NULL
 $ lp__  : num [1:4000(1d)] 6897 6896 6898 6898 6897 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL

　
数字とアルファベットがたくさんあって，イラッときますね。でも落ち着いて。普通にfit = stan(～～～)としたときは合計4000個のMCMCサンプルが入っています。$b のところを見ましょう。推定したパラメータ b はb1, b2, b3の3つある（「1：3」のところ）。4000行3列にまとまっている状態です。なので，b1 のヒストグラムはhist(ms$b1[,1])で書けます。最大事後(MAP)推定値も max(ms$b1[,1]) で出ます。

　
ggplot を使うときはちょっとやっかいです。そもそもms$b[,1]はデータフレーム形式ではないから、そのまま ggplot には使えません。いったんデータフレーム形式に変換します。data.frame()のカッコの中が大事。データフレームを作ったときの列名をつけておきます。

> b1 = data.frame(b1 = ms$b[,1])
> str(b1)
'data.frame':  4000 obs. of  1 variable:
 $ b1: num  3.48 3.57 3.56 3.52 3.62 ...
> head(b1)
        b1
1 3.484667
2 3.572416
3 3.559215
4 3.521891
5 3.615658
6 3.657821

　
b1という列名がついた4000行1列のデータフレームになりました。ggplot でヒストグラムを作るなら

> ggplot(data=b1) + geom_histogram(aes(x=b1)) + theme_bw()

　
というように、いつもどおりの方法でできます。本エントリーの内容がわかっていれば，アヒル本のコードも読み解けるはず。コツコツ取り組んでいくと、Stan への理解が深まるだけでなく、R の操作や作図能力もぐっと上がると思います。

最後に

周りに統計に詳しくてすぐに質問できる人がいない社会科学系の社会人学生の方はいっしょにがんばっていきましょう！　そして、積極的に読書会や勉強会に参加しましょう！　意外と年齢層高いから大丈夫。

Osaka.Stan は12月23日のようです（正式な募集はまだ）。

Hommachi.Stats は12月10日です。

*1:ほんとは Rstan 本なんでしょうけども，表紙がアヒルというのが衝撃だったので。

*2:これを解消するのに、glmmStan や rstanarm や brms を使えば大丈夫ですけど･･･。

*3:たしか、DataCamp でこんな説明が出てきた気がします。

*4:ひょっとしたら本の中にそういう記述があったかもしれない。

*5:RStudioなどは補完してくれますが、それでもつい消しちゃって足りないときもあるんです。

2016-11-03

4章まで読了。『StanとRでベイズ統計モデリング』

統計分析ベイズデータ分析ノンプログラマーのためのR入門

当初1〜8章まででいいかなと目次を見てつぶやいたら，著者の方から「10章も読んだほういい」というアドバイスをいただきました(感謝です。読むモチベーションもあがります)。ということで，1〜8章＋10章を読むことを目標にします*1。自分の復習として感想を記録しておきます。

環境構築

今年からMacを使うようになり，Stanを使える環境はまだ作っていませんでした。次のページのとおりにやれば，問題なく完了しました。
https://github.com/stan-dev/rstan/wiki/Installing-RStan-on-Mac-or-Linux
※「If using g++ version 4.9 or higher〜〜」のところは不要だった。

読む前に

これまで1人でコソコソ学んできて（Kobe.Rを中心に他のRの勉強会には参加したことがあります），久保先生の『データ解析のための統計モデリング入門（みどり本）』と豊田先生の『基礎からのベイズ統計学』をなんとか読み通しています。なので，「はじめに」に書いてある「前提とする知識」はなんとかクリアしていると思います。

1〜3章は在来線特急の中で読める。

今週の出張のお供として読んでいます。1〜3章までは実際にコードを書いてみるところはほぼないので，電車の中で揺られても読めます。読んでみた第1印象は「とってもていねい！」ということ。用語が1つ1つ説明されています。「thinning って何？」「lp__って何？」とこれまで疑問に思っていたのですが，解決することができました*2。

　
3章2節には統計モデリングの手順があります。こういうことって，たぶん身近に詳しい人がいれば体得できるんでしょうが，私のような独学のものには助かります。

辛うじて突破の4章

Stanコードの書き方がていねいに解説されています。「そういう順番（やルール）で書けばいいのか」と目からウロコでした。私はプログラマでもデータサイエンティストでもなく，コードを書くときのルールをよく知らないので，こういう内容は助かります。

　
この章から実習があります。ついつい先を読んでしまいますが，ただ通読しているだけでは，わかった気になるので，実際に打ち込んで行きました。N とすべきところを n としていたり，ささいなことでエラーになるのを試行錯誤しながら，取り組んでいくことになるので，経験値が増えそうです。あと，model4-4.stan を読んでいて「このN_newって何？　突然出てきた！」と戸惑いましたが，サンプルコードを読んで解決できました。

　
52ページの解説は自分のR力のなさで，ちょっとぼんやりとした理解になりました（apply が苦手）。練習問題も（3）まではできたけど，R力のなさで（4）ができなかった…。残念。答えを見て，こんな簡単なものを解けない自分に自己嫌悪。

　
5章から先は少しずつやります。

*1:来年1月か2月までに読めたらいいなぁ，というゆっくりペースで行きます。JAGSを使った分析結果の論文を年内に書くのが優先なので。

*2:でも，正しいlp__の読み方は何だろう？　勝手に「えるぴーちょめちょめ」と読んでます。

2016-10-29

統計的因果推論勉強会第5回を開催しました！

統計的因果推論統計分析ノンプログラマーのためのR入門

経営学系統計分析エンドユーザーのため統計的因果推論勉強会第5回を行いました。最近週休3日制を打ち出した某有名大企業の新しいオフィス内での開催で，とても楽しかったです。参加者は7名（クローズド）。

speakerdeck.com

　
今回は宮川本の第5章（の特にバックドア基準）と星野本第4章「共変量選択と無視できない欠測」が勉強範囲でした。目的は2つ。1つは強く無視できる割り当て条件とバックドア基準が似たようなものであることを理解すること。もう1つはRでの実習。Rnotebookを使うと便利でした。配布する資料もカンタンに作れるので。

試用したものはこちらにあります。

github.com

　
Rでの実習は傾向スコアマッチングを星野本の巻末付録に載っているもので写経（使用したデータはいろいろねつ造したフィクションのものです），IPW法を岩波データサイエンス第3巻のサポートページに載っているもので写経しました。

　
c統計量（このcは大文字のほうがいいのか，小文字なのかわかりません･･･）については，私たちはお手軽に rms パッケージの lrm() でやったほうがラクです。「C」と記載されているところがc統計量です。

　
f:id:hikaru1122:20161029205730p:plain