前回、ベイズ推定の基本について野球選手の打率を例にとって学びました。今回はもう一歩踏み込み、学校教育を例にとって 階層ベイズモデル についてがんばって考えていきたいと思います 💪

使用するデータ

今回は Eight Schools というデータを使用します。これはアメリカの 8 つの高校で実施された SAT−V スコア向上のための特別な教育プログラムの平均値と標準誤差（※1）のデータです。例えば、学校 1 の平均と標準誤差はそれぞれ 28, 15 となります。

今回やりたいこと

上記の学校ごとの平均を、そのまま学校ごとの教育プログラムの効果としてもいいです。しかし今回は、「教育プログラムの自体の効果」というものが存在し、各学校のプログラム効果は教育プログラム自体の効果に似るはずだという仮説をもとに階層ベイズモデルを使用していきたいと思います。

このように、ドメイン知識などから個々の効果が全体の効果に近いと考える場合に階層ベイズモデルを選択すると、より妥当だと思える結果が得られる可能性があります。

全体像を把握する

階層ベイズでモデリングするにあたって今回の全体像を把握していきましょう。

今回使うデータは 8 つの各学校に在籍する生徒の平均的なプログラム受講前後の「スコア差」とその「標準誤差」です。これがどのような流れで得られたデータなのか考えてみましょう。

まず、すべての学校に共通する教育プログラム自体の効果（$${μ}$$）が存在すると考えます。今回は 8 校が選ばれていますが、それぞれの学校で行われる教育プログラムの効果（$${θ_i}$$）は、教師や設備の質などにより教育プログラム自体の効果（$${μ}$$）からばらつくと考えます。

$$
θ_i = \mu + \epsilon_i , \quad \epsilon_i \sim N(0, \tau^2)
$$

今回主に推定するのは、この $${\mu}$$, $${\tau}$$, $${\theta_i}$$ の値です。

階層ベイズの事後分布を考える

stanコードを書く前に、事前分布と尤度を整理しておきます。

事前分布

教育プログラム自体の効果（$${\mu}$$）とばらつき（$${\tau}$$）に関して、今回は事前知識がないため一様分布を事前分布とします。

$$
\mu, \tau \sim \text{\footnotesize 一様分布}
$$

尤度

尤度として以下を定義します。

$$
\begin{align*}
\tag{1} &\prod_{i = 1}^{8} p(\,θ_i\, | \,μ, τ\,) \\
\tag{2} &\prod_{i = 1}^{8} p(\,y_i \,|\, θ_i\,)
\end{align*}
$$

事後分布

上記の事前分布と尤度から、今回の事後分布は以下のように求められます。

$$
\begin{align*}
\text{\footnotesize 事後分布} \,\,\,&\propto\,\,\, p(μ) \,\, p(τ) \,\,\times \prod_{i = 1}^{8} p(\,θ_i\, | \,μ, τ\,) \,\,\times \prod_{i = 1}^{8} p(\,y_i \,|\, θ_i\,)
\end{align*}
$$

rstanで推定してみる

モデリングと推定

では今まで話してきたことを rstan に書き起こしてみましょう。

# stanコード

stancode <- "
data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  eta ~ normal(0, 1);        // prior log-density
  y ~ normal(theta, sigma);  // log-likelihood
}
"

ここで transformed parameters ブロックに注目したいと思います。

今回 $${\theta_i}$$を求める式として $${\theta_i = \mu + \tau \times \eta_i}$$が出てきます。これは reparameterization といい、$${θ_i}$$のばらつきを$${τ}$$, $${η_i}$$で表現しています。こうすることにより、乱数生成による推定量の作成が安定します（※2）。

$$
\begin{align*} θ_i &= μ + ε_i\,\,\,\,\,&(\,ε_i \sim N(0, τ)\,) \\
&=μ + τ \times η_i \,\,\,\,\,&(\,η_i \sim N(0, 1)\,)\end{align*}
$$

$${τ}$$は学校間の変動の大きさ（標準偏差）を表し、$${η_i}$$は$${N(0,1)}$$に従う$${μ}$$からの偏差を表します。

# データの準備と推定

schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

fit <- stan(model_code = stancode, data = schools_dat)

学校ごとのデータを準備し、先程の stan コードと合わせて stan 関数に引数として渡し、$${μ}$$, $${τ}$$, $${\theta_i}$$ を推定します。

結果の考察

それでは fit の結果を見てみましょう。

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%   25%   50%   75% 97.5% n_eff Rhat
mu        7.83    0.13 5.09  -2.46  4.77  7.83 11.07 17.63  1469    1
tau       6.74    0.17 5.59   0.28  2.60  5.40  9.37 21.51  1037    1
eta[1]    0.40    0.02 0.94  -1.50 -0.22  0.43  1.03  2.22  3480    1
eta[2]    0.00    0.01 0.88  -1.77 -0.57 -0.01  0.58  1.74  4151    1
eta[3]   -0.21    0.02 0.91  -1.92 -0.83 -0.23  0.39  1.61  3010    1
eta[4]   -0.05    0.01 0.88  -1.78 -0.62 -0.05  0.53  1.67  4034    1
eta[5]   -0.35    0.02 0.87  -1.97 -0.95 -0.37  0.19  1.44  3004    1
eta[6]   -0.20    0.01 0.89  -1.90 -0.80 -0.22  0.39  1.58  3907    1
eta[7]    0.35    0.01 0.88  -1.42 -0.21  0.37  0.93  2.05  3640    1
eta[8]    0.08    0.02 0.94  -1.74 -0.55  0.08  0.69  1.97  3767    1
theta[1] 11.64    0.17 8.67  -1.73  5.99 10.38 15.79 33.22  2619    1
theta[2]  7.85    0.08 6.10  -4.29  3.89  7.84 11.78 20.21  5156    1
theta[3]  6.10    0.14 7.75 -11.29  1.80  6.55 10.83 20.54  3010    1
theta[4]  7.50    0.11 6.59  -6.08  3.49  7.57 11.62 20.30  3742    1
theta[5]  5.03    0.11 6.42  -9.40  1.26  5.52  9.39 16.23  3509    1
theta[6]  6.12    0.11 6.80  -8.36  2.22  6.46 10.52 18.64  3688    1
theta[7] 10.89    0.12 6.84  -0.84  6.32 10.15 14.73 26.66  3398    1
theta[8]  8.53    0.15 7.99  -6.78  3.82  8.20 12.95 25.77  2811    1
lp__     -4.80    0.08 2.60 -10.43 -6.39 -4.50 -2.93 -0.54  1164    1

Samples were drawn using NUTS(diag_e) at Wed Dec 20 06:32:07 2023.
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).

実際得られている値と比較して、$${\theta_i}$$の推定結果が$${\mu}$$に縮小された様子を図でも確認しておきます。

draws <- as.matrix(fit)
draws <- draws[,c("theta[1]","theta[2]","theta[3]","theta[4]","theta[5]","theta[6]","theta[7]","theta[8]")]

mcmc_recover_intervals(draws, schools$y) + hline_at(mu, colour = "darkred")

各校の効果（$${\theta_i}$$）の推定結果がかなり赤線（$${\mu}$$）に縮小しているのがわかります。全体の教育効果に各校の効果が似るはずという仮説をモデルに反映することができました 👏

さいごに

こちらは「統計・機械学習の数理 Advent Calendar 2023」アドベントカレンダー 22日目の記事です。他の記事もぜひ覗いてみてくださいね。

今年はベイズ統計に入門できてとっても楽しかったです。いつかちゃんと勉強してみたいと思ってたので嬉しく思います。引き続きもう少しいろいろ例を考えたいと思うし、また別の統計モデルも勉強してみたい！

来年も楽しい一年になりそうです。みなさんの一年もいい年になりますように🎄

注釈

※1 Eight Schools の標準誤差

今回提供されているデータは標準偏差ではなく標準誤差です。

母平均 $${\mu}$$, 母分散$${\sigma^2}$$ のとき、中心極限定理によって標本平均$${\displaystyle \bar{X} = \frac{X_1 + …+X_n}{n}}$$は近似的に$${\displaystyle N\left (\mu, \frac{\sigma^2}{n} \right)}$$に従います。

標本平均 $${\bar{X}}$$の標準偏差$${\displaystyle \frac{\sigma}{\sqrt{n}}}$$を標本平均の標準誤差（standard error）といい、標本平均の推定値が母平均からどの程度ばらつくかを測ることができます。

※2 reparameterization による乱数生成の安定性

試しに以下のように parameterization しない stan コードを実行した場合、うまく収束していないと警告が表示されます。

stancode <- "
data {
  int<lower=0> J;         // number of prefectures 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] theta;
}
model {
  //target += normal_lpdf(eta | 0, 1);       // prior log-density
  //target += normal_lpdf(y | theta, sigma); // log-likelihood

  theta ~ normal(mu, tau);
  y ~ normal(theta, sigma);
}
"

fit <- stan(model_code = stancode, data = schools_dat)

※3 θ1のEAP推定量

最初、fit の結果を見たとき、単純に mu, tau, eta[1] の mean の値で計算（$${E(\mu) + E(\tau) \times E(\eta_1)}$$）した結果と theta[1] の推定結果が異なるのはなんでだろうと疑問を抱きました。

$$
\begin{align*}
E(\mu) + E(\tau) \times E(\eta_1) = 7.83 + 6.74 \times 0.40 = 10.526
\end{align*}
$$

fit の結果を見ると $${\theta_1}$$の推定結果は 11.64 なので、10.526 とは乖離があります。これは今回$${\theta_i}$$の推定値としてEAP推定量を用いているから生じる事象です。期待値の性質から、確率変数 $${X, Y}$$の和の期待値は期待値の和と同じです。しかし積の期待値が期待値の積と一致するとは限りません。

$$
\begin{align*}
&E(X+Y) = E(X) + E(Y) \\[5pt]
&E(XY) \neq E(X) E(Y)
\end{align*}
$$

つまり、$${\theta_1}$$のEAP推定量は以下のように書き下すことができます。

$$
E(\theta_1) = E(\mu + \tau \times \eta_1) = E(\mu) + E(\tau \times \eta_1)
$$

$${E(\tau \times \eta_1)}$$部分は各確率変数が独立でなければ $${E(\tau) \times E(\eta_1)}$$ と分解することはできません。今回、$${\tau}$$と$${\eta_1}$$が相関しているため、MCMCサンプリング時に $${\theta_1}$$を算出する際 $${\tau}$$が高いと$${\eta_1}$$も高くなる傾向があり、単純に $${E(\mu) + E(\tau) \times E(\eta_1)}$$ と計算した場合と異なる結果になりました。

実際に手計算で $${E(\mu) + E(\tau \times \eta_1)}$$ の結果を確かめてみましょう。$${\tau, \eta_1}$$の共分散や各々の変数の期待値はMCMCサンプリングの結果から算出した下記表の値を用います。

$$
\begin{align*}
E(\tau \times \eta_1) &= Cov(\tau, \eta_1) + E(\tau) \times E(\eta_1) \\[3pt]
&= 1.1414 + 6.7353 \times 0.3967= 3.8133\\[10pt]
E(\theta_1) &= E(\mu) + E(\tau \times \eta_1) \\[3pt]
&= 7.831 + 3.8133 = 11.6443 \\
\end{align*}
$$

fit の theta[1] の結果である「11.64」とほぼ一致する値が得られました。

Special Thanks

今回も記事を書くにあたってたくさんお世話になりました！
utaka233