事後予測チェック
シミュレーションデータと観測データの分布を比較して、モデルの適合度を検証します。
「事後予測チェック」はCoddyKit上の無料R Academyレッスンです。 これはレッスン4/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
事後予測チェックとは
ベイズモデルをフィッティングした後は、「このモデルは、観測したデータに似たデータを生成するか」と問う必要があります。事後予測チェック（PPC）では、事後分布から複製データセットyrepをシミュレーションし、観測データyとグラフで比較することで、この問いに答えます。
事後サンプルの抽出
extract(fit, pars='mu')は名前付きリストを返します。$muには、そのパラメータのウォームアップ後の全サンプルを含む数値ベクトルが格納されます。4チェーンでウォームアップ後の反復が1000回なら、サンプルは4000個得られます。
# library(rstan)
# mu_samples <- extract(fit, pars = 'mu')$mu
# sigma_samples <- extract(fit, pars = 'sigma')$sigma
#
# cat('Samples drawn:', length(mu_samples), '
')
# cat('Posterior mean of mu:', mean(mu_samples), '
')
# cat('90% CI:', quantile(mu_samples, c(0.05, 0.95)), '
')事後サンプルからyrepを生成する
各事後抽出値(mu_s, sigma_s)について、元のデータと同じサイズの複製データセットをシミュレーションします。これらを行列yrepに保存します。各行が1つのシミュレーションデータセットを表します。
# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# n_obs <- length(y)
# S <- length(mu_samples) # 4000 posterior draws
#
# yrep <- matrix(NA, nrow = S, ncol = n_obs)
# for (s in seq_len(S)) {
# yrep[s, ] <- rnorm(n_obs, mean = mu_samples[s], sd = sigma_samples[s])
# }
# dim(yrep) # [4000, 8]ppc_dens_overlay() — 密度の比較
bayesplot::ppc_dens_overlay(y, yrep[1:50,])は、観測データのカーネル密度（濃い線）に、ランダムに選んだ50個のシミュレーションデータセットの密度（薄い線）を重ねて表示します。モデルの適合度が良ければ、濃い線は薄い線の集まりの中に収まります。
# library(bayesplot)
#
# ppc_dens_overlay(y, yrep[1:50, ])
#
# Interpretation:
# - Dark line (y_obs) surrounded by light lines (yrep): good fit
# - Dark line systematically outside the cloud: model misfit
# - Light lines much wider than dark: overdispersed model
# - Light lines much narrower than dark: underdispersed modelppc_stat() — 検定統計量の確認
ppc_stat(y, yrep, stat = 'mean')は、各シミュレーションデータセットで計算した検定統計量（例：平均）のヒストグラムを表示し、観測された統計量の位置に垂直線を引きます。観測値がヒストグラムの大部分に入っていれば、モデルはデータのその側面を適切に捉えています。
# library(bayesplot)
#
# ppc_stat(y, yrep, stat = 'mean') # does model capture the mean?
# ppc_stat(y, yrep, stat = 'sd') # does model capture spread?
# ppc_stat(y, yrep, stat = 'max') # does model capture extremes?
#
# If observed stat is in the tail of the histogram,
# the model fails to reproduce that statistic.ベイズp値
ベイズp値（事後予測p値）は、検定統計量が観測値よりも極端になるシミュレーションデータセットの割合です。0.5に近い値は適切なキャリブレーションを示し、0または1に近い値は、その統計量についてモデルがデータに適合していないことを示します。
# Bayesian p-value for the mean:
# obs_mean <- mean(y)
# rep_means <- apply(yrep, 1, mean)
# pval <- mean(rep_means >= obs_mean)
# cat('Bayesian p-value (mean):', round(pval, 3), '
')
# # 0.5 is perfect; < 0.05 or > 0.95 suggests misfitその他のbayesplot PPC関数
bayesplotには、密度の重ね合わせ以外にも多くのPPC可視化機能があります。
ppc_hist(y, yrep[1:8,])— ヒストグラムのグリッドppc_scatter_avg(y, yrep)— 観測値とyrepの平均値の散布図ppc_intervals(y, yrep)— 各観測値の周囲の不確実性区間ppc_rootogram(y, yrep)— カウントデータ用
# library(bayesplot)
#
# # Grid of 8 simulated histograms vs the observed
# ppc_hist(y, yrep[1:8, ])
#
# # Scatter: y_obs (x) vs mean of yrep (y) — should hug diagonal
# ppc_scatter_avg(y, yrep)
#
# # 50% and 90% posterior predictive intervals around each y_i
# ppc_intervals(y, yrep)PPCプロットの解釈 — モデルの不適合
よくある不適合のパターンと原因：
- yrepが広すぎる — 事前分布が広すぎるか、モデルが過分散です
- yrepがずれている — 尤度の分布族が適切ではありません（例：歪んだデータに正規分布を使用）
- yrepが多峰性を捉えられない — 混合モデルが必要です
- yrepが極端な値に対応できない — 裾の重い分布が必要です
# Example: if data has a long right tail but yrep does not,
# consider switching:
# y ~ normal(mu, sigma) => y ~ student_t(nu, mu, sigma)
#
# Or for count data:
# y ~ poisson(lambda) => y ~ neg_binomial_2(mu, phi) (overdispersion)
cat('PPCs guide model improvement by revealing specific failure modes
')StanのmodelブロックでPPCを行う
generated quantitiesブロックを使えば、Stan内で直接yrepを生成できます。これにより、Rでパラメータを再抽出する必要がなくなり、計算上は同等の処理を実行できます。
# Stan model with generated quantities:
# '
# generated quantities {
# array[N] real y_rep;
# for (n in 1:N) {
# y_rep[n] = normal_rng(mu, sigma);
# }
# }
# '
# Then extract in R:
# yrep <- extract(fit, pars = 'y_rep')$y_rep # [S, N] matrixLeave-One-Out交差検証
PPCに加えて、loo::loo(fit)を使うとleave-one-out交差検証を実行し、競合するモデルを比較できます。ELPD（期待対数予測密度）が高いモデルがより望ましいモデルです。loo::loo_compare(loo1, loo2)を使うと、モデルを順位付けできます。
# library(loo)
# loo1 <- loo(fit1) # normal model
# loo2 <- loo(fit2) # student-t model
#
# comparison <- loo_compare(loo1, loo2)
# print(comparison)
#
# Model with elpd_diff > 0 is preferred
# se_diff > |elpd_diff| means difference is not reliablePPCのベストプラクティス
厳密な事後予測チェックを行うには、次の方法に従ってください。
- 全体的な妥当性確認として、必ず
ppc_dens_overlay()から始めます - 分析目的に関連する、分野固有の統計量（
ppc_stat()）で追加確認します - 視覚的な確認には少なくとも50個のyrep抽出値を、ベイズp値には4000個すべてを使用します
- PPCの失敗はモデル改善の手がかりです — 失敗そのものではなく、診断結果です
# Workflow:
# 1. Fit model -> extract() -> generate yrep matrix
# 2. ppc_dens_overlay(y, yrep[1:50,]) -- visual global check
# 3. ppc_stat(y, yrep, stat='mean') -- check mean
# 4. ppc_stat(y, yrep, stat='sd') -- check spread
# 5. ppc_stat(y, yrep, stat='max') -- check tails
# 6. If misfit found -> revise model -> refit -> re-check確認問題：ベイズp値の解釈
最大値を検定統計量としたベイズp値が0.03だった場合、モデルについて何が分かりますか。
事後予測チェックのまとめ
RStanとbayesplotによるPPCのワークフロー：
- サンプルを抽出します：
extract(fit, pars='mu')$mu - yrepを生成します：事後抽出値ごとに
rnorm(n, mu_s, sigma_s)を呼び出すループを実行します - 全体的な確認：
ppc_dens_overlay(y, yrep[1:50,]) - 統計量の確認：
ppc_stat(y, yrep, stat='mean') - ベイズp値：
mean(apply(yrep,1,stat) >= stat(y))— 0.5に近ければ良好です - モデル比較には
loo::loo_compare()を使用します
よくある質問
「事後予測チェック」レッスンは無料ですか？
はい。「事後予測チェック」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「事後予測チェック」で何を学びますか？
シミュレーションデータと観測データの分布を比較して、モデルの適合度を検証します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか？
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン4/4です。
「事後予測チェック」レッスンにはどのくらい時間がかかりますか？
ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか？
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。