今回は「階層ベイズモデル」をRで推定する方法について解説します。
「階層ベイズモデル」という用語は、少なくとも社会学をやっている学生からすると馴染みが薄いですが、これは要するに、「マルチレベルモデル(Multilevel model)」あるいは「階層線形モデル(Hierarchical linear model)」あるいは「混合効果モデル(Mixed effect model)」と呼ばれているモデルを、ベイズ的な方法によって推定することを意味します。数式についてはのちほど。
ふだんStataやR、SPSSといった統計ソフトでマルチレベルモデルを推定する場合、基本的には最尤法(ML)あるいは制限付き最尤法(REML)といった方法でパラメータを推定します。しかし、ランダム効果などのパラメータが多くなってきたり、非線形モデルを推定したりするときには、ちゃんとパラメータが推定できなかったりします。
そこで、MCMCなど、たくさんのサンプルを生成して、そこからパラメータの分布を求めようというのが「階層ベイズモデル」のねらいです。
なお、本ページは麦山のウェブサイトに置かれています。例示に使用したデータはFinch, W. Holmes, Jocelyn E. Bolin, and Ken Kelley. 2014. Multilevel Modeling Using R. CRC Press.で使用されているデータであり、こちらからダウンロードすることができます。Languageというデータがそれです。
今回使用するパッケージは以下のとおりとなります。入っていない場合はインストールしてください。
#マルチレベルモデルを推定する
#install.packages("lme4")
#lme4は作者のこだわり?のためかp値が算出されない(信頼区間は出力される)ため、p値を表示したい場合はこのパッケージを入れる必要がある
#install.packages("lmerTest")
#ブートストラップ法を行う際に使用する
#install.packages("boot")
#いろいろと便利なプロットを出力できる
#install.packages("sjPlot")
#install.packages("sjmisc")
#install.packages("sjlabelled")
#Rstanを経由して階層ベイズモデルを推定する
#install.packages("brms")
#dplyrやtidyr, ggplot2など一連のパッケージ
#install.packages("tidyverse")
library(lme4)
library(lmerTest)
library(optimx)
library(boot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(brms)
library(tidyverse)
今回は1週間ごとにテストのスコアを測定した次のデータを使ってみます。いくつか変数が入っていますが、ここでは個人(id)、学校(school)、時点(time)、スコア(score)の4つの変数だけを使うことにします。
# データの読み込み
widedata <- read.csv("Language.csv", header = TRUE, sep=",")
widedata2 <- widedata[,c(1:8)]
##生徒数が10人未満の学校を除外し、かつ500ケースを無作為に抽出したデータを作成
set.seed(111)
subdata <- widedata2 %>%
subset(school != 1 & school != 7 & school != 17 & school != 24 & school != 28 & school != 29) %>%
sample_n(size = 500)
# wideデータをlongデータに変換し、変数"id"を作成
longdata <- subdata %>%
gather(time, score, LangScore1:LangScore6) %>%
rename("id" = ID) %>%
arrange(id)
# 文字列の削除
longdata$time <-
as.numeric(gsub("LangScore","",longdata$time))
dim(longdata)
## [1] 3000 4
kable(head(longdata, n = 15))
id | school | time | score |
---|---|---|---|
2 | 4 | 1 | 211 |
2 | 4 | 2 | 213 |
2 | 4 | 3 | 214 |
2 | 4 | 4 | 215 |
2 | 4 | 5 | 211 |
2 | 4 | 6 | 223 |
3 | 4 | 1 | 223 |
3 | 4 | 2 | 229 |
3 | 4 | 3 | 228 |
3 | 4 | 4 | 232 |
3 | 4 | 5 | 232 |
3 | 4 | 6 | 247 |
12 | 4 | 1 | 214 |
12 | 4 | 2 | 216 |
12 | 4 | 3 | 194 |
このデータは、同一個人が複数回測定されたデータとなっています。各個人のscoreの推移を確認してみましょう。
# 10ケースを無作為抽出
set.seed(123)
rand <- sample(unique(longdata$id), 10)
# 無作為抽出した各個体についてスコアの推移をみる
ggplot(data = longdata[longdata$id %in% rand,], aes(x = time, y = score)) +
geom_point() +
ylim(150, 250) +
geom_smooth(method = "lm", formula = y ~ x) +
facet_wrap("id")
これをみると、おおむねどの個人も、timeの経過にしたがってscoreが上向きになっていることがわかります。
この傾向を捉えるために、同一個人を複数回にわたって観察したというデータの構造を無視して、通常の線形モデルを当てはめてみましょう。モデルは以下のとおり。
# 線形モデルの推定
m1 <- lm(score ~ time, data = longdata)
# 結果を出力
summary(m1)
##
## Call:
## lm(formula = score ~ time, data = longdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.445 -10.296 1.938 11.704 43.321
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 193.1469 0.6803 283.92 <2e-16 ***
## time 3.3830 0.1747 19.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.34 on 2998 degrees of freedom
## Multiple R-squared: 0.1112, Adjusted R-squared: 0.1109
## F-statistic: 375.1 on 1 and 2998 DF, p-value: < 2.2e-16
このモデルもまぁまぁ悪くないかもしれませんが、同一個人を複数回観察した、というデータの特徴を活かすことで、もっと現実をうまく反映したモデルを作れるかもしれません。そこでマルチレベルモデルを使います。
マルチレベルモデルは、個々の観察値が上位の単位にネストしている、という構造を反映するためのモデルです。たとえば以下のようなモデルを考えることができます。
上の式を1つの式にまとめると、次のようになります(分布についての情報は省略)。
マルチレベルは基本的にはOLSのように解析的に解くことはできず、最尤法などの数値計算によって解くのがふつうです。ここではこの方法の詳細には立ち入りませんが、興味のあるかたはRaudenbush, Stephen W., and Anthony S. Bryk. 2002. Hierarchical linear models: Applications and data analysis methods. Sage. などを参照ください。
まず切片のみのモデルを推定してみましょう。
# 切片のみモデルの推定
model0 <- lmer(score ~ 1 + (1 | id),
data = longdata,
REML = FALSE,
na.action = na.omit)
# 推定結果を出力
summary(model0)
## summary from lme4 is returned
## some computational error has occurred in lmerTest
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + (1 | id)
## Data: longdata
##
## AIC BIC logLik deviance df.resid
## 23460.7 23478.7 -11727.3 23454.7 2997
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5066 -0.5052 0.0738 0.5800 3.1822
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 206.76 14.379
## Residual 93.43 9.666
## Number of obs: 3000, groups: id, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 204.9873 0.6668 307.4
切片のみのモデルを推定するとはすなわち,従属変数の総分散のうち,どの程度が生徒レベル(レベル2)の分散であり,どの程度が時点×生徒レベル(レベル1)の分散であるかを分割することを意味しています.
ちなみに,ICC(Intra-Class Correlation, 級内相関係数)を計算することもできます。ICCとは、総分散のうちレベル2の分散が占める割合がどの程度であるかを示す指標で、0以上1以下の値を取ります。これは以下のように定義されます。
このモデルについて、ICCは以下のように計算できます。
icc <- 206.76/(206.76 + 93.43)
icc
## [1] 0.6887638
ついで、説明変数を1つ追加したモデルは以下のとおりです。
# 独立変数timeを追加
model1 <- lmer(score ~ 1 + time + (1 | id),
data = longdata,
REML = FALSE,
na.action = na.omit)
# 推定結果を出力
summary(model1)
## Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
## approximations to degrees of freedom [lmerMod]
## Formula: score ~ 1 + time + (1 | id)
## Data: longdata
##
## AIC BIC logLik deviance df.resid
## 22063.0 22087.0 -11027.5 22055.0 2996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9193 -0.5075 0.0400 0.5668 5.0363
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 213.43 14.609
## Residual 53.38 7.306
## Number of obs: 3000, groups: id, 500
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 193.1469 0.7207 678.3000 268.01 <2e-16 ***
## time 3.3830 0.0781 2500.0000 43.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## time -0.379
以上では切片の分散のみを考慮していましたが、独立変数の効果(ここでは時点、言い換えればスコアの伸びかた)が生徒によって異なるという可能性があります。この場合、傾きの部分にもランダム効果を認めるモデルを推定することになります。式は以下のとおりとなります。