2017-12-21

極私的Ｒの使い方５選（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の使い方でした。もっといいものがあれがお教えください。よろしくお願いします。

2017-12-09

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

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

　
以前『データ分析をマスターする〜』の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ページ）ので，以下，ｔ分布を使った外れ値があるデータにも対応できるロバストな分析を行います。

　
ところで、今回のケースが外れ値といえるのか少し確かめたいと思います。外れ値を調べるために Smirnov-Grubbs 検定（すみるのふぐらぶす）があります。この検定はデータが正規分布から発生している仮定しています。しかし、今回はｔ分布から発生すると考えていきたいので、自由度46（=47-1）のｔ分布の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）。

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

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

mrunadon.github.io

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

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

github.com

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

　
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)

作者: 畑農鋭矢,水落正明
出版社/メーカー: 有斐閣
発売日: 2017/10/12
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

ロバスト統計: 外れ値への対処の仕方 (ISMシリーズ:進化する統計数理)

作者: 藤澤洋徳
出版社/メーカー: 近代科学社
発売日: 2017/07/31
メディア: 単行本
この商品を含むブログを見る

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

作者: John K. Kruschke,前田和寛,小杉考司
出版社/メーカー: 共立出版
発売日: 2017/07/22
メディア: 大型本
この商品を含むブログを見る

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

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (8件) を見る

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)

作者: Richard McElreath
出版社/メーカー: Chapman and Hall/CRC
発売日: 2016/02/19
メディア: ハードカバー
この商品を含むブログを見る

2017-11-26

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

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

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 先生によれば，次のとおりだそうです。

んー、これalpha版で多分test期間中なはず。
— kazutan v3.4.2 (@kazutan) 2017年11月26日

補足2

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

f:id:hikaru1122:20171126223710p:plain

2017-10-28

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

Kobe.R ベイズデータ分析

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

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

作者: John K. Kruschke,前田和寛,小杉考司
出版社/メーカー: 共立出版
発売日: 2017/07/22
メディア: 大型本
この商品を含むブログを見る

　
　
難しいです。懇切丁寧なのは確かですが，難しいです･･･。両章ともテーマは階層モデリングでした。参加者の方から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

2017-10-23

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

ベイズデータ分析統計分析

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

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

作者: 畑農鋭矢,水落正明
出版社/メーカー: 有斐閣
発売日: 2017/10/12
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

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

合計特殊出生率（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)

作者: 畑農鋭矢,水落正明
出版社/メーカー: 有斐閣
発売日: 2017/10/12
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

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

作者: John K. Kruschke,前田和寛,小杉考司
出版社/メーカー: 共立出版
発売日: 2017/07/22
メディア: 大型本
この商品を含むブログを見る

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

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (8件) を見る

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

作者: 清水裕士
出版社/メーカー: ナカニシヤ出版
発売日: 2014/10/01
メディア: 単行本
この商品を含むブログを見る

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

作者: Richard McElreath
出版社/メーカー: Chapman and Hall/CRC
発売日: 2016/02/19
メディア: ハードカバー
この商品を含むブログを見る

2017-09-24

第1回 Kobe.R annex（犬4匹本輪読会）開催しました！

ベイズデータ分析 Kobe.R 参加報告

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章の便利関数をパッケージ化されました！

参加させていただき、ありがとうございました。

8章の便利コマンドですが、パッケージとしてgithubに上げました。
devtools::install_github("kyusque/DBDA2E-utilities", args = "--no-multiarch")
— ヨゼ (@kyusque) 2017年9月24日

LT枠（分析ネタなんでも）

状態空間モデルについて
共変量バランシング傾向スコアを使った統計的因果推論について

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

2017-09-20

大きい 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本でA４・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（日本）