RStan:階層ベイズモデル

# シミュレーションのパラメータ設定
J <- 10 # 校舎の数
N_j <- 20 # 各校舎の生徒数
N <- J * N_j # 全生徒数

# 真のパラメータ
mu_true <- 60 # 全校舎の平均学力の平均（グランド平均）
sigma_g_true <- 10 # 校舎間の学力レベルのばらつき
sigma_s_true <- 15 # 生徒間の成績のばらつき（校舎内SD）

seed <- 20250701
set.seed(seed)

# 1. 各校舎の真の平均学力（theta）を生成
theta_true <- rnorm(n = J, mean = mu_true, sd = sigma_g_true)

# 2. 各生徒のデータを生成
school_id <- rep(1:J, each = N_j)
score <- numeric(N)
for (i in 1:N) {
  score[i] <- rnorm(n = 1, mean = theta_true[school_id[i]], sd = sigma_s_true)
}

# データフレームにまとめる
sim_data <- data.frame(
  student_id = 1:N,
  school_id = factor(school_id),
  score = score
)

cat("--- 生成したデータの一部を確認 ---\n")
head(sim_data)

# 校舎ごとの単純な平均点と真の値をプロット
library(ggplot2)
library(dplyr)
sim_data %>%
  group_by(school_id) %>%
  summarise(raw_mean_score = mean(score)) %>%
  mutate(true_theta = theta_true) %>%
  ggplot(aes(x = school_id)) +
  geom_point(aes(y = raw_mean_score, colour = "観測された平均点"), size = 4, alpha = 0.8) +
  geom_point(aes(y = true_theta, colour = "真の平均学力"), size = 4, shape = 17) +
  labs(
    title = "各校舎の観測平均点と真の平均学力",
    x = "校舎ID", y = "点数", colour = "種類"
  ) +
  scale_color_manual(values = c("観測された平均点" = "dodgerblue", "真の平均学力" = "tomato")) +
  ylim(0, 100)

# MCMCサンプルからthetaの推定値を取得
theta_est <- rstan::extract(fit, pars = "theta")$theta

# thetaの事後平均と95%信用区間を計算
theta_summary <- data.frame(
  school_id = 1:J,
  mean = apply(theta_est, 2, mean),
  q2.5 = apply(theta_est, 2, quantile, probs = 0.025),
  q97.5 = apply(theta_est, 2, quantile, probs = 0.975)
)

# 観測された単純平均と真の値を結合
plot_data <- sim_data %>%
  group_by(school_id) %>%
  summarise(raw_mean = mean(score)) %>%
  mutate(school_id = as.integer(school_id)) %>%
  left_join(theta_summary, by = "school_id") %>%
  mutate(true_theta = theta_true)

# 収縮プロット
ggplot(plot_data, aes(x = raw_mean, y = mean)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, color = "skyblue") +
  geom_point(color = "dodgerblue", size = 3) +
  geom_hline(yintercept = mean(rstan::extract(fit, pars = "mu")$mu), linetype = "dotted", color = "red") +
  annotate("text",
    x = 55, y = mean(rstan::extract(fit, pars = "mu")$mu) + 1,
    label = "全体の平均 (mu)", color = "red"
  ) +
  labs(
    title = "階層ベイズによる収縮効果",
    subtitle = "単純平均（横軸）がベイズ推定値（縦軸）に変換される様子",
    x = "観測された単純平均点",
    y = "ベイズ推定された校舎の平均学力（事後平均）"
  )

ggplot(plot_data) +
  geom_point(aes(x = school_id, y = raw_mean, colour = "観測平均"), size = 3) +
  geom_point(aes(x = school_id, y = true_theta, colour = "真の値"), size = 4) +
  geom_pointrange(aes(x = school_id, y = mean, ymin = q2.5, ymax = q97.5, colour = "ベイズ推定値"), size = 0.8) +
  scale_x_continuous(breaks = 1:J) +
  scale_color_manual(
    name = "種類",
    values = c("観測平均" = "gray", "真の値" = "tomato", "ベイズ推定値" = "dodgerblue"),
    guide = guide_legend(override.aes = list(
      linetype = c("blank", "blank", "solid")
    ))
  ) +
  labs(
    title = "各校舎の平均学力の真の値と推定値の比較",
    x = "校舎ID", y = "点数"
  )