分析の目的

8つの学校で行ったコーチングの効果に関する研究データから、学校ごとの“真の”コーチング成果を推定する。

データの概要

学校ごとに教師と生徒は異なるがコーチングの内容は同じ。Yはそれぞれの学校のコーチング成果の推定量。SIGMAはコーチング成果の標準誤差。

library("rstan")
options(mc.cores = parallel::detectCores()) # マルチコアCPUでメモリも十分にある場合はこのオプションを追加
rstan_options(auto_write = TRUE) # Stanのモデルが変更されていない場合に再コンパイルされないようにする。
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))
summary(schools_dat)
##       Length Class  Mode   
## J     1      -none- numeric
## Y     8      -none- numeric
## SIGMA 8      -none- numeric

モデル式

このデータから学校1の真のコーチング成果をどの程度に見積もるのがよいだろうか？
①学校1の成果の推定量は28なのでそれをそのまま学校1の真の成果と見積もる。
②8つの学校の成果の平均値は8.75なのでそれを学校1の真の成果と見積もる。

①の場合は同じコーチングをしているにもかかわらず学校1以外のデータを無視していることになる。 ②の場合は教師と生徒が学校ごとに異なるにもかかわらず全ての学校で同じ効果であったとみなすことになる。

このような極端なふたつの仮定の中間をとって次のように考える。
- 真の推定量\(\theta_{j}\)は平均\(\mu\)、標準偏差\(\tau\)の正規分布から得られる。
- データ\(Y\)は平均\(\theta_{j}\)、標準偏差\(\Sigma\)の正規分布から得られる。

データの生成メカニズムに関する上記の仮定を式で表すと次のようになる。
\[ \begin{align} \eta_{j} & \sim Normal(0, 1) \\ \theta_{j} & = \mu+\tau\eta_{j} \\ Y_{j} & \sim Normal(\theta_{j}, \Sigma) \end{align} \]

パラメータ推定実行

fit <- stan(file = 'EightSchools.stan', data = schools_dat, iter = 1000, chains = 4)

推定結果

print(fit)
## Inference for Stan model: EightSchools.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu         8.05    0.18 5.29  -2.37   4.79   7.99  11.39  18.61   837 1.00
## tau        6.79    0.21 5.47   0.26   2.74   5.53   9.49  20.62   670 1.01
## eta[1]     0.43    0.02 0.92  -1.48  -0.15   0.49   1.02   2.14  2000 1.00
## eta[2]     0.03    0.02 0.86  -1.71  -0.56   0.01   0.59   1.81  2000 1.00
## eta[3]    -0.18    0.02 0.95  -2.03  -0.84  -0.20   0.47   1.73  2000 1.00
## eta[4]    -0.05    0.02 0.89  -1.81  -0.64  -0.05   0.53   1.74  2000 1.00
## eta[5]    -0.38    0.02 0.87  -1.98  -0.97  -0.41   0.15   1.44  2000 1.00
## eta[6]    -0.23    0.02 0.88  -1.90  -0.82  -0.26   0.34   1.54  1770 1.00
## eta[7]     0.34    0.02 0.91  -1.56  -0.25   0.37   0.97   2.06  1797 1.00
## eta[8]     0.08    0.02 0.96  -1.86  -0.56   0.08   0.74   1.97  2000 1.00
## theta[1]  11.90    0.23 8.47  -1.77   6.10  10.72  16.19  32.51  1387 1.00
## theta[2]   8.02    0.14 6.38  -4.57   4.04   7.74  11.86  21.32  2000 1.00
## theta[3]   6.24    0.18 7.93 -12.50   2.05   6.78  11.01  20.49  2000 1.00
## theta[4]   7.61    0.15 6.68  -6.11   3.43   7.70  11.75  21.06  2000 1.00
## theta[5]   5.01    0.14 6.47  -9.40   1.22   5.40   9.24  16.92  2000 1.00
## theta[6]   6.12    0.15 6.82  -7.96   2.07   6.40  10.61  18.78  2000 1.00
## theta[7]  10.79    0.15 6.88  -1.17   6.17  10.27  14.86  25.53  2000 1.00
## theta[8]   8.83    0.22 8.25  -7.41   4.15   8.56  13.17  26.30  1381 1.00
## lp__     -39.56    0.11 2.70 -45.62 -41.28 -39.31 -37.61 -35.07   592 1.01
## 
## Samples were drawn using NUTS(diag_e) at Tue Jul 31 23:53:56 2018.
## 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).
# パラメータの事後分布
plot(fit, pars=c("mu", "tau", "eta","theta"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

# パラメータのペアプロット
pairs(fit, pars = c("mu", "tau", "lp__"))