stanで階層線形モデル(勉強版)

author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1

 

個人ごとの傾きと切片、集団レベルの傾きと切片を推定する

(上図は階層線形ベイズモデルを用いた集団レベルの切片と傾きの事後分布。下図は事後予測チェックです)

はじめに

今回は「階層線形モデル」をやってみたいと思います。

参考文献はDoing Bayesian Data Analysisの第4版。

chapter.17の、ややこしい手続き(stanの中でデータを標準化する,標準化しない場合の事前分布の値の計算など)を省略したものになります。

また、今回はわかりやすさを重視して、愚直に書き下したモデルコードをご紹介します。

目的

  1. 反復測定データの加工と変数生成がRでできるようになる。

  2. Rとrstanで階層線形モデルによる個人レベル・集団レベルの予測ができるようになる

想定するケース

階層線形モデルは、反復測定データ(e.g. 一人の人が複数回、独立変数と従属変数の組を持つ)に適用されます。

心理学の領域ですと、「経験サンプリング法(Experience Sampling Method)」で取得したデータが代表的でしょうか。

例えばですが、1日に5回、”怒り感情(従属変数)”と”腰痛(独立変数)”に回答したとします。

欠損がなければ、1人あたり5つのデータ組(従属変数と独立変数)が得られることになります。

今回はこの、腰痛(変数X)と怒り感情(変数Y)を頭のなかにイメージしながら、話を進めていくことにしましょう。

データは1日の調査で、5回調査に回答してもらった場合を想定します。参加者は5人とします。

<データセットの例、全部適当です>

RawDat.jpg

<以下略>

各個人に複数回の回答データがあることがわかります。

分析の目的

今回の研究疑問は「腰痛が強いときほど、怒り感情は高まるか?」というところに定めましょう。

要するに、「強く痛みを感じているときほど不機嫌なのか?」という話です。

考えていく順序としては、

  1. Aさんのデータだけで回帰分析をする(痛みの強さをX軸、怒り感情をY軸)。BさんCさん…Nさんについても同様に回帰分析する。

  2. 各参加者でバラバラな切片と傾きを、集団レベルの階層パラメータでまとめる(全体でひとつの切片と傾きを求める)

というイメージでしょうか。これを同時に行っていくモデルです。

そういう分析ができれば、「腰痛が強いときは、怒り感情がつよいか?」について、個人レベルと集団レベルで、腰痛と怒り感情の関係を推定することができます。

変数の変換と生成

実際の経験サンプリング研究では、いくつかの統制変数を追加的に作成したり、データ自体を加工(センタリングや標準化)するケースがあります。

変数Xを個人内偏差得点に変換する。

いわゆるwith-inセンタリングです。「各回答の腰痛得点 - その個人の腰痛平均値」を計算して、偏差得点をつくります。

これにより、「普段の自分の平均よりも腰痛がひどい時は・・・」について、考えていくことができます。

一次点前のデータXを統制変数に使う

例えば、「ある時点の腰痛→同時点の怒り」だけじゃなくて、1時点前の痛みが次時点の怒りを引き起こしているんじゃないの?という可能性を考えてみましょう。

そういった前時点からの影響を統制するためには、一次点前の腰痛得点を「1時点後(現時点)にずらして独立変数に投入(統制)する」ことが行われます。

今回はそれをLag変数とよびましょう。

独立変数の個人平均を統制変数に使う

「腰痛の強さが強い時ほど、怒り感情は高いか」を考えるときに、「そもそも痛みが平均的に強いひとほど、怒り感情も強いんじゃないの?」

という可能性も考えることができます。

腰痛がもともとつよいというだけで、怒り感情の高さが説明できる部分を統制してあげないと、

「腰痛が強い時ほど」の影響を抽出することはできません。

「個人の腰痛平均 - 全体の腰痛平均」を計算して(いわゆるbetweenセンタリング)、

「その人の平均的な腰痛レベルが集団の中で高いか低いか」の相対的な値を生成し、統制変数に入れることで、

平均的な痛みが怒り感情を説明する部分を統制することができます。

こういったデータの変換も、Rで行っていきましょう。tidyverseという強力な味方がありますので、ちょっと使ってみたいと思います。

出来上がり変数のイメージ(適当です。大体こんなかんじ)

CompDat.jpg

なお、どの変数をどういう風にセンタリングするかは研究の目的次第ですが、センタリングまたは正規化(scaling)は行ったほうが、推定がうまくいきます。

標準化しない傾きや切片を求めたい場合は、標準化して推定したものを元に単位に戻す作業が必要になります。

今回は、単位は変えず、センタリングのみを行って傾きと切片を個人ごと、集団レベルで求めたいと思います。

ライブラリの読み込み


#####################################
# Hierarchical Bayesian Linear Model
#####################################

invisible({rm(list=ls());gc();gc()})

#library and directory
library(rstan)
library(ggplot2)
library(shinystan)
library(reshape2)
library(loo)
library(tidyverse) 
library(plotly)
library(Rmisc)
library(GGally)
library(ellipse)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


#
# 目的: 反復測定データ(一人に複数回の回答があるデータ)に対して
#         個人ごとに切片と傾きを推定する。
#         それと同時に、切片や傾きの集団レベルの階層パラメータを推定する。
#

サンプルデータの作成

サンプルデータは本当に適当に作りました。

#サンプルデータの作成。適当に作りました。
ID <- c(1,1,1,1,2,2,2,2,2,2,
        3,3,3,3,3,3,3,3,3,3,
        4,4,4,4,4,4,4,4,4,4,
        5,5,5,5,5,5,5,5,5,5)

Y <- c(4,6,4,6,3,8,4,3,8,4,
       11,5,8,3,2,11,5,5,2,7,
       9,4,7,5,6,7,1,4,6,5,
       10,3,3,3,2,9,4,2,2,1)

X <- c(2,3,2,3,3,6,1,2,6,2,
       9,6,3,5,4,9,2,2,3,4,
       6,1,9,6,4,5,2,3,6,7,
       13,4,2,2,1,1,4,2,1,3)

longDat<-data.frame(cbind(ID,Y,X))

加工前のデータはこういう形です

print(longDat)

> print(longDat)
   ID  Y  X
1   1  4  2
2   1  6  3
3   1  4  2
4   1  6  3
5   2  3  3
6   2  8  6
7   2  4  1
8   2  3  2
9   2  8  6

<略>

データの変換、加工

続いて、変数の変換を行っていきます。

センタリングやLag変数を、こう、ごにょごにょして作ります。

Excelでやるより1000倍くらい早いはずです。

コメントアウトに作業内容を付しておきます。

なお、今回は途中で欠損を削っています。欠損まで推定したい場合は別の手続きが必要です。


#経験サンプリングで頻繁に使用する変数生成とセンタリング
# 1. 変数Xのwith-inセンタリング : withinXの生成
# 2. 一時点前の変数X得点を生成して、with-in でセンタリング : withinLagX
# 3. 個人ごとの平均得点を全体平均でbetweenセンタリング: betweenMeanX

#以下パッケージtidyverseを用いたデータ処理を行う。

dat<-longDat %>% #データフレームであるlongDatを使う。→ %>%で次の関数に渡す

    group_by(ID) %>%  #IDごとにグループとみなす → %>% でグループ化した状態を次の関数にわたす
    
    dplyr::mutate(meanX = mean(X)) %>% #IDごとにXの平均を算出し、mutate関数で算出したものを列に追加。
    
    dplyr::mutate(lagX=lag(X)) %>% #データを一時点後の行にずらす
    
    dplyr::mutate(withinX = (X-meanX),withinLagX = (lagX-meanX)) %>% #with-inでセンタリング。個人ごとの各データに対して、個人内平均を引く
    
    ungroup() %>% #グループ化を解除する。これで、データ全体の処理
    
    dplyr::mutate(betweenMeanX = (meanX-mean(X))) %>% #Xの個人平均が,全体の平均からどれくらい離れているかbetweenセンタリング
    
    dplyr::mutate(Y = (Y-mean(Y)))%>% #従属変数をbetweenでセンタリングする
    
    dplyr::select(ID,Y,withinX,betweenMeanX,withinLagX) %>% #分析に必要な列だけ取り出す
    
    tidyr::drop_na(withinLagX) #lag変数によって欠損が出た行を削除。なお、欠損を別途推定する場合はこの手続は不要。


#dplyrでデータを扱うと、普通のdata.frameと表示の異なるデータフレームになるので、不安な場合はこうする
compDat<-data.frame(dat) 

変換後の完成データは次のような形です。数行のコードで上手く行きました。

print(compDat)
   ID     Y   withinX betweenMeanX withinLagX
1   1  0.95  0.500000   -1.4750000 -0.5000000
2   1 -1.05 -0.500000   -1.4750000  0.5000000
3   1  0.95  0.500000   -1.4750000 -0.5000000
4   2  2.95  2.666667   -0.6416667 -0.3333333
5   2 -1.05 -2.333333   -0.6416667  2.6666667
6   2 -2.05 -1.333333   -0.6416667 -2.3333333
7   2  2.95  2.666667   -0.6416667 -1.3333333
8   2 -1.05 -1.333333   -0.6416667  2.6666667
9   3 -0.05  1.300000    0.7250000  4.3000000

<以下略>

つづいて、行数や参加者数を取得し、Stanに渡すデータとしてまとめます。


#ID番号の最大値を取って、参加者数とする(もともとfactor型でIDが入っていたりすると間違うので以下の手続きをとった)
Nsub <- max(as.integer(factor(compDat$ID)))

#データの行数を取得
Nrow<-length(compDat$ID)

#Stanに渡すデータの作成
#今回は、変数を見やすいように、全てベクトルで渡すことにする
#線形代数を使った行列演算は行わない

ID <- compDat$ID #ID列を一列で取り出す
Y <- compDat$Y #従属変数Yを一列で取り出す
withinX <- compDat$withinX #withinセンタリングした独立変数X
withinLagX<-compDat$withinLagX #withinセンタリングした一時点前のデータX
betweenMeanX<-compDat$betweenMeanX #betweenセンタリングした個人平均X

# stanに渡すlistデータ
datastan <- list(ID=ID,Nsub=Nsub,Nrow=Nrow,
                 Y=Y,
                 withinX=withinX,
                 withinLagX=withinLagX,
                 betweenMeanX=betweenMeanX)

stan設定、モデルと推定の実行

今回は、モデルを愚直に書き下します。

行列を使えばもっとスマートになりますが、今回はそれぞれの傾きのパラメータを一つ一つ宣言し、回帰式もきちんと書いていきます。


#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(6000)
#Warm-up数
warmup=c(3000)
#チェーン数
chains=c(3)
#thinning
thin =c(3)


parameters = c("interceptMu","interceptSigma", #個人ごとの切片の集団レベル平均と標準偏差
               "beta_wXmu","beta_wXsigma", #個人ごとwithinXの傾きに対する集団レベル平均と標準偏差              
               "beta_wLagXmu","beta_wLagXsigma",  #個人ごとwithinLagX傾きに対する集団レベル平均と標準偏差  
               "beta_bMeanXmu","beta_bMeanXsigma", #個人ごとbetweenMeanXの傾きに対する集団レベル平均と標準偏差
               "intercept", #個人ごとの切片
               "beta_wX", #wXの個人レベル傾き
               "beta_wLagX", #wLagXの個人レベル傾き
               "beta_bMeanX", #bMeanXの個人レベル傾き
               "sigma", #回帰式からデータまでの誤差(標準偏差)
               "mu", #線形モデルの予測値
               "log_lik" #対数尤度(情報量規準の算出に使用)
)




#StanModel
HBLM <- "
data {
int Nsub;
int Nrow;
int ID[Nrow];
real Y[Nrow];
real withinX[Nrow];
real withinLagX[Nrow];
real betweenMeanX[Nrow];
}

parameters {
//Hierarchical
real interceptMu;
real <lower=0> interceptSigma;
real beta_wXmu ;
real <lower=0> beta_wXsigma ; 
real beta_wLagXmu ;
real <lower=0> beta_wLagXsigma ; 
real beta_bMeanXmu ;
real <lower=0> beta_bMeanXsigma ; 

//regresson error
real <lower=0>sigma;

//Individual
real intercept [Nsub];
real beta_wX [Nsub];
real beta_wLagX [Nsub];
real beta_bMeanX [Nsub];
}

transformed parameters {
real mu [Nrow];
for (i in 1:Nrow){
mu[i] = intercept[ID[i]]+beta_wX[ID[i]]*withinX[i]+beta_wLagX[ID[i]]*withinLagX[i]+beta_bMeanX[ID[i]]*betweenMeanX[i];
}
}

model {
//Hierarchical
interceptMu ~ normal(0,100);
interceptSigma ~ cauchy(0,2.5);
beta_wXmu ~ normal(0,100);
beta_wXsigma ~ cauchy(0,2.5); 
beta_wLagXmu ~ normal(0,100);
beta_wLagXsigma ~ cauchy(0,2.5); 
beta_bMeanXmu ~ normal(0,100);
beta_bMeanXsigma ~ cauchy(0,2.5); 

//error
sigma ~ cauchy(0,2.5);


//Individual
for(i in 1:Nsub){
intercept[i] ~ normal(interceptMu,interceptSigma);
beta_wX[i] ~ normal(beta_wXmu,beta_wXsigma);
beta_wLagX[i] ~ normal(beta_wLagXmu,beta_wXsigma);
beta_bMeanX[i] ~ normal(beta_bMeanXmu,beta_bMeanXsigma);
}

//each data
for (j in 1:Nrow){
Y[j] ~ normal(mu[j],sigma);
}

}

generated quantities {
real log_lik[Nrow];
for (i in 1: Nrow){
log_lik[i] = normal_lpdf(Y[i]|mu[i],sigma);
}
}


"

#stan fit
fit <- stan(model_code = HBLM,
            data=datastan, 
            pars=parameters,
            iter=iterations, 
            chains=chains, 
            thin=thin,
            warmup = warmup)

上記モデルのポイントは2つです。

一つ目は、 mu[i] = ・・・の一番コアな部分。beta_wX[ID[i]]っていう使い方がミソです。

反復測定したデータのi行目にあるデータの予測値をmu[i]とし、そのi番目の予測値を、IDのi番目にある数字(つまり何番目の参加者か)のパラメータを使って定義しています。 こうすれば、for文を繰り返し使ったりしなくてもよくなります。

もうひとつのポイントは、参加者ごとに推定される切片と傾きが、集団レベルで平均値と標準偏差を共有するように組まれている(階層化)されている点です。

例えば、 i人分あるintercept[i]は、interceptMuを平均とし、interceptSigmaを標準偏差とする正規分布に従うことが仮定されています。

これによって、集団レベルの切片の平均と標準偏差を推定することができます。

また、このような階層パラメータを設定することにより、データ数の少なかった(欠損が多かった)参加者についても、

比較的安定した推定を行うことができます。

推定結果の確認

パラメータの上から15個だけ取り出して見てみます。

お、うまくいってそうですね。が、たぶん幾つかのパラメータで収束してないようです。(もう少し、iterationとthinを増やしてみてください。)

階層レベルの傾き(個人ごとの傾きの集団レベル推定平均)はどれも95%を区間に0を含んでいます。 beta_wXmu(腰痛がいつもより強い時は…のデータの集団レベル傾き)がおしいですね。

round(summary(fit)$summary[1:15,c(1,3,4,8,9,10)],digit=3)


                   mean    sd   2.5% 97.5%    n_eff  Rhat
interceptMu      -0.129 1.058 -2.584 1.739  532.354 1.011
interceptSigma    1.025 0.984  0.107 3.625  371.574 1.011
beta_wXmu         0.605 0.363 -0.143 1.320 1713.187 1.002
beta_wXsigma      0.438 0.451  0.027 1.736  281.610 1.019
beta_wLagXmu     -0.047 0.270 -0.535 0.443 1516.114 1.002
beta_wLagXsigma   0.278 0.338  0.009 1.175  259.526 1.013
beta_bMeanXmu     0.371 1.117 -1.860 2.763 1364.006 1.001
beta_bMeanXsigma  1.186 1.215  0.125 4.354  312.690 1.019
intercept[1]      0.036 1.471 -3.112 2.908  794.201 1.007
intercept[2]      0.109 1.081 -2.168 2.265  600.266 1.012
intercept[3]     -0.016 1.070 -2.484 2.174  871.782 1.001
intercept[4]     -0.185 1.163 -2.829 2.122  685.831 1.005
intercept[5]     -0.562 1.151 -3.152 1.380  139.486 1.023
beta_wX[1]        0.671 0.644 -0.530 2.123 2232.345 1.002
beta_wX[2]        0.745 0.327  0.165 1.470 1813.426 1.004

Gelman-Rubin testは、 shinystan起動 → DIAGNOSE →Neffとかのタブで確認してください。

全部 “ None “ ならOKです。

shinystanの起動コマンドはこちら

launch_shinystan(fit)

事後分布のプロット

長くなってきたので説明はつけませんが、事後分布をプロットしていきます。

今回は集団レベル(階層パラメータ)の推定値のみプロットします。


#パラメータのggplot----------------------------------------

#結果を確認したいパラメータ(ここでは、階層レベルの切片と傾き3つ)
parameters2<-c("interceptMu","beta_wXmu","beta_wLagXmu","beta_bMeanXmu")
#結果を確認したいパラメータに対応した名前
varName<-c("Intercept","SlopeWithinX","SlopeWithinLagX","SlopeBetweenMeanX")


#fitに格納したMCMCサンプルを取り出す
d.ext<-rstan::extract(fit) 

#関心のあるパラメータだけ抜き出す
df.ext<-data.frame(d.ext[c(parameters2)])

#名前を上記にしたがってつけ直す
names(df.ext) = varName


#整然データになるよう、meltする
df.comp<-melt(df.ext,idvar=c())


#ggplotの真っ白なキャンバスを用意
g<-ggplot()+theme_set(theme_bw(base_size=18,base_family="HiraKakuProN-W3"))

#使用するデータ、x軸、y軸に何の変数を使うか、色は何の変数で塗り分けるか、塗りつぶしfillは何の変数で行うか指定
g <- ggplot(df.comp,aes(x = value,y = ..density..,colour=variable,fill=variable))

#0の値に黒の縦線を入れる
g <- g+ geom_vline(xintercept = 0)

#ヒストグラムを上書きしていく
g <- g + geom_histogram(position = "identity",alpha = 0.5)

#確率密度曲線を重ね書き
g <- g + geom_density(stat = "density",position = "identity",alpha = 0.4)

#パラメータ(variable)ごとに図を分割
g<-g+facet_wrap(~variable,scales="free")

#線の色のパレットを一括指定
g<-g + scale_color_brewer(palette="Set2")
#塗りつぶし色のパレットを指定
g<-g + scale_fill_brewer(palette="Set2")

#タイトルを各
g<-g+ggtitle("Posteriors: Hierarchical Linear Parameters")


#軸の調整
g <- g + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
g <- g + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
g <- g + theme(legend.title = element_text(size=14), #legend size
               legend.text = element_text(size=10),  #legend text size 
               legend.key= element_rect(colour = "white")) # legend mark-box colour

#描画
g

お、一時点前のLag変数と(WithinLagX)、個人の腰痛平均(BetweenMeanX)を統制したときに、

「腰痛がいつもより強いときは、怒り感情も高まる」というwithinXの効果がでそうでしたね。 

惜しいですが、withinXの集団レベル傾きの95%区間には、0が含まれています。

MAP推定値

一応MAP推定値も出しておきます。


#MAP estimate-----------------------------------------------------------------

#MAP estimate

MAPs  = c(rep(NA,length(df.ext)))
for (i in 1:length(df.ext)){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
names(MAPs) <- c(colnames(df.ext))
#-----------------------------------------------------------------------------

結果です。

print(MAPs)

        Intercept      SlopeWithinX   SlopeWithinLagX SlopeBetweenMeanX 
       -0.3675120         0.6800010        -0.1042150         0.4731112 

視覚的事後予測チェック

説明無し。一気に行きます。

#Posterior Predictive Checks---------------------------------------------------
#いくつの事後予測分布を生成するか?
repN<-9

#今回は反復測定したデータだったので、データの行数分乱数を生成する
Nrow<-Nrow

#事後分布サンプルからの乱数生成##########################################
predY<-matrix(NA,Nrow,repN)
for(i in 1:repN){
    
    #modelにあわせて変更--------------------------------------
    predY[,i] <-rnorm(Nrow,d.ext$mu[i,],d.ext$sigma)
    #---------------------------------------------------------

    }

#########################################################################



#整然データになるようmeltする
m.Dat<-melt(predY,idvar=c())
#列名を付け替える
colnames(m.Dat)<-c("n", #何番目の乱数データか
                   "sample", #何番目の事後予測分布か
                   "pred_y") #予測値となる事後分布からの乱数の値

#事後予測分布の番号をfactor型にしておく
m.Dat$sample<-as.factor(m.Dat$sample) 



#ggplot------------------------------------------------------------------

#データYを使って、x軸にはYの値を指定、その後、ヒストグラムを書く。縦軸は確率密度。alphaは透過性。binwidthは一つの棒グラフに入れるYの値の区間
p<-ggplot(data.frame(Y),aes(x=Y))+ geom_histogram(aes(y=..density..),alpha=0.5,binwidth=1)

#theme and graph title size
p<-p+theme_bw(base_size=18,base_family="HiraKakuProN-W3") # Graph Title Size and Font
p<-p+ggtitle("Posterior Predictive Checks")

#事後予測分布の確率密度曲線を重ね書き。何番目の事後予測分布なのか(sample)によって色を分ける
p<-p+geom_line(data=m.Dat,aes(x=pred_y,y=..density..,colour=sample,group=sample),stat="density")

#軸の調整
p <- p + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
p <- p + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
p <- p + theme(legend.title = element_text(size=10), #legend size
               legend.text = element_text(size=8),  #legend text size 
               legend.key= element_rect(colour = "white")) # legend mark-box colour
#Colour Variations
#Sequential(連続カラー): Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
#Discrete(離散カラー)  : BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
p <- p + scale_color_brewer(palette="OrRd")

#描画
p

おりゃっ!

データが少ないにも関わらず、割りと予測がなせている方では無いでしょうか。

それでは、

 EnjoyStan!!

Written on January 13, 2017