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_{i j} = β_{0} + β_{1} X_{i j} + ε_{i j}, ε_{i j} \sim N (0, σ^{2})

$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_{i j} = β_{0 j} + β_{1 j} X_{i j} + r_{i j} β_{0 j} = γ_{0} + u_{0 j}, β_{1 j} = γ_{1} + u_{1 j},

$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{matrix} u_{0 j} \\ u_{1 j} \end{matrix}) \sim ((\begin{matrix} 0 \\ 0 \end{matrix}), (\begin{matrix} σ_{u_{0}}^{2} \\ σ_{u_{1}, u_{0}} & σ_{u_{1}}^{2} \end{matrix})) \end{aligned} r_{i j} \sim N (0, σ_{r}^{2})

$\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

$j$ によって異なること、および、変数

X_{i j}

$X_{ij}$ の傾きがグループ

j

$j$ によって異なることを認めたモデルとなります。ここで切片および傾きのばらつきは2変量正規分布にしたがうと仮定しておきます（この仮定は後から緩めることもできます）。この場合、新たに投入された

u_{0 j}, u_{1 j}

$u_{0j}, u_{1j}$ のことをランダム効果（Random effect）、前者をとくにランダム切片（Random intercept）、後者をとくにランダム傾き（Random slope）といいます。ここではランダム効果の間に相関を認めていますが、これを0に設定することも可能です。

上の式を1つの式にまとめると、次のようになります（分布についての情報は省略）。

Y_{i j} = γ_{0} + γ_{1} X_{i j} + u_{0 j} + u_{1 j} X_{i j} + r_{i j}

$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_{i j} = γ_{0} + u_{0 j} + r_{i j}, u_{0 j} \sim N (0, σ_{u}^{2}) r_{i j} \sim N (0, σ_{r}^{2})

$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以下の値を取ります。これは以下のように定義されます。

I C C = \frac{σ_{r}^{2}}{σ_{u}^{2} + σ_{r}^{2}}

$\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_{i j} = γ_{0} + γ_{1} X_{i j} + u_{j} + r_{i j}, u_{j} \sim N (0, σ_{u}^{2}) r_{i j} \sim N (0, σ_{r}^{2})

$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_{i j} = β_{0} + β_{1} X_{i j} + β_{2} X_{i j} + u_{0 j} + u_{1 j} X_{i j} + r_{i j}, \begin{aligned} (\begin{matrix} u_{0 j} \\ u_{1 j} \end{matrix}) \sim ((\begin{matrix} 0 \\ 0 \end{matrix}), (\begin{matrix} σ_{u_{0}}^{2} \\ σ_{u_{1}, u_{0}} & σ_{u_{1}}^{2} \end{matrix})) \end{aligned} r_{i j} \sim N (0, σ_{r}^{2})

$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)

# パラメータとパラメータの標準誤差をプロット
stanplot(model2b)

5.3.2 Estimation with Prior 事前分布を指定して推定

先ほどはパラメータの事前分布として一様分布を指定しましたが、事前分布を指定することももちろん可能です。事前分布については以下のような資料が参考になります（といいつつ自分はほとんど読んでいませんが……）。

Gelman, Andrew. 2006. “Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).” Bayesian analysis 1(3) : 515-534.

*https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations

# brm関数を使って推定
model2b <- brm(score ~ 1 + time + (1 + time | id), 
               data    = longdata,
               prior   = c(set_prior("normal(0,10)", class = "b")), # 係数に正規分布を指定。変数の係数ごとにべつべつに別々に指定することもできる。
               chains  = 4, 
               iter    = 1000,
               warmup  = 500
               )

## Compiling the C++ model

## Start sampling

## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).
## 
## Gradient evaluation took 0.000674 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 6.74 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: 52.2943 seconds (Warm-up)
##                132.228 seconds (Sampling)
##                184.523 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
## 
## Gradient evaluation took 0.000358 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 3.58 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: 41.5998 seconds (Warm-up)
##                13.08 seconds (Sampling)
##                54.6798 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
## 
## Gradient evaluation took 0.000364 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 3.64 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: 62.1324 seconds (Warm-up)
##                120.081 seconds (Sampling)
##                182.213 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
## 
## Gradient evaluation took 0.000372 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 3.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: 32.6484 seconds (Warm-up)
##                138.523 seconds (Sampling)
##                171.172 seconds (Total)

## Warning: There were 530 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.54      0.67    17.30    19.86        802    1
## sd(time)                1.64      0.11     1.42     1.84        648    1
## cor(Intercept,time)    -0.76      0.03    -0.82    -0.69       1190    1
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept   193.23      0.85   191.54   194.88        622 1.02
## time          3.38      0.10     3.19     3.58       1500 1.01
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     6.64       0.1     6.44     6.85       1005    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).

6 Three-level Multi-level Model 3レベル・マルチレベルモデル

さらに発展的な分析として、もう1つレベルを追加してみます。このデータは、生徒が複数回観察されているだけでなく、学校で階層化されているという構造を持っています。つまり、学校によって生徒のスコアの平均的な高さが違っていたり、あるいは伸び方が違っていたりする可能性を考えることもできます。先に推定したModel 2はさらに以下のように拡張することができます。

Y_{i j k} = β_{0 j k} + β_{1 j k} X_{i j k} + r_{i j k} β_{0 j k} = γ_{00 k} + u_{0 j k}, β_{1 j k} = γ_{10 k} + u_{1 j k}, γ_{00 k} = δ_{000} + e_{00 k}, γ_{01 k} = δ_{010} + e_{01 k}

$Y_{ijk} = \beta_{0jk} + \beta_{1jk}X_{ijk} + r_{ijk} \\ \beta_{0jk} = \gamma_{00k} + u_{0jk}, \\ \beta_{1jk} = \gamma_{10k} + u_{1jk}, \\ \gamma_{00k} = \delta_{000} + e_{00k}, \\ \gamma_{01k} = \delta_{010} + e_{01k} \\$

\begin{aligned} (\begin{matrix} u_{0 j k} \\ u_{1 j k} \end{matrix}) \sim N ((\begin{matrix} 0 \\ 0 \end{matrix}), (\begin{matrix} σ_{u_{0}}^{2} \\ σ_{u_{1}, u_{0}} & σ_{u_{1}}^{2} \end{matrix})) \end{aligned} \begin{aligned} (\begin{matrix} e_{00 k} \\ e_{10 k} \end{matrix}) \sim N ((\begin{matrix} 0 \\ 0 \end{matrix}), (\begin{matrix} σ_{e_{0}}^{2} \\ σ_{e_{1}, e_{0}} & σ_{e_{1}}^{2} \end{matrix})) \end{aligned} r_{i j k} \sim N (0, σ_{r}^{2})

$\begin{aligned} \begin{pmatrix} u_{0jk} \\ u_{1jk} \end{pmatrix} \sim N \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} \\ \begin{aligned} \begin{pmatrix} e_{00k} \\ e_{10k} \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_{e_0}^2 & \\ \sigma_{e_1, e_0} & \sigma_{e_1}^2 \end{pmatrix} \right) \end{aligned} \\ r_{ijk} \sim N(0, \sigma_{r}^2)\\$ まとめると以下のような式になります（誤差項の分布については省略）。

Y_{i j k} = δ_{000} + δ_{010} X_{i j k} + e_{00 k} + u_{0 j k} + e_{01 k} X_{i j k} + u_{1 j k} X_{i j k} + r_{i j k}

$Y_{ijk} = \delta_{000} + \delta_{010}X_{ijk} + e_{00k} + u_{0jk} + e_{01k}X_{ijk} + u_{1jk}X_{ijk} + r_{ijk}$

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

まず2レベルのときと同様にして、独立変数を投入しない切片のみのモデルを推定してみます。このモデルは以下の通りとなります。

Y_{i j k} = δ_{000} + e_{00 k} + u_{0 j k} + r_{i j k}, r_{i j k} \sim N (0, σ_{r}^{2}) u_{0 j k} \sim N (0, σ_{u}^{2}) e_{00 k} \sim N (0, σ_{e}^{2})

$Y_{ijk} = \delta_{000} + e_{00k} + u_{0jk} + r_{ijk}, \\ r_{ijk} \sim N(0, \sigma_r^2) \\ u_{0jk} \sim N(0, \sigma_u^2) \\ e_{00k} \sim N(0, \sigma_e^2)$

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

## 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) + (1 | school)
##    Data: longdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  23419.6  23443.6 -11705.8  23411.6     2996 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4760 -0.5109  0.0688  0.5921  3.1868 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 172.46   13.132  
##  school   (Intercept)  42.72    6.536  
##  Residual              93.43    9.666  
## Number of obs: 3000, groups:  id, 500; school, 29
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  204.859      1.435   142.7

これを見てみると、個人レベル（レベル2）だけでなく、学校レベル（レベル3）の分散もかなり大きいことがわかります。このことは、学校レベルをモデルに組み込むことの有用性を示しているといえます。

ついで独立変数を加えたモデルはこのようになります。

Y_{i j k} = δ_{000} + δ_{010} X_{i j k} + e_{00 k} + u_{0 j k} + r_{i j k}, r_{i j k} \sim N (0, σ_{r}^{2}) u_{0 j k} \sim N (0, σ_{u}^{2}) e_{00 k} \sim N (0, σ_{e}^{2})

$Y_{ijk} = \delta_{000} + \delta_{010}X_{ijk} + e_{00k} + u_{0jk} + r_{ijk}, \\ r_{ijk} \sim N(0, \sigma_r^2) \\ u_{0jk} \sim N(0, \sigma_u^2) \\ e_{00k} \sim N(0, \sigma_e^2)$

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

## Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
##   approximations to degrees of freedom [lmerMod]
## Formula: score ~ 1 + time + (1 | id) + (1 | school)
##    Data: longdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  22021.9  22051.9 -11006.0  22011.9     2995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9488 -0.5108  0.0418  0.5673  5.0397 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 179.13   13.384  
##  school   (Intercept)  42.72    6.536  
##  Residual              53.38    7.306  
## Number of obs: 3000, groups:  id, 500; school, 29
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  193.0185     1.4613   26.2000  132.09   <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.187

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

3レベルの場合はより複雑になってきますが、以下の3つのパターンを区別してそれぞれ推定して、適合度を比較してみましょう。

生徒レベルの傾きのランダム効果を認めるが、学校レベルの傾きのランダム効果は認めないモデル $Y_{i j k} = δ_{000} + δ_{010} X_{i j k} + e_{00 k} + u_{0 j k} + e_{01 k} X_{i j k} + r_{i j k}$ $Y_{ijk} = \delta_{000} + \delta_{010}X_{ijk} + e_{00k} + u_{0jk} + e_{01k}X_{ijk} + r_{ijk}$
学校レベルの傾きのランダム効果を認めるが、生徒レベルの傾きのランダム効果は認めないモデル
$Y_{i j k} = δ_{000} + δ_{010} X_{i j k} + e_{00 k} + u_{0 j k} + u_{1 j k} X_{i j k} + r_{i j k}$ $Y_{ijk} = \delta_{000} + \delta_{010}X_{ijk} + e_{00k} + u_{0jk} + u_{1jk}X_{ijk} + r_{ijk}$
学校レベルの傾きのランダム効果も、生徒レベルの傾きのランダム効果もどちらも認めるモデル

Y_{i j k} = δ_{000} + δ_{010} X_{i j k} + e_{00 k} + u_{0 j k} + e_{01 k} X_{i j k} + u_{1 j k} X_{i j k} + r_{i j k}

$Y_{ijk} = \delta_{000} + \delta_{010}X_{ijk} + e_{00k} + u_{0jk} + e_{01k}X_{ijk} + u_{1jk}X_{ijk} + r_{ijk}$

# Student-level random-slope model
model5 <- lmer(score ~ 1 + time + (1 + time | id) + (1 | school),
                     data      = longdata,
                     REML      = FALSE,
                     na.action = na.omit)
# School-level random-slope model
model6 <- lmer(score ~ 1 + time + (1 | id) + (1 + time | school),
                     data      = longdata,
                     REML      = FALSE,
                     na.action = na.omit)
# Student-level and School-level random-slope model
model7 <- lmer(score ~ 1 + time + (1 + time | id) + (1 + time | school),
                     data      = longdata,
                     REML      = FALSE,
                     na.action = na.omit)
anova(model4, model5, model7)

## Data: longdata
## Models:
## object: score ~ 1 + time + (1 | id) + (1 | school)
## ..1: score ~ 1 + time + (1 + time | id) + (1 | school)
## ..2: score ~ 1 + time + (1 + time | id) + (1 + time | school)
##        Df   AIC   BIC logLik deviance   Chisq Chi Df Pr(>Chisq)    
## object  5 22022 22052 -11006    22012                              
## ..1     7 21833 21875 -10910    21819 192.537      2  < 2.2e-16 ***
## ..2     9 21785 21839 -10883    21767  52.651      2   3.69e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(model4, model6, model7)

## Data: longdata
## Models:
## object: score ~ 1 + time + (1 | id) + (1 | school)
## ..1: score ~ 1 + time + (1 | id) + (1 + time | school)
## ..2: score ~ 1 + time + (1 + time | id) + (1 + time | school)
##        Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## object  5 22022 22052 -11006    22012                             
## ..1     7 21903 21945 -10944    21889 122.86      2  < 2.2e-16 ***
## ..2     9 21785 21839 -10883    21767 122.33      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

適合度を見てみると、学校レベル、生徒レベルの傾きのランダム効果まで含んだモデルが最も適合度がよいようです。最も適合度の良かったModel 7の結果は次のとおりとなります。

summary(model7)

## Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
##   approximations to degrees of freedom [lmerMod]
## Formula: score ~ 1 + time + (1 + time | id) + (1 + time | school)
##    Data: longdata
## 
##      AIC      BIC   logLik deviance df.resid 
##  21784.7  21838.8 -10883.4  21766.7     2991 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2584 -0.4875  0.0356  0.5258  5.3727 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 272.157  16.497        
##           time          1.794   1.339   -0.73
##  school   (Intercept)  89.688   9.470        
##           time          1.124   1.060   -0.86
##  Residual              43.971   6.631        
## Number of obs: 3000, groups:  id, 500; school, 29
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 192.6856     2.0205  25.3750   95.37  < 2e-16 ***
## time          3.4584     0.2287  26.0700   15.12 2.02e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## time -0.823

6.3 Baysian Estimation ベイズ推定

3レベルモデルもやはりベイズで推定することができます。先ほどの最終モデルであるModel 7をベイズ推定で求めてみましょう。

# ベイズモデルを推定
model7b <- brm(score ~ 1 + time + (1 + time | id) + (1 + time | school), 
               data    = longdata,
               prior   = c(set_prior("normal(0,10)", class = "b")), # 係数に正規分布を指定。
               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.000782 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 7.82 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: 147.776 seconds (Warm-up)
##                60.7502 seconds (Sampling)
##                208.526 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).
## 
## Gradient evaluation took 0.000558 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 5.58 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: 122.5 seconds (Warm-up)
##                41.4612 seconds (Sampling)
##                163.961 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).
## 
## Gradient evaluation took 0.000569 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 5.69 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: 188.055 seconds (Warm-up)
##                101.074 seconds (Sampling)
##                289.129 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).
## 
## Gradient evaluation took 0.000523 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 5.23 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: 190.773 seconds (Warm-up)
##                34.2823 seconds (Sampling)
##                225.055 seconds (Total)

## Warning: There were 41 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(model7b)

##  Family: gaussian(identity) 
## Formula: score ~ 1 + time + (1 + time | id) + (1 + time | school) 
##    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)          16.49      0.62    15.33    17.74        962 1.00
## sd(time)                1.33      0.11     1.11     1.54        909 1.01
## cor(Intercept,time)    -0.73      0.04    -0.81    -0.64       1451 1.00
## 
## ~school (Number of levels: 29) 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)           9.90      1.83     6.86    14.03        604 1.01
## sd(time)                1.10      0.20     0.76     1.56       1250 1.00
## cor(Intercept,time)    -0.81      0.10    -0.94    -0.57       1685 1.00
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept   192.71      2.07   188.60   196.68       1497    1
## time          3.46      0.23     3.02     3.91       1738    1
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     6.64       0.1     6.45     6.85       1276    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(model7b)

# パラメータとパラメータの標準誤差をプロット
stanplot(model7b)

7 Appendix: Plotting Results 補足：結果を図示する

先の最尤法での推定で得たランダム効果がどの程度の分散をともなうものであるのか、プロットして確認することができます。

sjp.lmer(model7, 
         type       = "re",
         facet.grid = FALSE,
         sort.est  = "sort.all",
         y.offset   = .4)

同じことは多分ベイズ推定でもできる気がするのですが、コマンドが今のところわかっていないので今後の課題です。

8 Acknowledgement 謝辞

本ページの作成にあたり、2017年ICPSR summer programで受講したMultilevel Model II: Advanced Topicsでの授業内容および演習で使用したコードが大変参考となりました。InstructorのJohn PoeおよびTAのMichael Giordanoに感謝申し上げます。

階層ベイズモデルを実装する：lme4とbrmsパッケージを用いたマルチレベルモデルの基礎

麦山亮太（東京大学大学院人文社会系研究科博士課程）

2017/9/20