Rによるカテゴリカルデータ分析事例(3) ～{rms}パッケージによる順序ロジスティック回帰～
●2016年1月3日 追記
本エントリー中で使っている makedummies() がさらに便利に使えるようになっています。
アップデート版の使い方はコメント欄を参照してください。
以上、追記終わり。
----
前回までに、あるサービスを受けた500人の顧客の満足度および変数ｘは、性別や年齢層というデモグラフィック的な要素によって違いを受けるとは言えないことがわかりました。唯一、変数ｘは性別によって違いがあることがわかりました。
前回のエントリーはこちら↓
hikaru1122.hatenadiary.jp
変数ｘというとわかりにくいので、顧客のサービス参加度（以下、参加度）としておきます。本当はもっと複雑で、Kobe.R では分析そっちのけでこの議論が少しあったのですが、ここでは簡潔に「参加度」としておきます。要は、あるサービスに関わろうとする度合いが強い顧客ほど満足度が高いのではないか、ということを調べます。
知りたいこと
顧客のサービス参加度は満足度（cs）に正の影響を与えているかどうか。
→参加度が高いほど、満足度も高くなるのでは？
分析の手順
１．まず相関をチェック。
２．順序ロジスティック回帰分析で、参加度の満足度への影響を調べる。
使用するデータ
使うデータは前回のエントリーと同じです。ダウンロードは↓から。
https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed
実証分析
０．データの確認 まずクロス表で参加度と満足度の状況を確認します。たしかに、参加度が大きいほど、満足度も高い人が多そうです。
> xtabs(~x + cs, data =d)
| 満足度1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| 参加度1 | 2 | 2 | 4 | 51 | 26 | 12 | 2 |
| 2 | 1 | 1 | 7 | 14 | 38 | 15 | 10 |
| 3 | 0 | 0 | 4 | 6 | 61 | 34 | 9 |
| 4 | 0 | 1 | 2 | 6 | 28 | 63 | 27 |
| 5 | 2 | 0 | 2 | 2 | 3 | 15 | 50 |
１．相関係数をチェック
回帰分析の前に、相関係数を調べておきたいと思います。特に相関がなければ、回帰分析をする必要がないと考えるからです。参加度は1～5の順序尺度、満足度は1～7の順序尺度です。連続値ではないので、厳密には普通の相関係数（ピアソンの積率相関係数）で調べることはできません。
5つや7つの順序尺度なら、連続値とみなして別にいいんじゃないかという意見もあります*1。私もそうしたいです、楽だし。でも、今回は学びのために、少し厳密になってみます。
3つ以上カテゴリを持つ順序カテゴリカル変数には、ポリコリック相関係数を調べます。cs とｘを順序カテゴリカル変数に変換しておき、データの構成を確認します。
> d$cs2 <- ordered(d$cs) > d$x2 <- ordered(d$x) > str(d) 'data.frame': 500 obs. of 6 variables: $ gender : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 2 1 1 2 ... $ age_seg: Factor w/ 5 levels "1","2","3","4",..: 3 4 4 4 3 1 2 4 3 3 ... $ cs : int 5 7 4 6 6 6 5 6 7 5 ... $ x : int 3 2 1 1 3 4 2 5 2 3 ... $ cs2 : Ord.factor w/ 7 levels "1"<"2"<"3"<"4"<..: 5 7 4 6 6 6 5 6 7 5 ... $ x2 : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 3 2 1 1 3 4 2 5 2 3 ...
ポリコリック相関係数を出すには、{polycor}パッケージのpolychor()を使います。
> library(polycor) > polychor(d$cs2, d$x2) [1] 0.5913401
まずまずの値でしょうか。hetcor()だとp値を知ることができます。出てくる数値は同じく、0.59です。今回は有意（p値 1.665e-12）のようです。
> hetcor(d$cs2, d$x2) Two-Step Estimates Correlations/Type of Correlation: d$cs2 d.x2 d$cs2 1 Polychoric d.x2 0.5913 1 Standard Errors: d$cs2 d.x2 0.02967 Levels: 0.02967 n = 500 P-values for Tests of Bivariate Normality: d$cs2 d.x2 1.665e-12
２．順序ロジスティック回帰
正の相関があることが確認できました。顧客の参加度（1～5）がサービスの満足度（1～7）に影響を与えているかどうかを調べます。どちらも間隔尺度なので、本来は普通の回帰分析を使うことはできません*2。応答変数が順序カテゴリカル変数の場合、順序ロジスティック回帰が適切です。説明変数もカテゴリカルデータです。ダミー変数を作って、分析をしてみることにします。
まずデータの準備です。少し回りくどいですが、初めにcs2 ではなく、cs を抜き出したのは、makedummies() が順序カテゴリカル変数をすべてダミー変数化してしまうためです*3。顧客の参加度のみダミー変数を作りたいため、順序尺度ではない cs を抜き出しました。
> library(dplyr) ## データｄから cs と x2 を抜き出します。 > d2 <- select(d, cs, x2) ## 顧客の参加度x2 のダミー変数を作ります。 > library(makedummies) > d2 <- makedummies(d2) ## 最後に 満足度 cs を順序カテゴリカルにします。 > d2$cs <- ordered(d2$cs) > str(d2) 'data.frame': 500 obs. of 5 variables: $ cs : Ord.factor w/ 7 levels "1"<"2"<"3"<"4"<..: 5 7 4 6 6 6 5 6 7 5 ... $ x2_2: num 0 1 0 0 0 0 1 0 1 0 ... $ x2_3: num 1 0 0 0 1 0 0 0 0 1 ... $ x2_4: num 0 0 0 0 0 1 0 0 0 0 ... $ x2_5: num 0 0 0 0 0 0 0 1 0 0 ...
ダミー変数のベースとなるのは「顧客の参加度 1」です。では、順序ロジスティック回帰分析を行います。使用するパッケージは{rms}で、関数はorm()です*4。
> orm(cs ~ ., data = d2) ## ↑の ~. は cs 以外の全部の変数を使う指示の省略形です。 Logistic (Proportional Odds) Ordinal Regression Model orm(formula = cs ~ ., data = d2) Frequencies of Responses 1 2 3 4 5 6 7 5 4 19 79 156 139 98 Model Likelihood Discrimination Rank Discrim. Ratio Test Indexes Indexes Obs 500 LR chi2 213.58 R2 0.364 rho 0.569 Unique Y 7 d.f. 4 g 1.555 Median Y 5 Pr(> chi2) <0.0001 gr 4.736 max |deriv| 1e-05 Score chi2 215.96 |Pr(Y>=median)-0.5| 0.302 Pr(> chi2) <0.0001 Coef S.E. Wald Z Pr(>|Z|) y>=2 3.4341 0.4644 7.39 <0.0001 y>=3 2.8328 0.3564 7.95 <0.0001 y>=4 1.6431 0.2283 7.20 <0.0001 y>=5 -0.1027 0.1860 -0.55 0.5809 y>=6 -1.9960 0.2191 -9.11 <0.0001 y>=7 -3.8202 0.2575 -14.83 <0.0001 x2_2 1.1229 0.2786 4.03 <0.0001 x2_3 1.6875 0.2598 6.50 <0.0001 x2_4 2.7115 0.2729 9.94 <0.0001 x2_5 4.4556 0.3537 12.60 <0.0001
係数 x2_2 ～ x2_5 を確認すると、顧客の参加度1をベースにして、2～5と参加度が上がるごとに、係数が増えています。よって、顧客の参加度は満足度に正の影響を与えていると言えるので、顧客の参加を促す施策を行ったほうがよいことがわかります。
今回、参考にしたもの
- 作者: 藤井良宜,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2010/04/22
- メディア: 単行本
- クリック: 13回
- この商品を含むブログ (3件) を見る
- 作者: 太郎丸博
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2005/07
- メディア: 単行本
- 購入: 1人 クリック: 26回
- この商品を含むブログ (12件) を見る
- 作者: 山本勲
- 出版社/メーカー: 中央経済社
- 発売日: 2015/10/21
- メディア: 単行本
- この商品を含むブログ (1件) を見る
相関係数の補足
なお、ノンパラメトリックなスピアマンの順位相関係数を知りたい場合は、次のようにします。ポリコリック相関係数のほうが高い値が出やすいように思います。今回も少し高かったです（今回のポリコリック相関係数は0.59）。
> cor(d$cs, d$x, method = "spearman") [1] 0.5690423
順序ロジスティック回帰の補足
順序ロジスティック回帰は{MASS}パッケージのpolr()がよく例として出されます。以下、polr() でやってみます。なお、polr() の結果はp値が出力されません*5。そのため、{texreg}パッケージのscreenreg()でp値を確認します。
> library(MASS) > kekka <- polr(cs ~ ., data = d2) ## ↑の ~. は cs 以外の全部の変数を使う指示の省略形です。 > library(texreg) > screenreg(kekka) ## カッコ内は標準誤差です。 Re-fitting to get Hessian =========================== Model 1 --------------------------- x2_2 1.12 *** (0.28) x2_3 1.69 *** (0.26) x2_4 2.71 *** (0.27) x2_5 4.46 *** (0.35) --------------------------- AIC 1345.59 BIC 1387.73 Log Likelihood -662.79 Deviance 1325.59 Num. obs. 500 =========================== *** p < 0.001, ** p < 0.01, * p < 0.05
結果は同じです。こっちのほうがシンプルな結果が出るので、やはりpolr()が便利なのでしょうか。
Rによるカテゴリカルデータ分析事例(2) ～(1) の補足：等分散性を確認していたか？～
前回のエントリーを書いたあと、現在は手に入りにくい本『Rによる統計解析の基礎』（私は買いました）の最新版『Rによる保健医療データ解析演習』を確認し、さらに「マイナーだけど最強の統計的検定 Brunner-Munzel 検定」を読むと、ウィルコクソンの順位和検定（マン・ホイットニーのU検定）*1とクラスカル・ウォリス検定には等分散性の確認が必要のようです。
Rによる統計解析の基礎 (Computer in Education and Research)
- 作者: 中澤港
- 出版社/メーカー: ピアソンエデュケーション
- 発売日: 2003/10
- メディア: 単行本
- 購入: 2人 クリック: 28回
- この商品を含むブログ (7件) を見る
『Rによる保健医療データ解析演習』は↓の「ダウンロード」セクションでゲットできます。とても勉強になります。
http://minato.sip21c.org/msb/index.html
念のため、分散の同質性を確認します。この確認にはフリグナー・キリーン検定を行います。デフォルトで使えるfligner.testか {coin}パッケージの fligner_test を使います。せっかくなので、前回と同様、{coin} パッケージを使ってみます。データは前回のエントリーと同じものです。帰無仮説は「各群（性別または年齢層）の分散は同じ」。有意水準を5%とします。
次の分析結果から、帰無仮説を棄却することはできませんでした。よって、前回のエントリーの分析を行ってOKです。
> fligner_test(cs ~ gender, data =d) Asymptotic Two-Sample Fligner-Killeen Test data: cs by gender (0, 1) Z = -1.0815, p-value = 0.2795 alternative hypothesis: true ratio of scales is not equal to 1 > fligner_test(cs ~ age_seg, data =d) Asymptotic K-Sample Fligner-Killeen Test data: cs by age_seg (1, 2, 3, 4, 5) chi-squared = 5.6545, df = 4, p-value = 0.2265 > fligner_test(x ~ gender, data =d) Asymptotic Two-Sample Fligner-Killeen Test data: x by gender (0, 1) Z = 0.00056251, p-value = 0.9996 alternative hypothesis: true ratio of scales is not equal to 1 > fligner_test(x ~ age_seg, data =d) Asymptotic K-Sample Fligner-Killeen Test data: x by age_seg (1, 2, 3, 4, 5) chi-squared = 6.177, df = 4, p-value = 0.1863
*1:等分散性が確認できれなければ、 Brunner-Munzel 検定を使います。
Rによるカテゴリカルデータ分析事例(1) {coin}パッケージによるウィルコクソンの順位和検定とクラスカル・ウォリス検定
先日（2015年12月26日土曜日）に開催された Kobe.R 第23回で発表を行いました。当日はちょっと風邪気味だったため、あまり準備ができず、テキトーな内容になってしまいました…。そのため、ここで補足をしたいと思います。
発表タイトルは「Rによるカテゴリカルデータ分析事例」。現在、私が分析を進めているデータを使った内容です。扱ったデータの構成は
- 性別（0が男性、1が女性の名義尺度）
- 年齢層（1が20代、2が30代、3が40代、4が50代、5が60代の名義尺度）
- あるサービスを受けた満足度（1～7の順序尺度）
- 変数x（1～5の順序尺度）
というもの。実際のデータセットはもっと多くて50個ほど変数があります。サンプルサイズは500です。今回使用するデータは↓においてあるので、興味がある方はやってみてください。 https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed
分析によって調べたいことは次の3つです。
- 満足度および変数xは性別によって異なるか。
- 満足度および変数xは年齢層によって異なるか。
- 変数xは満足度に正の影響を与えているか（変数xが高いほど、満足度も高いか）。
今回は１と２を扱います。３は次回のエントリーに回します。では、やってみます。有意水準は5%とします。
０．使用するデータの準備と確認
データの構成は次のようになっています。
> d <- blogyou.data.20151229 > str(d) 'data.frame': 500 obs. of 4 variables: $ gender : int 0 0 0 0 0 1 1 0 0 1 ... $ age_seg: int 3 4 4 4 3 1 2 4 3 3 ... $ cs : int 5 7 4 6 6 6 5 6 7 5 ... $ x : int 3 2 1 1 3 4 2 5 2 3 ...
性別と年齢層を因子型に変えて、もう1度構成を確認します。
d$gender <- as.factor(d$gender) d$age_seg <- as.factor(d$age_seg) > str(d) 'data.frame': 500 obs. of 4 variables: $ gender : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 2 1 1 2 ... $ age_seg: Factor w/ 5 levels "1","2","3","4",..: 3 4 4 4 3 1 2 4 3 3 ... $ cs : int 5 7 4 6 6 6 5 6 7 5 ... $ x : int 3 2 1 1 3 4 2 5 2 3 ...
１．満足度および変数xは性別によって異なるか？
性別は2値の名義尺度、満足度と変数xは順序尺度です。まずはクロス表を作ってデータを眺めてみます。
●満足度と性別の関係
> xtabs(~gender + cs, data=d)
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| 男性(0) | 1 | 3 | 6 | 38 | 66 | 59 | 35 |
| 女性(1) | 4 | 1 | 13 | 41 | 90 | 80 | 63 |
ん～、どうでしょう？ 女性のほうが高い満足度にある人が多いような…。図示してみます。
> library(lattice) > histogram(~factor(cs)|gender, data = d) ## cs は離散値のため、factor() で因子型に直しています。以下、x も同様です。
では、実際に性別によって満足度の差があるか検定します。上図より正規分布に従っているようには見えないので、ノンパラメトリック検定のウィルコクソンの順位和検定（マン・ホイットニーのU検定）を行います。帰無仮説は「性別によって差がない（差が０）」、対立仮説が「性別によって差がある（差が０でない）」です。
さて、ここで{coin}パッケージを使います。ウィルコクソンの順位和検定を行うには、wilcox.test や {exactRankTests} の wilcox.exact を使うのですが、{exactRankTests}をインストールして使おうとしたら「{coin}を使え！」とコンソール画面で怒られたので、{coin} を使うことにします。
## {coin}パッケージがない場合は、インストールしてください。 library(coin) > kekka <- wilcox_test(cs ~ gender, distribution = "exact", data = d) ## distribution = "exact"とすると、パソコンによっては結果が出るまで時間がかかります。 > kekka Exact Wilcoxon-Mann-Whitney Test data: cs by gender (0, 1) Z = -1.0953, p-value = 0.2737 alternative hypothesis: true mu is not equal to 0
帰無仮説は棄却できませんでした。したがって、性別によって満足度に違いがあるとは言えませんでした。次に、変数ｘと性別の関係と調べます。
●変数ｘと性別の関係
先ほどと同じ検定を行いますが、その前にクロス表とヒストグラムで確認します。
> xtabs(~gender + x, data =d)
| 1 | 2 | 3 | 4 | 5 | |
|---|---|---|---|---|---|
| 男性(0) | 52 | 30 | 51 | 53 | 22 |
| 女性(1) | 47 | 56 | 63 | 74 | 52 |
histogram(~factor(x)|gender, data =d)
女性は特に変数ｘが5の人が多いように見えます。検定します。
> kekka2 <- wilcox_test(x ~ gender, distribution = "exact", data = d) > kekka2 Exact Wilcoxon-Mann-Whitney Test data: x by gender (0, 1) Z = -2.2005, p-value = 0.02769 alternative hypothesis: true mu is not equal to 0
変数x は性別によって違いがあることがわかりました。
２．満足度および変数xは年齢層によって異なるか。
年齢層は5つの水準がある名義変数です。年齢層と満足度および変数xの関係をクロス表とヒストグラムで確認します。
> xtabs(~age_seg + cs, data =d)
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| 20代(1) | 0 | 1 | 4 | 18 | 32 | 25 | 20 |
| 30代(2) | 4 | 2 | 3 | 18 | 34 | 22 | 17 |
| 40代(3) | 0 | 1 | 6 | 19 | 30 | 25 | 19 |
| 50代(4) | 0 | 0 | 0 | 14 | 31 | 37 | 18 |
| 60代(5) | 1 | 0 | 6 | 10 | 29 | 30 | 24 |
> xtabs(~age_seg + x, data =d)
| 1 | 2 | 3 | 4 | 5 | |
|---|---|---|---|---|---|
| 20代(1) | 24 | 23 | 17 | 19 | 17 |
| 30代(2) | 20 | 22 | 24 | 19 | 15 |
| 40代(3) | 25 | 10 | 24 | 29 | 12 |
| 50代(4) | 18 | 16 | 21 | 29 | 16 |
| 60代(5) | 12 | 15 | 28 | 31 | 14 |
histogram(~factor(cs)|age_seg, data =d) histogram(~factor(x)|age_seg, data =d)
どうも年齢層による違いは満足度、変数ｘともに見られそうにありません。検定します。水準が3つ以上なので、先ほどのウィルコクソンの順位和検定（マン・ホイットニーのU検定）ではなく、クラスカル・ウォリス検定を行います。
## kruskal_test() は {coin} パッケージにあります。 kekka3 <- kruskal_test(cs ~ age_seg, data = d) kekka4 <- kruskal_test(x ~ age_seg, data = d)
満足度および変数xの年齢層による違いは次の結果から差があるとはいえませんでした。
> kekka3 Approximative Kruskal-Wallis Test data: cs by age_seg (1, 2, 3, 4, 5) chi-squared = 8.0373, p-value = 0.0895 > kekka4 Asymptotic Kruskal-Wallis Test data: x by age_seg (1, 2, 3, 4, 5) chi-squared = 5.4405, df = 4, p-value = 0.245
次回に続きます。今回、参考にしたものはこちらの2冊です。
- 作者: 藤井良宜,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2010/04/22
- メディア: 単行本
- クリック: 13回
- この商品を含むブログ (3件) を見る
- 作者: 太郎丸博
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2005/07
- メディア: 単行本
- 購入: 1人 クリック: 26回
- この商品を含むブログ (12件) を見る
Rでダミー変数をつくる。makedummies 関数で。
●2015年12月28日 追記
makedummies 関数の作者の方から、コメントで情報をいただきました。
以下のように、 install_github で
簡単に makedummies が使えるようになっています。
## devtools がなければインストール→ install.packages("devtools") library(devtools) install_github("toshi-ara/makedummies")
●2016年1月3日 追記
makedummies がさらに便利になりました。詳細は↓です。
http://hikaru1122.hatenadiary.jp/entry/2015/12/30/213032
以上、追記おわり。
今回のモチベーション
最近、ずっとカテゴリカルデータをいじっていて、ダミー変数を作る機会が出てきました。Rでダミー変数を作るのには、パッケージがあったはずだと思い、調べていたら、こちらの2つがみつかりました。
dummyVars についてはこちら。
estrellita.hatenablog.com
dummy.data.frame についてはこちら。
qiita.com
makedummies 関数
ただ、これらはダミー変数を作るときに基準となる変数もそのまま作ります。なので、基準となる変数のダミーはいりません。勝手に基準となるダミーを作らない関数がないかと探ししたところ、ありました。
それがこちらの makedummies 関数です。
github.com
作ってくださった方に感謝です。データフレームしか受け付けないので、その点だけ注意すればダミー変数が簡単につくれます。github のソースをコピーして、適当なファイルに保存して、適当なフォルダに入れておきます。今回は、makedummies.R というファイルにして、getwd() で表示されるフォルダに入れています。
source("makedummies.R") > d <- read.delim("clipboard") > str(d) 'data.frame': 18 obs. of 2 variables: $ y: int 5 7 4 6 6 6 6 6 7 5 ... $ x: int 1 2 3 1 2 3 1 2 3 1 ... > head(d) y x 1 5 1 2 7 2 3 4 3 4 6 1 5 6 2 6 6 3
x のダミー変数を作りたいと思います。x が2のときに「1」となるダミー、x が3のときに「1」となるダミーを作ります。1は基準なのでダミーはいりません。ダミー変数としたいものは、因子型である必要があるので、因子型に直します。
> d$x <- factor(d$x) > str(d) 'data.frame': 18 obs. of 2 variables: $ y: int 5 7 4 6 6 6 6 6 7 5 ... $ x: Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
順序付きのカテゴリカル変数でもOKです。あとは makedummies 関数に入れるだけ。自動的に変数名がつきます。無事、基準となる変数はダミーが作られていません。ちなみに、複数の因子型変数がある場合、全部ダミー変数を作ります。
> d2 <- makedummies(d) > head(d2) y x_2 x_3 1 5 0 0 2 7 1 0 3 4 0 1 4 6 0 0 5 6 1 0 6 6 0 1
私には重宝する関数です。改めて感謝。
ちなみに、カテゴリカルデータの分析はこちらの2冊が勉強なりました。
- 作者: 藤井良宜,金明哲
- 出版社/メーカー: 共立出版
- 発売日: 2010/04/22
- メディア: 単行本
- クリック: 13回
- この商品を含むブログ (3件) を見る
- 作者: 太郎丸博
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2005/07
- メディア: 単行本
- 購入: 1人 クリック: 26回
- この商品を含むブログ (12件) を見る
Rによる級内相関係数（ICC）の2つの出し方
2～3ヶ月前から、マルチレベル分析*1を学んでいます。分析自体はできます（便利なことに、コード1行でできるので）。しかし、手持ちのデータでマルチレベル分析を使うべきなのかどうか、迷います。そんなときは級内相関係数（ICC）を確認するのが一般的なようです*2。例によって、ウェブ上で情報収集しました*3。
使うデータ
あるサービスの500人分の満足度を１～７で集計したもの。性別（男性0、女性1）も記録している。
> cs <- read.delim("clipboard") > head(cs) gender manzoku 1 0 5 2 0 6 3 0 4 4 0 5 5 0 4 6 1 6
ICCの出し方
１．lmerTest または lme4 パッケージのlmerを使う。
切片だけの変量効果を調べます。summary()を使って、Random effectsのVarianceを計算してICCを出します。すなわち、
gender の Variance ÷（gender の Variance＋ Residual の Variance）
です。今回はとっても低い値になります。マルチレベル分析をする必要はなさそうです。
> library(lmerTest) > icc_by_lmerTest <- lmer(manzoku ~ 1 + (1|gender), data = cs) > summary(icc_by_lmerTest) ※前略 Random effects: Groups Name Variance Std.Dev. gender (Intercept) 2.497e-14 1.58e-07 Residual 1.742e+00 1.32e+00 ※後略
２．ICCパッケージのICCestを使う。
x がグループ、y が対象となるデータ、data がデータフレームです。x は因子型でないと怒られます。
> library(ICC) > icc_by_ICCest <- ICCest(x = gender, y = manzoku, data = cs) Warning message: In ICCest(x = gender, y = manzoku, data = cs) : 'x' has been coerced to a factor > icc_by_ICCest <- ICCest(x = factor(gender), y = manzoku, data = cs) > icc_by_ICCest $ICC [1] -0.002186884
なんだか、lmerTest を使った計算方法と結果が違います。でも、とても小さい値なので、やはりマルチレベル分析をする必要はなさそうです。なお、ICC の数値が違います。計算方法が違うので、どちらが正しいものかはわかりません…。今回はどちらも小さな値になったので、よかったものの、微妙な値だったらどうしようか…。
番外 psychパッケージのICCを使う。
psych パッケージの説明を読んだんですが、データの形式がちょっと違っていたので、面倒そうでやめました。
参考にしたページ
こちらのブログのエントリーには多くの資料のリンクがまとめてあり、助かりました。 http://mizumot.com/lablog/archives/179
なお、マルチレベル分析の基本的なところはこちらがたいへん勉強になります。
www.slideshare.net
読了 『完全独習 ベイズ統計学』（2015）
個人的に大きなイベントが終わり、うまくいきそうなので、学びと研究を再開です。
リハビリとして、今年11月に出版された『完全独習 ベイズ統計学入門』を読んで、ブログを書きます。本の詳細は下をクリックしてください。
- 作者: 小島寛之
- 出版社/メーカー: ダイヤモンド社
- 発売日: 2015/11/20
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (3件) を見る
『図解・ベイズ統計「超」入門』より、もう少し計算などを使ってベイズ統計を学んでいく本のように感じました。したがって、『図解・ベイズ～』などのベイズ統計を紹介する本を1冊でも読んだことがある人は苦しまずに読みきれると思います。
図解・ベイズ統計「超」入門 あいまいなデータから未来を予測する技術 (サイエンス・アイ新書)
- 作者: 涌井貞美
- 出版社/メーカー: SBクリエイティブ
- 発売日: 2013/12/18
- メディア: 新書
- この商品を含むブログ (13件) を見る
『完全独習 ベイズ～』は面積図などを使って、あの手この手でわかりやすく説明をしてくれます。『基礎からのベイズ統計学』がさっぱりわからないという方は、この『完全独習 ベイズ～』を読んでみるといいのではないかと思いました。逆に『基礎からの～』を読んだ方は必要がないというか、『完全独習 ベイズ～』の第2部を立ち読みでくらいでいいんじゃないかと。
基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門
- 作者: 豊田秀樹
- 出版社/メーカー: 朝倉書店
- 発売日: 2015/06/25
- メディア: 単行本
- この商品を含むブログ (3件) を見る
今回の本でもっともよかったのは最後のブックガイドです。ベイズ統計への理解を深めるために次はこういう本を読めばいいのだな、と教えてくれます。その中の1冊にチャレンジしたいと思います（RやStan関連の本の紹介はありません）。
とにかくこれから勉強したい方にこういうのがあるよ、とオススメできる本が増えたのは嬉しいです。
統計学エンドユーザーのための JAGS 入門
岩波データサイエンス vol.1 のどこかのページに「初心者はJAGS（BUGS）」と書いてあったので*1、こつこつ学び始めています。世の中は Stan 大人気ですが、気にしない…。まず JAGS に慣れてから Stan をやります…。本エントリーは、前回のエントリーと多くがかぶっています。
今回は、次のページを参考にさせていただいています。
生態学のデータ解析 - R で JAGS (rjags 編)
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html
生態学のデータ解析 - JAGS 雑
http://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html
分析に使用したのは、JAGS 4.0.0 および rjags 4.4 です。
JAGS と rjags のインストール
JAGS と rjags はそれぞれ別にインストールします。JAGS 3.0 系 をインストールしていて、rjags 4.4 にしていた場合は動きません。JAGS も４にする必要があります。
JAGS はこちらから .exe をダウンロードできます。
http://sourceforge.net/projects/mcmc-jags/files/
rjags はいつもの install.packages でOKです。
分析の手順
１．データの準備
今回は先ほどの「RでJAGS（rjags 編）」のデータを使用しました。上記ページにある csv ファイルから、ページのとおりにデータをセットします。JAGSに渡すデータはリスト形式である必要があります。
２．モデルを作る
「RでJAGS（rjags 編）」にあるモデルを自分なりにいじってみました。いまのJAGSマニュアルによると、行末にセミコロンは不要のようです。モデルを作る前に図解しました。これも岩波データサイエンスにそうしろ、と書いてあったように思います。分布の形はテキトーです*2。
さて、モデルです。ロジスティック回帰でリンク関数はロジットリンクです。詳しくはみどりぼん第６章を参照してください。
Model <-" model{ for (i in 1:N) { Positive[i] ~ dbin(p[i], Total[i]) logit(p[i]) <- z[i] z[i] <- beta[1] + beta[2] * Logdose[i] } beta[1] ~ dnorm(0, tau) beta[2] ~ dnorm(0, tau) tau ~ dgamma(1.0E-6, 1.0E-6) } "
追記2015年11月20日
そういえば、JAGS の dnorm の分散は逆数でした。
tau <- 1/(sd*sd)
sd ~ dgamma(1.0E-3, 1.0E-3)
がより正しい書き方なのかな…？
また、twitter にてアドバイスをいただきました。
一様分布や半コーシー分布のほうがよい、とのことです。
ありがとうございます。
いちばんの変更点は tau の事前分布です。一度、同じように dgamma(1.0E-3, 1.0E-3) でやってみたのですが、うまく行きませんでした。代わりに、dgamma(1.0E-6, 1.0E-6) を使いました。
その理由は2つあります。
（１）JAGS 4.0 マニュアルに載っていた例で dgamma(1.0E-6, 1.0E-6) が使われていたこと
（２）前回のエントリーで dgamma(1.0E-6, 1.0E-6) を使ってやり直したら、Stan とほぼ同じ結果が得られたこと
です。1.0E-3 と 1.0E-6 で結果が変わるのってシビア。事前分布って難しい。
モデルを書いたら、書き出します。
writeLines(Model, con ="testModel.txt") # 拡張子は特に不要ですが、いちおう。
３．MCMCする
モデルを jags オブジェクトにしていきます。
# rjags のロード。 library(rjags) # jags モデルを作る。 testModel <- jags.model(file = "testModel.txt", data = list.data, n.chain = 3)
jags.model の file = は別ファイルである必要があります。なので、先ほど書きだした次第です。次にバーンイン、そして本格的なサンプリングを行います。n.chain は２以上でないと、後で Gelman & Rubin 収束診断指標のグラフが描けないようなので、３にしています。初期値は設定していません。まだ理解できていないので。
では、本格的なMCMCサンプリングを行います。
# バーンイン期間のための試運転。 update(testModel, n.iter = 1000) # 本番。MCMC サンプルを kekka に収める。 kekka <- coda.samples( testModel, variable.names = c("beta[1]", "beta[2]", "tau"), n.iter = 20000)
coda.samples の variable.names では結果で見たいパラメータを文字列として指定します。これで、下ごしらえは終わりました。rjags パッケージで使うのは、jags.model, update, coda.samples の３つだけです。writeLines は普通に使えます。
結果を見る。
収束診断をします。codamenu() でメニューを選びながら、確認します。引数はいらないです。メニューで「2→kekka（MCMCサンプリングを格納したオブジェクト）→2→2」とするとGelman & Rubin の指標が見られます。1.1 より下なので、収束しているようです。他にもcodamenu()でいろいろ指標を選べます。
Potential scale reduction factors: Point est. Upper C.I. beta[1] 1.01 1.02 beta[2] 1.01 1.02 tau 1.00 1.00
図でも確認します。plot(kekka)です。
推定値は summary(kekka) で確認します。
> summary(kekka) Mean SD Naive SE Time-series SE beta[1] -5.99179 1.10916 0.0045281 0.053419 beta[2] 1.18754 0.19647 0.0008021 0.009398 tau 0.06042 0.07121 0.0002907 0.001165
「RでJAGS（rjags 編）」ページに書いてある「真の値」は、 beta[1] = -5.0、beta[2] = 1.0 です。けっこう近づきました。tau だけはダメでした。図示します。
少しずつ、JAGS に慣れてきたようです。今後、みどりぼんや岩波データサイエンスの例にも取り組んでいきたいと思います。
参考文献
JAGS マニュアル http://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/
rjags マニュアル https://cran.r-project.org/web/packages/rjags/rjags.pdf
生態学のデータ解析 - R で JAGS (rjags 編) http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html
生態学のデータ解析 - JAGS 雑 http://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html
岩波データサイエンス
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2015/10/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (5件) を見る
みどりぼん
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (25件) を見る
DBDA２
Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan
- 作者: John Kruschke
- 出版社/メーカー: Academic Press
- 発売日: 2014/11/17
- メディア: ハードカバー
- この商品を含むブログを見る
