Knowledge As Practice

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

purrr の map を使って検定を繰り返す

ずっとKJ法とかポストイットとかを使ってアンケートデータと格闘しているのでまったく統計分析はしていませんでした。少し目処が立ったので、次の研究のための予備調査を兼ねて昔のデータを見直すことにしました。

 
こんなデータがあるとします。

f:id:hikaru1122:20180918210302p:plain

 
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になっても大丈夫。もっと上手な(オシャレな)方法があるのだと思います。あれば教えてください。

参考(やっぱりあった)

https://sebastiansauer.github.io/looping-purrr/

【追記あり】マルコフ連鎖の収束診断の基準と書き方メモ

マルコフ連鎖の収束って、どうやって書けばいいのかわからないので、ざっと調べることにしました。事の発端は、ShinyStanで収束診断(DIAGNOSE)をするとき3つの指標が示されるのですが、3つの指標全部が None じゃないといけないか、と r-wakalang で質問したことです。回答もいたただけて、もう少し書き方も含めて調べてみようとなりました。

 
r-wakalang はこちら。 github.com

 
ShinyStan の使い方は↓の33ページ以降が役立ちます。

www.slideshare.net

 
cinii で検索したり、twitter で流れてきたものをかいつまんで読んでいます。全部で4つです。

 
●浅田他(2014)では「収束診断方法は複数あり~」「Gelman-Rubin統計量(R-hat)は1.1未満となっているときに、連鎖が定常状態に収束していると判断した」「Gelman-Rubin統計量(R-hat)の点推定値は~~のいずれも1.1以下であった」とあり、R-hatの計算方法として古谷(2010)、マッカーシー(2007)を参考せよとある。

 
●Namba et al.(2018)には「The value of Rhat for all parameters equalled 1.0, indicating convergence across the four chains.」とあります。

 
●清水(2017)には「サンプリング回数を1万回、バーンイン期間を5000回、マルコフ連鎖の数を4に指定して推定したところ、すべてのパラメータのR_hatが1となり、収束が確認された」とあり、脚注に「1に近いと収束したと判断される指標。1.05以下が目安とされる(Gelman, Carlin, Stern, &Rubin, 2004)」とBDAの第2版が引用されていました。

 
ちなみに、こちらではBDA第2版で「1.1より下」という基準が引用されています。ん~、とりあえず出張から帰ったら最新版のBDA3を確認する。 takehiko-i-hayashi.hatenablog.com

 
●渡辺他(2017)では「計算の結果,Gelman-Rubin 統計量はすべてのパラメータに対して1.01以下となり,各マルコフ連鎖が定常状態に収束していることを目視で確認した」とあり、BDA第3版が引用されていました。

 
ざっと調べたところ書き方も含め、1.1未満、1.0になってる、1.05以下と思ったよりいろいろありました。みどり本も帰ったら確認してみよう。

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

  • 作者: Andrew Gelman,John B. Carlin,Hal S. Stern,Donald B. Rubin,David B. Dunson
  • 出版社/メーカー: Chapman and Hall/CRC
  • 発売日: 2011/12/15
  • メディア: ハードカバー
  • この商品を含むブログを見る

 

参考文献

浅田正彦, 長田穣, 深澤圭太, & 落合啓二. (2014). 状態空間モデルを用いた階層ベイズ推定法によるキョン (Muntiacus reevesi) の個体数推定. 哺乳類科学, 54(1), 53-72.

Namba, S., Kabir, R. S., Miyatani, M., & Nakao, T. (2018). Dynamic displays enhance the ability to おりdiscriminate genuine and posed facial expressions of emotion. Frontiers in Psychology, 9, 672.

清水裕士. (2017). 二者関係データをマルチレベル分析に適用した場合に生じる諸問題とその解決法. 実験社会心理学研究, 56(2), 142-152.

渡辺 憲, 高麗 秀昭, 小林 功, 柳田 高志, 鳥羽 景介, 三井 幸成, 階層ベイズモデルを用いた丸太の天然乾燥における乾燥時間の推定および丸太の諸形質が乾燥性に及ぼす影響の評価, 木材学会誌, 公開日 2017/03/30, Online ISSN 1880-7577, Print ISSN 0021-4795, https://doi.org/10.2488/jwrs.63.63, https://www.jstage.jst.go.jp/article/jwrs/63/2/63_63/_article/-char/ja

 

追記(2018年6月16日)

BDA3を確認。Rhat ついては287ページに「閾値として1.1の設定で一般的によしとしてきた」と書いてあった。犬4匹本では,Rhat について(Gelman-Rubin統計量として)184ページで言及されている。みどり本では,206ページで言及されている。犬4匹本もみどり本もRhatの目安は1.1。豊田先生編著『基礎からのベイズ統計学』には190ページでGelman(1996)を引用して,Rhatが1.1ないしは1.2より小さいことが収束の判断として紹介されている。

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

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

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

 
というわけで,収束の判断のRhatについては1.1より小さいことが目安でよさそう(もちろんRhatだけで判断していいかどうかは別問題ということで)。

博士後期課程に入って2年経ちました

「もうブログ終わっちゃったんですか?」ときかれたので書きます。書かなかったのは、ここ最近は質的研究に取り組んでいて、特に書くネタがなかったからです。たぶん、かなりRの操作を忘れています・・・。

 
さて、JAIST 東京の社会人博士後期課程に入ってから丸2年経ちます。私は情報科学ではなくて、知識科学系(ザックリいうと経営学のほう)に属しています。あまり社会人学生の情報がない、という声も聞くので、すこし感想を書いておきます。

 
1.研究はどう進めているのか
2.ここ2年間の成果報告

の2点です。

 

1.研究はどう進めているのか

 
基本は平日の夜と土日にやっています。毎日、論文書いたり研究できたりするわけではありません。しかし、1段落でも論文を読んだり、1文字でも何か書く努力はしています。社会人院生は研究が進まないときに「仕事に逃げる」ということができ(仕事が夜遅くなるときがあるし、バタバタするのも日常茶飯事)、「仕事が忙しくて研究が進まなかった」と言い訳が立ちます。

 
こういうことはなるべく避けたいなと思い、防止策として今年初めから社会人学生どうしで「研究進捗部」を作り、みんなでコツコツ進捗を出しています。

 
f:id:hikaru1122:20180330123152j:plain
↓のアイデア丸パクリ参考にしました。

www.buzzfeed.com

 

2.ここ2年間の成果報告

 
1年目(2016年)は特に成果なし。授業を取ったり(5~6科目を取らないといけない)、研究に関係ありそうな書籍や論文にあたってみたり、新生活への対応に追われていた気がします。月に1回のゼミ前に有志で行う統計分析の勉強会で、半年間ほど統計的因果推論を学んだことがいちばん印象に残っています。

 
※資料はこのへんにあります。

https://www.slideshare.net/hikaru/1-62476733

Presentations by Hikaru Goto // Speaker Deck

 
2年目(2017年)は国際会議(大会? conference)に4つ投稿して、1つめはOK、2つめはNGが出ました。3つめと4つめは結果待ち、またはあさってに投稿です。OKが出たところは小規模ながら、質の高い場所で参加してよかったと思います。「研究がんばったで賞」もいただくことができ(なぜかその場で賞金400ドルまで出た)、黒字になって帰国しました。

f:id:hikaru1122:20170803141519j:plainf:id:hikaru1122:20170803183130j:plain
海外の学会ってお金持ちなんでしょうか。

https://www.irssm8atyonsei.com

 
NGが出たところは査読コメントに「理解してくれてない」というネガティブな印象を持ってしまったものの、英語力不足・準備不足があいまった結果だと思っています。

https://www.ieseg.fr/en/faculty-and-research/research-events/servsig-2018/

 
国内の大会でも1つポスター発表しました。個人的にはポスター発表のほうが好きです。質問の受け答えやフォローが口頭発表より濃くできるので。

第5回国内大会(2017/3/27-28) - サービス学会(Society for Serviceology)

 
あと、上海にある大学でも発表を行いました。話すことで自分のアイデアも固まるので、こういうのって大事ですね。

f:id:hikaru1122:20180330113121j:plain
上海、勢いある。

他に、すうがくぶんかの統計分析のセミナーに出たり、国際公共経済学会のサマースクールに1日だけ参加したり、行動計量学会春の合宿セミナーに1日だけ行ったり、マーケティング学会の研究会をのぞいてみたりしました。

 
結果、次の修了要件のうち、ほぼ半分は達成できたかと思います。大学には長期履修(目標は4年)を申請済みなので、残り期間もあと半分です。

 
国際会議で発表→○
授業11単位→△(あと1つ)
査読付論文→×
副テーマ研究→×

 
3年目の2018年度が正念場ですね。

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 アカウントでいつも接続しています。

https://rstudio.cloud/

f:id:hikaru1122:20171225223852j:plain

 
以前作ったプロジェクトを開いてみると……

f:id:hikaru1122:20171225232751j:plain

 
できた!!

 
まあ,当たり前ですね。

 
しかし,困ったことがあります。いざスクリプトを書こうとすると,キーボードの表示によって,画面半分がおおわれてしまいます。そのため,実行時にコンソールが見えず,いちいちキーボードを隠せねばなりません。こんな感じ。

f:id:hikaru1122:20171225224136j:plain

 
物理のキーボードがあればなんとかなるんじゃないかと思い,Smart Keyboard を買ってきました。今日もアップルストア心斎橋は混雑してましたね。セッティングしたらこんな感じです。

f:id:hikaru1122:20171225224810j:plain

 
これで通常の使い方(?)に近づいた感があります。気をつけたい点は,スクリプト実行のショートカットは 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」です。パワーパックにするとスニペットの機能が使えます。

 
【あとで動画を載せる】

 

www.alfredapp.com

 

2.BetterTouchTool で Rmarkdown を使いやすく

 
Rmarkdownのコードチャンク(下のようなもの)を出すときのショートカットも同じように楽な設定をしています。具体的には、指3本で下に動かすとチャンクが挿入されます。

 
【あとで動画を載せる】

 
Cmd + Option + I はなんとか覚えられそうですが、このショートカットを使うときは、いちいちキーボードに目を落とせねばならず(タッチタイプがそのレベル)、それがストレスだったのです。

 
もちろん、これもスニペットツールでやればいいのですが、なぜか私は BetterTouchTool のジェスチャーでチャンクが出るようにしています。どうしてそうしたのは思い出せません。

 

www.boastr.net

 

3.いちいち library()するのが面倒なときにWコロン

 
「このパッケージのこの関数だけを1回だけ使いたい」ということがありますよね。たぶん。そういうときに::(Wコロン)です。

 
どのパッケージかはっきりさせるため ::+関数があるのでしょうが、私の場合はどちらかというと不純な動機で「ちょっと1度だけ、各変数の平均と標準偏差をざっと見たい」なんてときに psych::describe(hogehoge) とかやっています。

 

4.パッケージ管理は GUI

 
install.packages() はほぼ使わなくなりました。基本、RStudio でパッケージ管理しています。

 
f:id:hikaru1122:20171221195943j:plain

 
「Install」ボタンを押して、インストールしたいパッケージを頭から数文字入れるとサジェスチョンを出してくれます。だから、スペルミスが減りました。あと、だいたい月に1回くらい「Update」ボタンを押します。一気にアップデートされるときは気持ちいいです。

 
実はここ最近、Jupyter Notebook も使おうとしていたのですが、パッケージのインストール・管理の機能がないので、使うのをやめたくらい。

 

5.データを渡すときは .xlsx のファイルで

 
なんとか周りにRを使う人を増やそうと、今年は月1のゼミ前に1時間ほど有志で統計分析の勉強会を3月から11月まで毎月行いました。経営学に近い分野の社会人学生は文系が多いので、あまり量的分析を好みません。しかし、それではいけないので、みんなでがんばりました。

 
はじめは例題として分析するデータをCSVやテキストファイルで渡していました。しかし、Windows の人や Mac の人が混在しているので、それぞれ対応する文字コードにしたファイルを渡すことになり、まず「文字コード???」となってつまずきます。

 
だったら、「もうエクセルでいいじゃない、RStudio のGUIだとエクセルファイル読み込めるし(正確にはread_excel()を使っているけど)、文字コードの問題も起きない」ということでエクセルで配布するようになりました。

f:id:hikaru1122:20171221211241j:plain

 
それ以来、文字コードのトラブルは起きず、インポート前に自主的にエクセルファイルの中身を見て(だいたい Office が入っていて、無意識にダブルクリックするみたい)、分析前に確認する習慣も付きそうです。

 
以上、ご意見があるかと思いますが、極私的なRの使い方でした。もっといいものがあれがお教えください。よろしくお願いします。

階層モデルでロバストな推定を試す

【追記】
「外れ値の確認」の項目を追記しました。


 
以前『データ分析をマスターする〜』の9章を階層モデルを使って分析しました。しかし,ちょっと気になることがあったので,再分析してみます。なお,本エントリーは Stan Advent Calendar 2017 の9日目です。

 
今回主に使うパッケージは次の2つです。

library(tidyverse)
library(brms)
library(rethinking)
# devtools::install_github("rmcelreath/rethinking",ref="Experimental")  # 必要であれば。

 
使用したデータは『データ分析をマスターする〜』のサポートサイトにありますが,使って見るならちゃんと買ったほうがいいかなぁと思います(私は買って読んでます)。  
 

以前の分析の課題 〜外れ値の存在〜

 
この前は都道府県別の女性労働力率と合計特殊出生率のデータを使って,シンプソンのパラドックスベイズ推定を使った階層モデル(マルチレベル分析)でどう解消するか考えました。↓がこの前のエントリーです。

 

hikaru1122.hatenadiary.jp

 
データの特徴は「全体で見ると負の相関だけど,年代別(1980年と2000年)で見ると正の相関が出る」というものです。散布図を再掲します。

 
f:id:hikaru1122:20171209021854p:plain

 
前回,傾きはほぼ同じだから,変量効果は切片だけ考えました。それはいいのですが,気になっていて見なかったことにしようとしたことが1つ。外れ値の存在です。

 

ggplot(d, aes(x=factor(year), y=tfr)) + geom_boxplot(aes(fill=factor(year))) + 
geom_jitter() + coord_flip() + 
theme(legend.position = "none")

 
f:id:hikaru1122:20171209021910p:plain

 
上図を見ると,各年代で合計特殊出生率の極端に高い都道府県と極端に低い都道府県が存在します。下表から、どの年代も最大値は沖縄(都道府県コード47)で,最小値は東京都(都道府県コード13)です。

 

library(knitr) # kable() で下の表を作るためにマークダウン出力用
d %>% group_by(year) %>% 
summarise(max(tfr), which.max(tfr), min(tfr), which.min(tfr)) %>% kable()
year min(tfr) which.min(tfr) max(tfr) which.max(tfr)
1980 1.44 13 2.38 47
2000 1.07 13 1.82 47

 

外れ値の確認

 
外れ値が存在するときは t 分布を使うとよいと犬4匹本に書いてあります(466ページ)ので,以下,t分布を使った外れ値があるデータにも対応できるロバストな分析を行います。

 
ところで、今回のケースが外れ値といえるのか少し確かめたいと思います。外れ値を調べるために Smirnov-Grubbs 検定(すみるのふぐらぶす)があります。この検定はデータが正規分布から発生している仮定しています。しかし、今回はt分布から発生すると考えていきたいので、自由度46(=47-1)のt分布の95%信頼区間を見て判断してみます。

t = abs(qt(0.05/2, 47-1))
SE1980 = sd(d$tfr[d$yearID==1])/sqrt(47)
SE2000 = sd(d$tfr[d$yearID==2])/sqrt(47)

> # 1980年の下側
> mean(d$tfr[d$yearID==1]) - t*SE1980
[1] 1.789009
> # 1980年の上側
> mean(d$tfr[d$yearID==1]) + t*SE1980
[1] 1.868012
> # 2000年の下側
> mean(d$tfr[d$yearID==2]) - t*SE2000
[1] 1.433888
> # 2000年の上側
> mean(d$tfr[d$yearID==2]) + t*SE2000
[1] 1.512069

 
1980年の95%信頼区間は1.79〜1.87,2000年の95%信頼区間は1.43〜1.51でした。よって,外れ値が存在すると言えそうです。上記の計算はこちらを参考にしました。

www2.hak.hokkyodai.ac.jp

 

t 分布を使ったロバストな推定

 

ロバストとは?

 
ロバストとは「頑健」という意味です。外れ値があっても,データのメインボディの中心やばらつきなどを妥当に推定する方法を「ロバスト推定」と言います(藤澤 2017)。

 
t分布には3つのパラメータがあります。自由度ν(にゅー),中心μ(みゅー),スケールσ(しぐま)です。μやσは平均と標準偏差の記号として使われますが,ここではt分布なので平均ではなく中心 center,標準偏差ではなくスケール scale と書きました(BDA3, 437ページ)。

 
t 分布の感じをつかむには Unadon さんのこちらのエントリーがわかりやすいです。

mrunadon.github.io

 
Stan のおすすめ事前分布を見ると,t(3, 0, 1)のようですが,上の箱ひげ図を見ると幅広そうなので,今回はできるだけ幅が広く,裾が長いt分布の t(1, 0, 10)でやってみます。なお,BDA3には「自由度が1や2のt分布はあんまり現実的ではない」とありますが,今回は実験ということで目をつぶります(437ページ)。

 
Stan おすすめ事前分布についてはこちらの説明があります。

github.com

 
Stan のコードは自前では書けませんので,今回も brms パッケージを使います。変量効果は前と同じ(半コーシー分布に従っている)。単回帰分析なので,切片とflpの傾きがt分布に従っていると考えます。モデルのイメージは下のような感じ(かっこいいグラフィカルな表現はできなかった)。

 
f:id:hikaru1122:20171209022057j:plain

 
コードは次のとおりです。

 

library(brms)
options(mc.cores = parallel::detectCores())

prior2 = c(set_prior("cauchy(0, 3)", class = "sd", group = "prefID"),
           set_prior("cauchy(0, 3)", class = "sd", group = "yearID"),
           set_prior("student_t(1,0,10)", class = "b"),
           set_prior("student_t(1,0,10)", class = "Intercept"))

kekka2 = brm(tfr ~ 1 + flp + (1 | prefID) + (1 | yearID), data = d, iter = 5000, 
             family = "student", prior = prior2,
             control = list(adapt_delta = 0.9, max_treedepth = 15))

 
ちゃんと指定した事前分布になっているかコードを確認します。

 

> stancode(kekka2)
(中略)

  // priors including all constants 
  target += student_t_lpdf(b | 1,0,10);  ←ココ
  target += student_t_lpdf(temp_Intercept | 1,0,10);  ←ココ
  target += student_t_lpdf(sigma | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10); 
  target += gamma_lpdf(nu | 2, 0.1)
    - 1 * gamma_lccdf(1 | 2, 0.1); 
  target += cauchy_lpdf(sd_1 | 0, 3) ←ココ
    - 1 * cauchy_lccdf(0 | 0, 3); 
  target += normal_lpdf(z_1[1] | 0, 1); 
  target += cauchy_lpdf(sd_2 | 0, 3) ←ココ
    - 1 * cauchy_lccdf(0 | 0, 3); 
  target += normal_lpdf(z_2[1] | 0, 1); 

(後略)

 

分析結果

 
分析結果は次のとおりです。無事収束して,心なしか有効サンプルサイズが増えたようです。また前回と違って切片の固定効果は0をまたがなくなり,最後の警告の数は減りました。都道府県の変量効果と傾きの固定効果がほぼ0なのは同じです。

 

> kekka2
 Family: student(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.09     0.17       1183 1.01

~yearID (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.96      0.82     0.15     3.11       1162 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     1.67      0.69     0.22     3.16       1287 1.00
flp          -0.00      0.00    -0.01     0.01        829 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.04      0.01     0.03     0.05       1538 1.00
nu       13.89     11.73     2.69    46.47       5153 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 7 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 
> 

 
警告メッセージについては,divergent transitions の数が小さいのと Rhat を見ると収束しているので無視します(アヒル本176ページ)。

 
傾きは前回と同じように変量効果も固定効果も0に近いです。切片の変量効果が減ってしまっています(前のモデルが1.22)が,切片の固定効果の95%信用区間は0をまたがなくなりました(切片の有効サンプルサイズも増えています)。図で確認します(太い線が50%信用区間,細い線が90%信用区間です。点は事後中央値)。

 
f:id:hikaru1122:20171209021942p:plain

 

モデル比較

 
WAIC は前回のモデルに及びませんでした。残念。

 

> waic(kekka1, kekka2)
                   WAIC    SE
kekka1          -268.04 16.62
kekka2          -260.81 18.53
kekka1 - kekka2   -7.23  3.47

 

まとめ

 
今回は「外れ値があるので,ロバストな推定をしよう」とチャレンジしてきました。結果はイマイチで,たった94行のデータに適用するのはあんまり意味がないのでしょう。しかし,いつもと違った事前分布に触れることができました。また Stan コードが書けなくても,パッケージで比較的に Stan に触れられる例にもなったかなと思います。

 
間違いなどがあれば,ぜひお教えください🙏
いつかお会いする機会があれば🍣を進呈いたします。

 

参考にした書籍

以下,今回,参考にした書籍です。どれも難しいですね…。

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

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

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

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

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

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

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

  • 作者: Andrew Gelman,John B. Carlin,Hal S. Stern,David B. Dunson,Aki Vehtari,Donald B. Rubin
  • 出版社/メーカー: Chapman and Hall/CRC
  • 発売日: 2013/11/05
  • メディア: ハードカバー
  • この商品を含むブログを見る

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)

RStudio Cloud による社会人学生への統計教育の実践

Hijiyama.R the Final で LT をしてきました。LT のタイトルは「RStudio Cloud による社会人学生への統計教育の実践」です。RStudio Cloud を使えば,みんな同じ環境でRに触れることができるよ,という内容でした。

    f:id:hikaru1122:20171126125144j:plain

 
月1ペースで社会人学生で統計分析の勉強会をしているのですが,次のような理由で分析する環境を整えるのが結構大変だと思っていました。

  • PC環境がバラバラで環境構築であきらめることも。

  • アカウント名が日本語になっていることがよくある。

  • RStudioやパッケージを全員がちゃんとインストールするのが難しい…。

 
一人ひとりフォローすればなんてことはないのですが,勉強会は月1で1時間程度で時間がないのです。どうしたもんかと悩んでいたところに救世主が!

 
f:id:hikaru1122:20171126220106p:plain

 
使い方はとても簡単。右上の Login から「Log in with Google」で利用できるようになります。ログイン後,しばらく待っているとスタンバイ完了。New Project でプロジェクトを作って,いつものようにRStudioを使うだけ。

 
f:id:hikaru1122:20171126220613p:plain

 
動作は重いです。データは右下のペインの「Upload」からアップします。R スクリプトを書いて保存して,URLを共有すれば,すぐに他の人も使えます。違うアカウントでログインすると,右下に作成したスクリプトやデータがありますし,共有前にパッケージをインストールしていたら,そのパッケージも使えます(しかし,rstan とはダメです)。

 
f:id:hikaru1122:20171126220841p:plain

 

まとめ

  • R を使ってこんなふうに分析できるよという紹介用レベルには使える。
  • 重い。

以上です。

補足

 
kazutan 先生によれば,次のとおりだそうです。

補足2

 
今日はいつもと違って,Jupyter Notebook + RISE でプレゼンしました。

f:id:hikaru1122:20171126223710p:plain

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