1 Introduction はじめに

今回は「階層ベイズモデル」を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というデータがそれです。

2 Preparing Packages パッケージの準備

今回使用するパッケージは以下のとおりとなります。入っていない場合はインストールしてください。

#マルチレベルモデルを推定する
#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)

3 Preparing Data データの準備

今回は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

4 Standard Linear Model 通常の線形モデル

このデータは、同一個人が複数回測定されたデータとなっています。各個人の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が上向きになっていることがわかります。

この傾向を捉えるために、同一個人を複数回にわたって観察したというデータの構造を無視して、通常の線形モデルを当てはめてみましょう。モデルは以下のとおり。

\[ Y_{ij} = \beta_0 + \beta_1X_{ij} + \varepsilon_{ij}, \ \ \ \varepsilon_{ij} \sim N(0, \sigma^2) \]

# 線形モデルの推定
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

このモデルもまぁまぁ悪くないかもしれませんが、同一個人を複数回観察した、というデータの特徴を活かすことで、もっと現実をうまく反映したモデルを作れるかもしれません。そこでマルチレベルモデルを使います。

5 Multilevel Model マルチレベルモデル

5.1 What is Multilevel Model? マルチレベルモデルとは

マルチレベルモデルは、個々の観察値が上位の単位にネストしている、という構造を反映するためのモデルです。たとえば以下のようなモデルを考えることができます。

\[ Y_{ij} = \beta_{0j} + \beta_{1j}X_{ij} + r_{ij} \\ \beta_{0j} = \gamma_{0} + u_{0j}, \\ \beta_{1j} = \gamma_{1} + u_{1j}, \\ \]

\[ \begin{aligned} \begin{pmatrix} u_{0j} \\ u_{1j} \end{pmatrix} \sim \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_{u_0}^2 & \\ \sigma_{u_1, u_0} & \sigma_{u_1}^2 \end{pmatrix} \right) \end{aligned} \\ r_{ij} \sim N(0, \sigma_{r}^2) \] これは、切片がグループ\(j\)によって異なること、および、変数\(X_{ij}\)の傾きがグループ\(j\)によって異なることを認めたモデルとなります。ここで切片および傾きのばらつきは2変量正規分布にしたがうと仮定しておきます(この仮定は後から緩めることもできます)。この場合、新たに投入された\(u_{0j}, u_{1j}\)のことをランダム効果(Random effect)、前者をとくにランダム切片(Random intercept)、後者をとくにランダム傾き(Random slope)といいます。ここではランダム効果の間に相関を認めていますが、これを0に設定することも可能です。

上の式を1つの式にまとめると、次のようになります(分布についての情報は省略)。 \[ Y_{ij} = \gamma_{0} + \gamma_{1}X_{ij} + u_{0j} + u_{1j}X_{ij} + r_{ij} \]

5.2 Maximum Likelihood 最尤法(最適化)での推定

マルチレベルは基本的にはOLSのように解析的に解くことはできず、最尤法などの数値計算によって解くのがふつうです。ここではこの方法の詳細には立ち入りませんが、興味のあるかたはRaudenbush, Stephen W., and Anthony S. Bryk. 2002. Hierarchical linear models: Applications and data analysis methods. Sage. などを参照ください。

5.2.1 Random-Intercept Model ランダム切片モデル

まず切片のみのモデルを推定してみましょう。

\[ Y_{ij} = \gamma_{0} + u_{0j} + r_{ij}, \\ u_{0j} \sim N(0, \sigma_u^2) \\ r_{ij} \sim N(0, \sigma_r^2) \]

# 切片のみモデルの推定
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以下の値を取ります。これは以下のように定義されます。 \[ \mathrm{ICC} = \frac{\sigma_r^2}{\sigma_u^2 + \sigma_r^2} \]

このモデルについて、ICCは以下のように計算できます。

icc <- 206.76/(206.76 + 93.43)
icc
## [1] 0.6887638

ついで、説明変数を1つ追加したモデルは以下のとおりです。 \[ Y_{ij} = \gamma_{0} + \gamma_1X_{ij} + u_j + r_{ij}, \\ u_j \sim N(0, \sigma_u^2) \\ r_{ij} \sim N(0, \sigma_r^2) \]

# 独立変数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

5.2.2 Random-Intercept and Random-Slope Model ランダム切片・ランダム傾きモデル

以上では切片の分散のみを考慮していましたが、独立変数の効果(ここでは時点、言い換えればスコアの伸びかた)が生徒によって異なるという可能性があります。この場合、傾きの部分にもランダム効果を認めるモデルを推定することになります。式は以下のとおりとなります。

\[ Y_{ij} = \beta_{0} + \beta_1X_{ij} + \beta_2 X_{ij} + u_{0j} + u_{1j} X_{ij} + r_{ij}, \\ \begin{aligned} \begin{pmatrix} u_{0j} \\ u_{1j} \end{pmatrix} \sim \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_{u_0}^2 & \\ \sigma_{u_1, u_0} & \sigma_{u_1}^2 \end{pmatrix} \right) \end{aligned} \\ r_{ij} \sim N(0, \sigma_r^2) \] このモデルでは、切片のランダム効果と、傾きのランダム効果の間に相関を認めています。例えば切片の値が大きいほど傾きが大きい(あるいは小さい)といった関連があるとき、この項は有意な値を示すことになります。もちろんこれを0に設定することも可能です。

実際のコードは以下のとおりです。

# timeの係数にランダム効果を認める
model2 <- lmer(score ~ 1 + time + (1 + time | id),
                     data      = longdata,
                     REML      = FALSE,
                     na.action = na.omit)

# 結果を出力
summary(model2)
## Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
##   approximations to degrees of freedom [lmerMod]
## Formula: score ~ 1 + time + (1 + time | id)
##    Data: longdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  21857.7  21893.7 -10922.8  21845.7     2994 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2522 -0.4737  0.0197  0.5210  5.2730 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 343.714  18.540        
##           time          2.687   1.639   -0.76
##  Residual              43.971   6.631        
## Number of obs: 3000, groups:  id, 500
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 193.1469     0.8739 500.0000  221.03   <2e-16 ***
## time          3.3830     0.1020 500.0000   33.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## time -0.716

ここまで、3つのモデルを推定してきましたが、どのモデルが最もデータの特徴をうまく捉えているかを見るために、モデルの適合度を比較してみます。これは、anovaなどの関数を使うことで行うことができます。

# モデル間でLoglikelihood testをして適合度を検定
anova(model0, model1, model2)
## Data: longdata
## Models:
## object: score ~ 1 + (1 | id)
## ..1: score ~ 1 + time + (1 | id)
## ..2: score ~ 1 + time + (1 + time | id)
##        Df   AIC   BIC logLik deviance   Chisq Chi Df Pr(>Chisq)    
## object  3 23461 23479 -11727    23455                              
## ..1     4 22063 22087 -11028    22055 1399.67      1  < 2.2e-16 ***
## ..2     6 21858 21894 -10923    21846  209.31      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

なお、切片と傾きの相関を0に設定したモデルを推定する場合は、以下のように記述します。

model2nocor <- lmer(score ~ 1 + time + (1 | id) + (-1 + time|id),
               data      = longdata,
               REML      = FALSE,
               na.action = na.omit)

summary(model2nocor)
## Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
##   approximations to degrees of freedom [lmerMod]
## Formula: score ~ 1 + time + (1 | id) + (-1 + time | id)
##    Data: longdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  22045.2  22075.2 -11017.6  22035.2     2995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0925 -0.4872  0.0313  0.5426  5.0025 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 235.08   15.332  
##  id.1     time          1.00    1.000  
##  Residual              48.55    6.968  
## Number of obs: 3000, groups:  id, 500
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 193.14693    0.74453 556.10000  259.42   <2e-16 ***
## time          3.38297    0.08689 556.10000   38.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## time -0.300

5.3 Baysian Estimation ベイズでの推定

ではこのモデルをベイズで推定してみます。このとき、C++のコンパイラが必要となります。OSによってどのようなアプリケーションを入れておくかが違います。

  • Mac OSを使用している場合:Xcodeというアプリケーションをインストール。
  • Windowsを使用している場合:Rtoolsというアプリをインストール。

5.3.1 Estimation without Prior 事前分布を指定しないで推定

では、先ほどのModel 2をベイズ推定で求めてみましょう。brmsパッケージに含まれているbrmという関数を使います。まずコードは以下のとおりです。モデル部分およびデータの指定の部分は先のlme4と同じで、わかりやすいです。計算の負荷を減らすために繰り返しの回数をちょっと少なめに設定しています。

# brm関数を使って推定
model2b <- brm(score ~ 1 + time + (1 + time | id), 
               data    = longdata,
               prior   = NULL, #事前分布を指定。NULLと記述した場合は一様分布
               chains  = 4, #chainの回数を指定
               iter    = 1000,  #繰り返しの回数を指定
               warmup  = 500 #ウォームアップの回数を指定
               )
## Compiling the C++ model
## Start sampling
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
## 
## Gradient evaluation took 0.000672 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 6.72 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 50.5034 seconds (Warm-up)
##                155.523 seconds (Sampling)
##                206.026 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
## 
## Gradient evaluation took 0.000426 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 4.26 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 36.8659 seconds (Warm-up)
##                16.2036 seconds (Sampling)
##                53.0696 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
## 
## Gradient evaluation took 0.00039 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 3.9 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 60.3281 seconds (Warm-up)
##                139.548 seconds (Sampling)
##                199.876 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
## 
## Gradient evaluation took 0.000356 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 3.56 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 39.3246 seconds (Warm-up)
##                12.0175 seconds (Sampling)
##                51.3421 seconds (Total)
## Warning: There were 511 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
# 結果の出力
summary(model2b)
##  Family: gaussian(identity) 
## Formula: score ~ 1 + time + (1 + time | id) 
##    Data: longdata (Number of observations: 3000) 
## Samples: 4 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 2000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Group-Level Effects: 
## ~id (Number of levels: 500) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)          18.55      0.67    17.25    19.92        803    1
## sd(time)                1.64      0.10     1.43     1.84        994    1
## cor(Intercept,time)    -0.76      0.03    -0.82    -0.69       1566    1
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept   193.17      0.86   191.53   194.92       1104    1
## time          3.38      0.10     3.18     3.57       1764    1
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     6.64      0.11     6.43     6.84       1142    1
## 
## 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).

パラメータの分布、パラメータが収束しているか、パラメータおよびパラメータの標準誤差などは以下のコマンドを使って確認することができます。

# パラメータの分布を確認
plot(model2b)