はじめに

このノートは「StanとRでベイズ統計モデリング」の内容を写経することで、ベイズ統計への理解を深めていくために作成している。

基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、「StanとRでベイズ統計モデリング」を読み進めるための自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。

今回は第8章「階層モデル」のチャプターを写経していく。

8.1 階層モデルの導入

この章で使用するデータの説明しておく。大手企業４社\(KID\)の年齢\(X-23\)、年収\(Y\)が記録されている40人のデータを使用する。\(X\)から23が引かれているのは解釈がしやすくするため。

options(max.print = 999999)
library(dplyr)
library(ggplot2)
library(rstan)

d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap08/input/data-salary-2.txt')
d

##     X   Y KID
## 1   7 457   1
## 2  10 482   1
## 3  16 518   1
## 4  25 535   1
## 5   5 427   1
## 6  25 603   1
## 7  26 610   1
## 8  18 484   1
## 9  17 508   1
## 10  1 380   1
## 11  5 453   1
## 12  4 391   1
## 13 19 559   1
## 14 10 453   1
## 15 21 517   1
## 16 12 553   2
## 17 17 653   2
## 18 22 763   2
## 19  9 538   2
## 20 18 708   2
## 21 21 740   2
## 22  6 437   2
## 23 15 646   2
## 24  4 422   2
## 25  7 444   2
## 26 10 504   2
## 27  2 376   2
## 28 15 522   3
## 29 27 623   3
## 30 14 515   3
## 31 18 542   3
## 32 20 529   3
## 33 18 540   3
## 34 11 411   3
## 35 26 666   3
## 36 22 641   3
## 37 25 592   3
## 38 28 722   4
## 39 24 726   4
## 40 22 728   4

ここでも分析の問題は、年功序列で賃金は上昇する傾向にあるが、会社によって「新卒時点の基本年収」や「年齢の伴う昇給額は異なる」と考えられるため、その会社によるグループ差を検討したい。

8.1.2 グループ差を考えない場合

グループ差を考えない場合は、これまで通り通常の回帰分析を行えばよく、下記のモデルを元にパラメタを推定する。

モデル8-1

\[ \begin{eqnarray} Y[n] &=& y_{base}[n] + \epsilon[n] \\ y_{base} &=& a + bX[n] \\ \epsilon[n] &\sim& Normal(0, \sigma_{Y}) \end{eqnarray} \]

4社にわたる40人で共通の直線式があり、そこに他人ごとに独立にノイズが加わって\(Y\)が生成される、と仮定している。モデル式は下記の通りにも書ける。

\[ \begin{eqnarray} Y[n] &\sim& Normal(a + bX[n], \sigma_{Y}) \end{eqnarray} \] ただ、このデータを可視化してみるとわかるが、実態として、会社ごとに「新卒時点の基本年収」や「年齢の伴う昇給額は異なる」可能性が高く、このモデルはデータをあまり上手くは説明できていない。

d$KID <- as.factor(d$KID)
res_lm <- lm(Y ~ X, data = d)
coef <- coef(res_lm)
ggplot(d, aes(X, Y, shape = KID)) +
  theme_bw(base_size = 15) +
  geom_abline(intercept = coef[1], slope = coef[2], linewidth = 2, alpha = 0.3) +
  facet_wrap(~ KID) +
  geom_line(stat = 'smooth', method = 'lm', se = FALSE, color = 'black', linetype = '31') +
  geom_point(size = 3) +
  scale_shape_manual(values = c(16, 2, 4, 9)) +
  labs(x = 'X', y = 'Y')

8.2.2 グループごとに切片と傾きを持つ場合

会社ごとに「新卒時点の基本年収」や「年齢の伴う昇給額は異なる」可能性が高いのでモデルを改良する。つまり、会社ごとに\(a,b\)を推定すればよさそうに思える。会社の数を\(K(1\sim4)\)とすると、\(a[1],a[2],a[3],a[4],b[1],b[2],b[3],b[4]\)を推定する。

モデルは下記の通り。このモデルでは、ノイズの大きさはすべて共通で、\(KID=1\)のデータで\(a[1],b[1]\)を推定し、\(KID=2\)のデータで\(a[2],b[2]\)を推定することを繰り返す。

モデル8-2

\[ \begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \end{eqnarray} \]

データを準備しておく。これ以降\(KID\)などグループを表現する変数の動かし方が重要なので、そのあたりを注意しておくと個人的にはモデルを理解しやすくなると思われる。

N <- nrow(d)
K <- 4
data <- list(N = N, K = K, X = d$X, Y = d$Y, KID = as.numeric(d$KID))
data

## $N
## [1] 40
## 
## $K
## [1] 4
## 
## $X
##  [1]  7 10 16 25  5 25 26 18 17  1  5  4 19 10 21 12 17 22  9 18 21  6 15  4  7
## [26] 10  2 15 27 14 18 20 18 11 26 22 25 28 24 22
## 
## $Y
##  [1] 457 482 518 535 427 603 610 484 508 380 453 391 559 453 517 553 653 763 538
## [20] 708 740 437 646 422 444 504 376 522 623 515 542 529 540 411 666 641 592 722
## [39] 726 728
## 
## $KID
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4
## [39] 4 4

ここでは、stan_model()関数で最初にコンパイルしておいてから、

model82 <- stan_model('note_ahirubayes07-82.stan')

sampling()関数でサンプリングする。

fit <- sampling(object = model82, data = data, seed = 1989)

推定結果を見ると、\(a[1],a[2],a[3],a[4],b[1],b[2],b[3],b[4]\)が推定されていることがわかる。ただ、\(b[4]=-1.2\)と推定されているが、年齢が上がると年収が下がることになるため、本来はありえない。このように推定されてしまうのは、データのサンプルサイズが小さいため。\(KID=4\)はサンプルサイズが3しか無いため、データに対して素直に推定すると、このような推定をしてしまう。これでは分析したい目的とは異なるので、さらにモデルを改良する必要がある。また、会社ごとの推定値は得られたものの、\(a,b\)に関する仮定がないので、会社による差はまだわからないままであり、このデータに含まれていない会社の切片と傾きについては答えることができない。

print(fit, prob = c(0.025, 0.5, 0.975), digits_summary = 1)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean    sd   2.5%    50%  97.5% n_eff Rhat
## a[1]  387.1     0.3  14.3  359.6  386.9  416.1  2404    1
## a[2]  329.2     0.3  16.5  297.4  328.8  362.4  2979    1
## a[3]  312.5     0.6  34.5  244.5  311.9  380.5  2841    1
## a[4]  753.5     3.5 161.1  434.5  751.7 1075.8  2064    1
## b[1]    7.5     0.0   0.9    5.8    7.5    9.2  2495    1
## b[2]   19.8     0.0   1.2   17.4   19.9   22.2  2811    1
## b[3]   12.5     0.0   1.7    9.2   12.5   15.9  2842    1
## b[4]   -1.2     0.1   6.5  -14.1   -1.1   11.6  2165    1
## s_Y    27.4     0.1   3.6   21.4   27.0   35.9  2214    1
## lp__ -148.2     0.1   2.5 -154.2 -147.8 -144.5  1085    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan  3 13:09:20 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

次に進む前にモデルの挙動を確認しておく。下記のStanモデルの通り、\(KID=1\)のデータで\(a[1],b[1]\)を推定し、\(KID=2\)のデータで\(a[2],b[2]\)を推定することを繰り返したわけではあるが、Stanではどのようにデータをやり取りしているのか。

model {
  for (n in 1:N)
    Y[n] ~ normal(a[KID[n]] + b[KID[n]]*X[n], s_Y);
}

a[KID[n]]の部分を使って理解を深めておく。ここではあまり図解する必要はないかもしれないが、階層モデルを複雑にしていくと、個人的には直感的によく分からなくなるので、図解しておく。

    X   Y KID
1   7 457   1 | a[KID[n]]: KID[ 1] -> 1 -> a[1]
2  10 482   1 | a[KID[n]]: KID[ 2] -> 1 -> a[1]
[snip]
15 21 517   1 | a[KID[n]]: KID[15] -> 1 -> a[1]
16 12 553   2 | a[KID[n]]: KID[16] -> 2 -> a[2]
[snip]
27  2 376   2 | a[KID[n]]: KID[27] -> 2 -> a[2]
28 15 522   3 | a[KID[n]]: KID[28] -> 3 -> a[3]
[snip]
37 25 592   3 | a[KID[n]]: KID[37] -> 3 -> a[3]
38 28 722   4 | a[KID[n]]: KID[38] -> 4 -> a[4]
39 24 726   4 | a[KID[n]]: KID[39] -> 4 -> a[4]
40 22 728   4 | a[KID[n]]: KID[40] -> 4 -> a[4]

8.1.4 階層モデル

ここからはグループ差に仮定を入れた階層モデルを考える。メカニズムとしては、\(a[1],a[2],a[3],a[4]\)を「すべての会社で共通の全体平均」と「会社差を表す項」に分解する。会社差は平均0標準偏差\(\sigma_{a}\)の正規分布から生成されると考える。\(b[1],b[2],b[3],b[4]\)も同様で、年齢に伴う昇給額は「すべての会社で共通の全体平均」と「会社差を表す項」に分解する。

このような制約を入れることで、新卒時点の年収の会社差のばらつきは\(\sigma_{a}\)ぐらいと解釈できるようになる。\(\sigma_{a}\)には無情報事前分布を仮定する。このように階層的に事前分布が設定されているため階層モデルと呼ばれる。

モデル8-3

\[ \begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \\ a[k] &=& a_{全体平均} + a_{会社差}[k] \\ a_{会社差}[k] &\sim& Normal(0, \sigma_{a}) \\ b[k] &=& b_{全体平均} + b_{会社差}[k] \\ b_{会社差}[k] &\sim& Normal(0, \sigma_{b}) \\ \end{eqnarray} \]

このモデルでは、データから下記のパラメタを推定する。また、\(a_{全体平均},b_{全体平均},\sigma_{a},\sigma_{b},\sigma_{Y}\)の事前分布に正規分布を仮定する。

\(\sigma_{Y}\)
\(a_{全体平均}\)
\(a_{会社差}[1],a_{会社差}[2],a_{会社差}[3],a_{会社差}[4]\)
\(\sigma_{a}\)
\(b_{全体平均}\)
\(b_{会社差}[1],b_{会社差}[2],b_{会社差}[3],b_{会社差}[4]\)
\(\sigma_{b}\)

書籍にもあるとおり、シュミレーションデータを作ることで、データ生成過程への理解を深める。

set.seed(1)
N <- 40
K <- 4
N_k <- c(15, 12, 10, 3)
a0 <- 350　
b0 <- 12
s_a <- 60
s_b <- 4
s_Y <- 25
X <- sample(x = 0:35, size = N, replace = TRUE)
KID <- rep(1:4, times = N_k)
a <- a0 + rnorm(K, mean = 0, sd = s_a)
b <- b0 + rnorm(K, mean = 0, sd = s_b)
d <- data.frame(X = X, KID = KID, a = a[KID], b = b[KID])
d <- transform(d, Y_sim = rnorm(N, mean = a + b*X, sd = s_Y))
# 各会社の最初と最後のレコードだけ表示している。
d[c(1, 15, 16, 27, 28, 37, 38, 40), ]

##     X KID        a        b    Y_sim
## 1   3   1 333.2392 8.671827 388.1681
## 15 24   1 333.2392 8.671827 541.0869
## 16 33   2 455.4742 7.333718 673.9706
## 27 22   2 455.4742 7.333718 613.8368
## 28  5   3 383.6448 7.737638 438.9363
## 37 21   3 383.6448 7.737638 540.6473
## 38 13   4 322.8330 5.744872 386.8960
## 40  5   4 322.8330 5.744872 376.4820

この部分でデータを生成する際のパラメタを設定している。レコード数は\(N=40\)、会社の数\(K=4\)やばらつき\(N_{k}=(15, 12, 10, 3)\)、全体平均\(a_{0}=350, b_{0}=12\)、標準偏差\(\sigma_{a}=60, \sigma_{b}=4,\sigma_{Y}=25\)である。

set.seed(123)
N <- 40
K <- 4
N_k <- c(15, 12, 10, 3)
a0 <- 350　
b0 <- 12
s_a <- 60
s_b <- 4
s_Y <- 25
X <- sample(x = 0:35, size = N, replace = TRUE)
KID <- rep(1:4, times = N_k)

そして、\(a,b\)ともに、例えば\(a_{1}\)であれば、全体平均\(a_{0} + a[1], \ a[1] \sim Normal(0, \sigma_{a})\)の関係から生成される。

a <- a0 + rnorm(K, mean = 0, sd = s_a)
a
[1] 342.2354 403.2042 340.9162 369.7875

b <- b0 + rnorm(K, mean = 0, sd = s_b)
b
[1] -0.9092913  8.9128329 13.1461943  7.1179521

最終的に、線形モデルの各パラメタとして利用され、ノイズ\(\sigma_{Y}\)が加わって年収\(Y\)が生成される。

Y_sim = rnorm(N, mean = a + b*X, sd = s_Y)

d$KID <- as.factor(d$KID)
ggplot(d, aes(X, Y_sim, shape = KID)) +
  theme_bw(base_size = 15) +
  facet_wrap(~KID) +
  geom_line(stat = 'smooth', method = 'lm', se = FALSE, linewidth = 1, color = 'black', linetype = '31') +
  geom_point(size = 3) +
  scale_shape_manual(values = c(16, 2, 4, 9)) +
  labs(x = 'X', y = 'Y')

データは先程と同じものを利用する。

d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap08/input/data-salary-2.txt')
N <- nrow(d)
K <- 4
data <- list(N = N, K = K, X = d$X, Y = d$Y, KID = as.numeric(d$KID))
data

## $N
## [1] 40
## 
## $K
## [1] 4
## 
## $X
##  [1]  7 10 16 25  5 25 26 18 17  1  5  4 19 10 21 12 17 22  9 18 21  6 15  4  7
## [26] 10  2 15 27 14 18 20 18 11 26 22 25 28 24 22
## 
## $Y
##  [1] 457 482 518 535 427 603 610 484 508 380 453 391 559 453 517 553 653 763 538
## [20] 708 740 437 646 422 444 504 376 522 623 515 542 529 540 411 666 641 592 722
## [39] 726 728
## 
## $KID
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4
## [39] 4 4

階層モデルのStanファイルはこちら。

data {
  int N;
  int K;
  real X[N];
  real Y[N];
  int<lower=1, upper=K> KID[N];
}

parameters {
  real a0;
  real b0;
  real ak[K];
  real bk[K];
  real<lower=0> s_a;
  real<lower=0> s_b;
  real<lower=0> s_Y;
}

transformed parameters {
  real a[K];
  real b[K];
  for (k in 1:K) {
    a[k] = a0 + ak[k];
    b[k] = b0 + bk[k];
  }
}

model {
  s_a ~ normal(0, 1e5);
  s_b ~ normal(0, 1e5);
  s_Y ~ normal(0, 1e5);
  
  for (k in 1:K) {
    ak[k] ~ normal(0, s_a);
    bk[k] ~ normal(0, s_b);
  }

  for (n in 1:N)
    Y[n] ~ normal(a[KID[n]] + b[KID[n]]*X[n], s_Y);
}

ここでは、stan_model()関数で最初にコンパイルしておいてから、

model83 <- stan_model('note_ahirubayes07-83.stan')

sampling()関数でサンプリングする。

fit <- sampling(object = model83, data = data, seed = 1989)

推定結果を見ると、「新卒給与」の会社差のばらつき\(\sigma_{a}\)はおよそ169万円、「年齢昇給額」の会社差のばらつき\(\sigma_{b}\)はおよそ17万円、ノイズの大きさは\(\sigma_{Y}\)はおよそ28万円とわかる。

\(a_{全体平均}\): 383.4[148.7~655.5]
\(b_{全体平均}\): 14.3 [-8.9~45.5]
\(\sigma_{a}\): 169.0[11.8~728.8]
\(\sigma_{b}\): 17.3[3.2~101.0]
\(\sigma_{Y}\): 28.4[22.2~36.7]

print(fit, prob = c(0.025, 0.5, 0.975), digits_summary = 1)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean    sd   2.5%    50%  97.5% n_eff Rhat
## a0     383.4     6.9 136.2  148.7  368.7  655.5   393  1.0
## b0      14.3     2.7  20.7   -8.9   12.5   45.5    60  1.0
## ak[1]    0.1     6.9 136.1 -274.3   11.1  243.1   394  1.0
## ak[2]  -48.5     6.9 136.9 -325.7  -31.7  184.1   399  1.0
## ak[3]  -58.7     7.0 138.9 -340.6  -39.0  174.9   388  1.0
## ak[4]  105.0     7.1 167.7  -97.0   65.2  519.7   554  1.0
## bk[1]   -6.5     2.7  20.7  -37.4   -4.7   16.6    60  1.0
## bk[2]    5.1     2.7  20.7  -25.4    6.7   28.6    60  1.0
## bk[3]   -2.3     2.7  20.8  -33.2   -0.7   21.0    60  1.0
## bk[4]   -4.8     2.7  21.1  -39.0   -1.8   14.2    62  1.1
## s_a    169.0    10.4 232.1   11.8  102.0  728.8   495  1.0
## s_b     17.3     3.5  40.0    3.2    8.5  101.0   129  1.0
## s_Y     28.4     0.1   3.7   22.2   28.0   36.7  1646  1.0
## a[1]   383.5     0.3  14.9  353.6  383.8  412.4  2331  1.0
## a[2]   334.9     0.3  17.3  301.9  335.0  368.9  3128  1.0
## a[3]   324.7     0.6  33.2  257.8  325.9  385.5  3103  1.0
## a[4]   488.4     5.9 139.9  298.7  454.9  823.1   555  1.0
## b[1]     7.7     0.0   0.9    5.9    7.7    9.6  3021  1.0
## b[2]    19.4     0.0   1.3   16.8   19.4   21.9  3125  1.0
## b[3]    11.9     0.0   1.6    8.9   11.9   15.2  3392  1.0
## b[4]     9.5     0.2   5.6   -3.9   10.8   17.2   552  1.0
## lp__  -173.4     0.2   3.7 -181.9 -172.9 -167.4   229  1.0
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan  3 13:09:45 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

1つ前のモデルでは、b[4] = -1.2と推定されていたが、階層モデルではb[4] = 9.5となっている。これは、\(b[1],b[2],b[3],b[4]\)の年齢に伴う昇給額は「すべての会社で共通の全体平均\(b_{0}\)」と「会社差を表す項\(b_{k}\)」に分解され、会社差\(b_{k}\)には平均0標準偏差\(\sigma_{b}\)の正規分布から生成されると考える、という仮定を置いている効果がここに出ている。

この階層モデルでもモデルの挙動を確認しておく。

a[KID[n]]の部分を使って理解を深めておく。この部分でnormal(0, s_a)に従うak[1]-ak[4]が生成される。

for (k in 1:K) {
  ak[k] ~ normal(0, s_a);
  bk[k] ~ normal(0, s_b);
}

// image 
ak[1] → AK1
ak[2] → AK2
ak[3] → AK3
ak[4] → AK4

そして、この部分でさきほどのak[1]-ak[4]が利用され、全体平均a0が加えられることで、a[1]-a[4]が生成される。

for (k in 1:K) {
    a[k] = a0 + ak[k];
    b[k] = b0 + bk[k];
  }

// image 
a0 + AK1 -> a[1] -> A1 
a0 + AK2 -> a[2] -> A2
a0 + AK3 -> a[3] -> A3
a0 + AK4 -> a[4] -> A4

生成されたa[1]-a[4](同様にb[1]-b[4])を利用して、正規分布の部分に情報が渡され\(Y\)が生成される。

for (n in 1:N){
  Y[n] ~ normal(a[KID[n]] + b[KID[n]]*X[n], s_Y);
}

// aのイメージ、bは省略

    X   Y KID
1   7 457   1 | a[KID[n]]: KID[ 1] -> 1 -> a[1] -> A1 | Normal(A1 + B1*X[ 1], s_Y)
2  10 482   1 | a[KID[n]]: KID[ 2] -> 1 -> a[1] -> A1 | Normal(A1 + B1*X[ 2], s_Y)
[snip]
15 21 517   1 | a[KID[n]]: KID[15] -> 1 -> a[1] -> A1 | Normal(A1 + B1*X[15], s_Y)
16 12 553   2 | a[KID[n]]: KID[16] -> 2 -> a[2] -> A2 | Normal(A2 + B2*X[16], s_Y)
[snip]
27  2 376   2 | a[KID[n]]: KID[27] -> 2 -> a[2] -> A2 | Normal(A2 + B2*X[27], s_Y)
28 15 522   3 | a[KID[n]]: KID[28] -> 3 -> a[3] -> A3 | Normal(A3 + B3*X[28], s_Y)
[snip]
37 25 592   3 | a[KID[n]]: KID[37] -> 3 -> a[3] -> A3 | Normal(A3 + B3*X[37], s_Y)
38 28 722   4 | a[KID[n]]: KID[38] -> 4 -> a[4] -> A4 | Normal(A4 + B4*X[38], s_Y)
39 24 726   4 | a[KID[n]]: KID[39] -> 4 -> a[4] -> A4 | Normal(A4 + B4*X[39], s_Y)
40 22 728   4 | a[KID[n]]: KID[40] -> 4 -> a[4] -> A4 | Normal(A4 + B4*X[40], s_Y)

8.1.6 階層モデルの等価な表現

モデル8-4はモデル8-3と等価である。すべての会社の平均が\(a_{全体平均}\)で、各社の\(a[k]\)が加わって生成されるという形から、平均\(a_{全体平均}\)で、標準偏差\(\sigma_{a}\)の正規分布から生成される形となっている。

モデル8-4

\[ \begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \\ a[k] &\sim& Normal(a_{全体平均}, \sigma_{a}) \\ b[k] &\sim& Normal(b_{全体平均}, \sigma_{b}) \\ \end{eqnarray} \]

Stanのモデルは下記の通り。

data {
  int N;
  int K;
  real X[N];
  real Y[N];
  int<lower=1, upper=K> KID[N];
}

parameters {
  real a0;
  real b0;
  real a[K];
  real b[K];
  real<lower=0> s_a;
  real<lower=0> s_b;
  real<lower=0> s_Y;
}

model {
  s_a ~ normal(0, 1e5);
  s_b ~ normal(0, 1e5);
  s_Y ~ normal(0, 1e5);
  
  for (k in 1:K) {
    a[k] ~ normal(a0, s_a);
    b[k] ~ normal(b0, s_b);
  }

  for (n in 1:N)
    Y[n] ~ normal(a[KID[n]] + b[KID[n]]*X[n], s_Y);
}

StanとRでベイズ統計モデリング07