この記事は Stan Advent Calendar 2016 の17日目の記事です。目下勉強しながらの執筆のため，もし間違いやもっと効率的な方法があれば，ご指摘をいただければ幸いです。

#はじめに
媒介分析については，こちらのサイトをご参照ください。このサイトでは，野球選手の体重（weight）が重くなるほど，ホームラン数（HR）が増えた結果，年収（salary）が増えるという，媒介過程を例に説明を行っています。Stanコードも公開されています。
なお本記事では、以下のように変数を略記することにします。

X：説明変数
M：媒介変数
Y：応答変数

#階層性がある媒介モデル（マルチレベル媒介モデル）
上記の野球のデータでは，各変数（体重，ホームラン数，年収）の値は，個人につき一つでした。しかし，これらのデータを複数回測定すれば，ある変数について個人の中に複数の値が存在することになり，階層性が生まれます（例：体重の計測1回目70kg，2回目71kg，3回目69kg...）。

あるいは，球団と選手の関係を考えてみると，選手はどこかの球団に所属しているので，やはり階層的なデータとみなすことが出来ます（例：同一球団内の選手A，選手B，選手C...）

通常の回帰分析は，残差が独立であることを仮定しますが，データに階層性が存在する場合は，その前提が崩れる恐れがあります。そこで，マルチレベルモデルや階層線形モデル，線形混合モデルと呼ばれる方法により，階層の上位ごと，下位ごとの違いもモデリングします。選手のデータを数回計測する場合を例に挙げると，計測間での違いと選手間の違いを，それぞれモデリングすることになります。そのうえで，以下二つの効果の積として間接効果を導出します。

XがMに及ぼす影響
MがYに及ぼす影響

#実践
実は，Stanでマルチレベル媒介モデルによる分析を行うための，bmlmというRパッケージが存在します。もうこれを使えばいいですね。
　　
...では話が終わってしまうので，このパッケージからサンプルデータを拝借して，自分でStanコードを書いてみます。使用するのは，BLch9というデータセットです。100名に対して経時的なデータを21回にわたり測定した、階層性のあるデータです。個人を識別する変数'id'に加え，'fwkstrs'，'fwkdis'，'freldis'，という三つの変数を使用します。変数名が似ているので，以下のように改名しておきます。

X（元はfwkstrs）：仕事におけるストレッサーの数
M（元はfwkdis）：仕事への不満の大きさ
Y（元はfreldis）：人間関係への不満の大きさ

また，元データでは個人IDは101~200の100名分存在しますが，分かりやすさのため，1~100に変えてしまいます。

sample_data

library(bmlm)
library(dplyr)

data(BLch9)
dat <- BLch9 %>%
  dplyr::select(id, fwkstrs, fwkdis, freldis) %>% 
  dplyr::rename(X=fwkstrs, M=fwkdis, Y=freldis)
dat$id <- dat$id - 100

##データ可視化
まずデータの特徴を確認します。ひとまず個人差を考慮せず，2100行のデータを用いて散布図を描いてみます。

visualization2

library(GGally)

dat %>% 
  dplyr::select(-id) %>% 
  GGally::ggpairs()

以降の分析では，各変数は正規分布するものと仮定することにします。
次に，XがMに及ぼす影響，MがYに及ぼす影響には，個人差が存在するかどうかを視覚化します。100名分のデータを描画すると煩雑になるので，ID1~10の10名を抜粋します。

visualization2

library(ggplot2)
library(cowplot)

temp <- dat %>% dplyr::filter(id<=10)
temp$id <- as.factor(temp$id)

g1<- ggplot(data=temp, aes(y=M, x=X, color=id))+
  stat_smooth(method="lm", se=FALSE)
g2<- ggplot(data=temp, aes(y=Y, x=M, color=id))+
  stat_smooth(method="lm", se=FALSE)

g <- cowplot::plot_grid(g1, g2, labels=c("X->M model", "M->Y model"), align="h") 
print(g)

切片も傾きも，個人ごとにばらばらです。そこで，切片と傾きに個人差を仮定したモデリングを行います。

##変量効果間に相関を仮定”しない”場合のStanコード
タイトルはひとまず無視して読み進めてください。
本記事では，StanとRでベイズ統計モデリングの第8章（階層モデル）を参考にさせていただきました。この章では，単回帰分析において，切片と傾きに個人差を仮定した場合の分析例が紹介されています（'model8-4.stan'）。

multilevel_mediation1

data {
  int N;                               
  int K;                            //グループ数
  real X[N];                        //説明変数
  real Y[N];                        //応答変数
  int<lower=1, upper=K> KID[N];     //グループを識別するID
}

parameters {
  real a0;                          //切片の全体平均
  real b0;                          //Xの係数の全体平均
  real a[K];                        //グループごとの切片
  real b[K];                        //グループごとの傾き
  real<lower=0> s_a;                //切片のグループ間でのばらつき
  real<lower=0> s_b;                //Xの傾きのグループ間でのばらつき
  real<lower=0> s_Y;                //応答変数のばらつき
}

model {
  for (k in 1:K) {
    a[k] ~ normal(a0, s_a);         //切片のグループ差は，この正規分布から発生
    b[k] ~ normal(b0, s_b);         //Xの傾きのグループ差は，この正規分布から発生
  }

  for (n in 1:N)                    //応答変数に関する確率モデル
    Y[n] ~ normal(a[KID[n]] + b[KID[n]]*X[n], s_Y);
}

このコードをマルチレベル媒介モデルに適用するためには，以下2点の追記が必要になります。

Xが複数存在する，重回帰分析のコードを加筆
回帰係数を掛け合わせて，間接効果を求めるコードを加筆

まず1点目について。
媒介分析は，X->Yの直接的な関係を，X->M，M->Yという間接的な関係として表現できるかどうかを知るための分析でした。ここでX->Mの分析は，まさに単回帰分析に相当します。したがって，上記のStanコードをそのまま拝借することが出来ます。ただし，Xで予測したいのはYではなくMですから，その点を書き換える必要があります。また，このコード内ではX->Mモデルにおけるパラメータを推定したいことが分かるように，パラメータ名の最後に'__M'と加筆しておきます。

M->Yモデルも一見すると単回帰分析に見えますが，Mの背後にはXが存在するので，XとMでYを予測する，重回帰分析になります。とはいえ，説明変数の数だけ，推定する傾きパラメータ（及びそのために必要なパラメータ）を増やせば済みます。M->Yモデルにおけるパラメータを推定したいことが分かるように，パラメータ名の最後に'__Y'と加筆しておきます。

次に2点目について。
こちらのサイトで紹介されているように，generated quantitiesブロックにおいて，以下の二つを掛け合わせるだけです。

X->MモデルにおけるXの傾き
M->YモデルにおけるMの傾き

これらをStanコードで表したものが以下です。

multilevel_mediation1

data{
  int N;                                 //総データ数
  int K;                                 //グループ数
  vector[N] X;                           //説明変数
  vector[N] M;                           //媒介変数
  vector[N] Y;                           //応答変数
  int<lower=1, upper=K> KID[N];          //グループを識別するID
}     
      
parameters {
  // X->M model
  real a0__M;                            //切片の全体平均
  real b0__M;                            //Xの係数の全体平均
  real a__M[K];                          //グループごとの切片 
  real b__M[K];                          //グループごとのXの傾き
  real<lower=0> s_a__M;                  //切片のグループ間でのばらつき
  real<lower=0> s_b__M;                  //Xの傾きのグループ間でのばらつき
  real<lower=0> s_Y__M;                  //媒介変数のばらつき
  
  // M->Y model      
  real a0__Y;                            //切片の全体平均
  real b0__Y1;                           //Xの係数の全体平均
  real b0__Y2;                           //Mの係数の全体平均
  real a__Y[K];                          //グループごとの切片
  real b__Y1[K];                         //グループごとのXの係数
  real b__Y2[K];                         //グループごとのMの係数
  real<lower=0> s_a__Y;                  //切片のグループ間でのばらつき
  real<lower=0> s_b__Y1;                 //Xの係数のグループ間でのばらつき
  real<lower=0> s_b__Y2;                 //Mの係数のグループ間でのばらつき
  real<lower=0> s_Y__Y;                  //応答変数のばらつき
}

model {
  for (k in 1:K) {
    // X->M model
    a__M[k] ~ normal(a0__M, s_a__M);     //切片のグループ差は，この正規分布から発生
    b__M[k] ~ normal(b0__M, s_b__M);     //Xの傾きのグループ差は，この正規分布から発生
    
    // M->Y model
    a__Y[k] ~ normal(a0__Y, s_a__Y);     //切片のグループ差は，この正規分布から発生
    b__Y1[k] ~ normal(b0__Y1, s_b__Y1);  //Xの傾きのグループ差は，この正規分布から発生
    b__Y2[k] ~ normal(b0__Y2, s_b__Y2);  //Mの傾きのグループ差は，この正規分布から発生
  }

  for (n in 1:N) {                       //MとYに関する確率モデル
    M[n] ~ normal(a__M[KID[n]] + b__M[KID[n]]*X[n], s_Y__M);
    Y[n] ~ normal(a__Y[KID[n]] + b__Y1[KID[n]]*X[n] + b__Y2[KID[n]]*M[n], s_Y__Y);
  }
  
  //prior
  s_a__M ~ cauchy(0, 2.5);               //推定の安定性のため，ばらつきの事前分布に半コーシー分布を指定
  s_b__M ~ cauchy(0, 2.5);
  s_a__Y ~ cauchy(0, 2.5);
  s_b__Y1 ~ cauchy(0, 2.5);
  s_b__Y2 ~ cauchy(0, 2.5);
}

generated quantities{
  real indirect;                           
  indirect = b0__M * b0__Y2;             //傾きの積により間接効果を計算
}

上記のStanコードを適当な名前で保存し，以下のRコードでサンプリングを試してみます。ここでは'BMM.stan'という名前で保存しました。

R_code

library(rstan)

stanmodel <- stan_model('BMM.stan')
data=list(N=nrow(dat), K=max(dat$id), Y=dat$Y, M=dat$M, X=dat$X, KID=dat$id)
fit <- sampling(stanmodel,
                data=data,
                iter=20000,
                warmup=3000,
                thin=5,
                seed=1234)

print(fit, pars=c("a0__M","b0__M","s_a__M","s_b__M","s_Y__M","a0__Y","b0__Y1","b0__Y2","s_a__Y","s_b__Y1","s_b__Y2","s_Y__Y","indirect"))

結果は以下の通りです。なおこれらは標準化された値ではありません。

今回は，iteration=20000，warmup=3000という大きな回数でサンプリングしているので，Rhatはどれも許容範囲に収まっていますが，'s_b__Y1'というパラメータのn_effが他と比べて小さいことが分かります。トレースプロット等のグラフは省略しますが，確認してみると自己相関が高そうでした。事前にサンプリングしてみて，こうなることが分かっていたので，thin=5でthinningを行ったのですが，あまり改善されていないみたいです。thinをもっと大きくするか，StanとRでベイズ統計モデリングの第4章にあるように，そもそものモデルを改良する必要があるかもしれません。とりあえず今回はこのまま進めます。

今回一番関心があるのは，間接効果（indirect）です。以下のグラフを見る限り，収束には問題なさそうです。分布は左右に大きく歪んでいるわけではないので，事後平均0.03を代表値として採用します。95%確信区間は0をまたいでいないので，「仕事におけるストレッサーが増えるほど，仕事への不満が大きくなり，その結果として人間関係への不満も大きくなる」という媒介過程が存在すると，確信を持てそうです。

##変量効果間に相関を仮定”する”場合のStanコード
ここでもう一度，以下のグラフに注目してみます。特に左のX->Mモデルが分かりやすいと思うのですが，何となく，切片が小さいほど傾きが大きい傾向があるように思われます。

もともと仕事に対する不満が低ければ，少しのストレッサーの増加でも，不満は上昇しやすいと考えられます。これはこのサンプルデータに限った話ではなく，ホームラン数と年収の関係にも同じように，切片が小さい個人ほど傾きが大きくなりやすいという，負の相関が認められると思われます。伸びしろですね。

というわけで，Stanコードを書き換えてみます。切片と傾きの個人差は，それぞれ独立に正規分布から発生したと仮定していましたが，さらにこれらの間に相関関係を想定します。そこでこちらのスライドの140枚目以降に記載されているように，多変量正規分布を仮定した書き方に変更します。

multilevel_mediation2

data{
  int N;                                 //総データ数
  int K;                                 //グループ数
  vector[N] X;                           //説明変数
  vector[N] M;                           //媒介変数
  vector[N] Y;                           //応答変数
  int<lower=1, upper=K> KID[N];          //グループを識別するID
}     

parameters {
  vector[5] f;                           //固定効果（平均）
  vector[5] r[K];                        //変量効果
  vector<lower=0>[5] tau;                //変量効果のSD
  corr_matrix[5] rho;                    //変量効果の相関行列
  real<lower=0> sigma__M;                //Mの残差SD
  real<lower=0> sigma__Y;                //Yの残差SD
}

transformed parameters{                  //共分散行列は指定が難しいためパラメータとして宣言せず，SDと相関行列に分けてこれらをパラメータで宣言する
  cov_matrix[5] taumx;                   //ランダムパラメータの共分散行列
  taumx = quad_form_diag(rho, tau);      //相関行列と（変量効果の）SDを共分散行列に変換
}

model {
  //MとYに関する確率モデル
  for (n in 1:N) {                       
    M[n] ~ normal(r[KID[n]][1] + r[KID[n]][2]*X[n], sigma__M);
    Y[n] ~ normal(r[KID[n]][3] + r[KID[n]][4]*X[n] + r[KID[n]][5]*M[n], sigma__Y);
  }
  
  //ランダムパラメータの分布　多変量正規分布
  r ~ multi_normal(f, taumx);

  //事前分布
  tau ~ cauchy(0, 2.5);
  rho ~ lkj_corr(1); 　　　　　　　　　　　//相関行列の無情報事前分布
  sigma__M ~ cauchy(0, 2.5);
  sigma__Y ~ cauchy(0, 2.5);
}

generated quantities{
  real indirect;                           
  indirect = f[2]*f[5];  　　　　　　　　 //傾きの積により間接効果を計算
}

上記のStanコードを'BMM2.stan'という名前で保存し，以下のRコードを走らせます。

R_code2

library(rstan)

stanmodel_2 <- stan_model('BMM2.stan')
data=list(N=nrow(dat), K=max(dat$id), Y=dat$Y, M=dat$M, X=dat$X, KID=dat$id)
fit2 <- sampling(stanmodel_2,
                data=data,
                seed=1234)

print(fit2, pars=c("f", "tau", "indirect"), digits=3)

サンプリング回数はデフォルトの設定にしましたが，収束はそんなに悪くなさそうです。間接効果の事後平均は，変量効果間に相関を仮定してもしなくても，あまり変わりませんでした。データ次第では結果が変わるかもしれませんが，このあたりは未検証です。

#検算
##bmlmパッケージ
間接効果の推定がうまくいっているか，bmlmパッケージを使って計算した結果と比較してみます。Stanで書いた場合と設定を統一してサンプリングしてみます。

bmlm_package

library(bmlm)

fit3 <- mlm(dat,
            id="id", x="X", m="M", y="Y",
            priors=list(tau_dy=2.5, tau_dm=2.5, tau_a=2.5, tau_b=2.5),
            iter=20000,
            warmup=3000,
            thin=5,
            seed=1234)

mlm_summary(fit3)
mlm_pars_plot(fit3, type="hist")

結果は以下の通りです。パラメータaがX->Mモデルにおける回帰係数，パラメータbがM->Yモデルにおける回帰係数を表しています。これらの事後平均の積を求めると，a * b = 0.19 * 0.15 = 0.0285なので，Stanで書いた場合の間接効果とほぼ一致します。ただ，パラメータabが間接効果を表しているのですが，事後平均は0.05になっていますね...。

bmlmでは，回帰係数の積に回帰係数間の共分散（covab）を加算した値が，間接効果（ab）として出力されます。ここでは小数点以下2桁までしか表示されていないので，丸められてしまって分かりづらいですが，abからcovabを減算すると，0.05 - 0.03 = 0.02ということで，上記の通りStanで書いた場合とおおよそ一致します。

##Rでシミュレーション
StanとRでベイズ統計モデリングの第8章（階層モデル）に，Rでシミュレーションするのは大切なステップであると明記されています。2変量正規分布に従う乱数を発生させ，以下二通りの方法で結果を比較してみます。

各変数の積を平均化する場合
各変数の平均の積を求める場合

simulation

library(mvtnorm)

sigma <- matrix(c(4,2,2,3), ncol=2)
x <- rmvnorm(n=500, mean=c(1, 2), sigma=sigma) #2変量正規分布の生成
colMeans(x) #各変数の平均
var(x) #共分散行列

product1 <- mean(x[,1] * x[,2]) #各変数の積を平均
print(product1)

product2 <- colMeans(x)[1] * colMeans(x)[2] #各変数の平均の積
print(product2)

print(product1 - product2) #差分

結果は以下の通りです。各変数の平均の積を計算した場合は，各変数の積を平均した場合よりも値が小さいです。そしてその差分は，変数間の共分散にほぼ一致します。

上記のStanコードでも，generated quantitiesブロックにおける間接効果の計算を以下のように変更すれば，bmlmの結果とおおよそ一致します。

indirect

generated quantities{
  real indirect;                           
  indirect = f[2]*f[5] + taumx[2, 5];  //傾きの積により間接効果を計算
}

というわけで，どうやらマルチレベル媒介モデルにおいて間接効果を導出する際は，各回帰係数の積に共分散を加算する必要があるようです。
ただその理由については，執筆者はまだ理解できておりません...。ぜひご教授ください...。

#おわりに
執筆者はいわゆる参加者内デザインで実験を計画することが多いです。このデザインでは，各実験参加者が複数の条件を経験することになります。個人の中に各条件が存在する，階層構造を持ったデータといえます。例えば以下の例では，各参加者が実験条件（experiment）と統制条件（control）を両方経験しています。

参加者内デザインの実験で媒介モデルを主張したい場合，マルチレベル媒介モデルを適用することになります。しかし，それを実現する方法がよく分からなかったので，Stanで書いてみようと思ったのが本記事の出発点でした。~~その過程でbmlmパッケージの存在を知り「なんだこれでいいじゃん」と一瞬思ったものの，~~試行錯誤しながら，ご助言をいただきながらStanで書いてみて，色々と勉強になりました。

もし本記事に間違いがありましたら，ご指摘を頂ければ幸いです。
最後に，本記事の執筆にあたり，ご助言をくださった方々に，心から御礼申し上げます。