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

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

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

2017/9/20