このページの目的

このページでは、BayesRTMB を使い始めるための最小限の流れを確認します。

ここで扱うのは、次の内容です。

1. rtmb_code() で小さなモデルを書く
2. rtmb_model() でモデルオブジェクトを作る
3. optimize() と sample() で推定する
4. MCMC の事後分布を可視化する
5. ラッパー関数で重回帰と交互作用を扱う
6. classic() による頻度主義的な t 検定を行う
7. JZS prior と MCMC によって Bayes factor を計算する

詳しいモデルコードの書き方や、混合モデル、GLMM、モデル比較の詳細は、別の vignette で扱います。

0. インストールと環境確認

BayesRTMB は CRAN からインストールできます。

install.packages("BayesRTMB")

開発版は GitHub からインストールできます。 pak を使う場合は、次のようにします。

install.packages("pak")
pak::pak("norimune/BayesRTMB")

remotes を使う場合は、次のようにします。

install.packages("remotes")
remotes::install_github("norimune/BayesRTMB")

pkgbuild::check_build_tools(debug = TRUE)

1. 最小モデルを書く

まず、もっとも小さな例として二項モデルを書きます。 ここでは、10 回の試行のうち成功が 6 回観測された状況を考えます。

library(BayesRTMB)

Trial <- 10
Y <- 6

dat <- list(Trial = Trial, Y = Y)

code <- rtmb_code(
  parameters = {
    theta <- Dim(lower = 0, upper = 1)
  },
  model = {
    Y ~ binomial(Trial, theta)
    theta ~ beta(1, 1)
  }
)

parameters ブロックでは推定するパラメータを宣言します。 ここでは成功確率 theta を、0 から 1 の範囲を持つパラメータとして定義しています。

model ブロックでは、観測データの分布と事前分布を書きます。 Y ~ binomial(Trial, theta) は、成功数 Y が二項分布に従うことを表しています。

2. モデルオブジェクトを作る

rtmb_model() にデータとモデルコードを渡すと、推定用のモデルオブジェクトが作られます。

mdl <- rtmb_model(dat, code)

## Pre-checking model code...
## Checking RTMB setup...

この段階では、まだ推定は行われていません。 mdl は、モデル定義とデータを保持した RTMB_Model オブジェクトです。

3. MAP 推定を行う

点推定をすばやく得たいときは、optimize() を使います。 BayesRTMB では、事前分布がある場合は MAP 推定、flat prior の場合は最尤推定に近い推定として扱えます。

fit_map <- mdl$optimize()
fit_map

## Starting RTMB optimization...
## 
## 
## Call:
## MAP Estimation via RTMB
## 
## Negative Log-Posterior: 1.38
## Approx. Log Marginal Likelihood (Laplace): -2.33
## 
## Point Estimates and 95% Wald CI:
## variable  Estimate  Std. Error  Lower 95%  Upper 95% 
## theta      0.60000     0.15492    0.29740    0.84166 

この例では、成功確率 theta の推定値はちょうど 0.60 です。

4. MCMC で事後分布を見る

事後分布全体を見たいときは、sample() を使います。 クイックスタートでは短い設定にしていますが、実際の分析ではより多くの反復数を使ってください。

set.seed(1)

fit_mcmc <- mdl$sample(
  sampling = 200,
  warmup = 200,
  chains = 2
)

fit_mcmc$summary()

## Starting sequential sampling (chains = 2)...
## chain 1 started... 
## chain 1: iter 100 warmup 
## chain 1: iter 200 warmup 
## chain 1: iter 300 sampling 
## chain 1: iter 400 sampling 
## chain 2 started... 
## chain 2: iter 100 warmup 
## chain 2: iter 200 warmup 
## chain 2: iter 300 sampling 
## chain 2: iter 400 sampling 
## variable   mean    sd    map   q2.5  q97.5  ess_bulk  ess_tail  rhat 
## lp        -3.30  0.69  -2.88  -5.25  -2.80       132       188  1.00 
## theta      0.58  0.13   0.60   0.32   0.83       165       177  1.01 

MCMC の結果では、平均、標準偏差、事後分位点、ESS、R-hat などを確認できます。 R-hat が 1 に近く、ESS が十分に大きいほど、MCMC の診断としては安心しやすくなります。

install.packages(c("future", "future.apply", "progressr"))

fit_mcmc <- mdl$sample(
  sampling = 1000,
  warmup = 1000,
  chains = 4,
  parallel = TRUE
)

5. MCMC の結果を可視化する

MCMC の結果は、数値だけでなく図でも確認します。 draws() で対象パラメータのサンプルを取り出し、密度、トレース、自己相関、区間推定を確認できます。

theta_draws <- fit_mcmc$draws("theta")

plot_dens(theta_draws)
plot_trace(theta_draws)
plot_acf(theta_draws)

それぞれの図は、次の目的で使います。

• plot_dens() は、事後分布の形を確認します。
• plot_trace() は、チェインがよく混ざっているかを確認します。
• plot_acf() は、自己相関が強すぎないかを確認します。

実際には、たとえば次のような図として確認できます。

Posterior density plot

Trace plot

Autocorrelation plot

6. ラッパー関数で重回帰を行う

標準的な分析では、rtmb_code() を自分で書かずにラッパー関数から始められます。

ここでは、debate データを使い、満足度 sat を talk、perf、およびその交互作用で説明する重回帰モデルを推定します。

data(debate)

mdl_lm <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  gmc = "all",
  prior = prior_normal()
)

fit_lm <- mdl_lm$optimize(
  se_method = "sampling",
  num_samples = 1000,
  seed = 1
)
fit_lm

## Starting RTMB optimization...
## 
## Using simulation-based error propagation (1000 samples)...
## 
## 
## Call:
## MAP Estimation via RTMB
## 
## Negative Log-Posterior: 402.24
## Approx. Log Marginal Likelihood (Laplace): -414.02
## 
## Point Estimates and 95% Sampling-based CI:
##     variable  Estimate  Std. Error  Lower 95%  Upper 95% 
## Intercept_c    3.43325     0.05232    3.33253    3.53604 
## b[talk]        0.26572     0.05305    0.15585    0.36839 
## b[perf]        0.14054     0.02924    0.08274    0.19857 
## b[talk:perf]   0.13022     0.03063    0.06780    0.19111 
## sigma          0.87135     0.03629    0.80181    0.94375 
## Intercept      3.41638     0.05254    3.31562    3.51660 

係数を一覧したいときは、plot_forest() が便利です。

fit_lm$draws(c("b[talk]", "b[perf]", "b[talk:perf]")) |>
  plot_forest(point_estimate = "MAP")

Regression coefficient forest plot

7. 交互作用を図で確認する

交互作用は、係数表だけでは解釈しにくいことがあります。 その場合は、conditional_effects() を使って予測値を図示します。

ce <- conditional_effects(fit_lm, effect = "talk:perf")
plot(ce)

Conditional effect plot

effect = "talk:perf" と書くと、talk の効果が perf の値によってどのように変わるかを確認できます。

より詳しく調べたい場合は、単純傾斜を計算する simple_effects() も使えます。

simple_effects(fit_lm, effect = "talk:perf")

8. t 検定を頻度主義的に行う

BayesRTMB のラッパー関数は、頻度主義的な分析にも使えます。 prior_flat() を使った t 検定に対して classic() を呼ぶと、通常の t 検定に近い形で結果を表示できます。

mdl_t <- rtmb_ttest(
  sat ~ cond,
  data = debate,
  prior = prior_flat()
)

fit_t_classic <- mdl_t$classic()
fit_t_classic

## Pre-checking model code...
## Checking RTMB setup...
## Starting RTMB optimization...
## 
## 
## Call:
## Classical estimation via ttest 
## 
## Log-Likelihood: -421.320, AIC: 848.640, BIC: 842.640
## 
## Point Estimates and Confidence Intervals:
##            Estimate Std. Error Lower 95% Upper 95%  df  t value     Pr
## diff       -0.37333    0.11297  -0.59564  -0.15102 298 -3.30484 .00107  **
## delta      -0.38161    0.11652  -0.61092  -0.15230 298       NA
## total_mean  3.43333    0.05648   3.32218   3.54449 298       NA
## sd          0.97831    0.04007   0.90254   1.06044 298       NA
## mean0       3.24667    0.07988   3.08947   3.40386 298       NA
## mean1       3.62000    0.07988   3.46280   3.77720 298       NA

ここでは、diff が2群の平均差、delta が標準化効果量を表します。 classic() は、BayesRTMB のモデルを頻度主義的な推定として確認したいときに使えます。

9. JZS prior で Bayes factor を計算する

同じ t 検定でも、JZS prior を使うと、効果量 delta に Cauchy 事前分布を置いた Bayes factor を計算できます。

mdl_t_jzs <- rtmb_ttest(
  sat ~ cond,
  data = debate,
  prior = prior_jzs()
)

set.seed(2)

fit_t_jzs <- mdl_t_jzs$sample()

bf <- fit_t_jzs$bayes_factor(fixed = list(delta = 0))
bf

## --- Bayes Factor Analysis (Bridge Sampling) ---
## Bayes Factor (BF12) : 21.4323
## Log Bayes Factor    : 3.0649 (Approx. Error = 0.0022)
## Evidence            : Strong evidence for Model 1 
## Comparison model    : Parameters fixed at list(delta = 0) 

fixed = list(delta = 0) は、効果量を 0 に固定した帰無モデルと比較する指定です。 この例では、Model 1、つまり効果量を推定するモデルのほうが支持されています。

実際の分析では、Bayes factor の安定性を確認するために、ここで示した例よりも多い MCMC サンプル数を使うことをおすすめします。

次に読むページ

このページでは、BayesRTMB の入口だけを扱いました。 目的に応じて、次のページに進んでください。

1. ラッパー関数の使い方
   回帰、GLM、混合モデル、t 検定、相関、因子分析、IRT など、標準的な分析をラッパー関数で行う方法を確認できます。

2. 階層モデル・GLMM・分散分析
   rtmb_glmer() を使った階層モデル、GLMM、残差相関、条件付き効果の可視化を詳しく確認できます。

3. モデルコードの書き方
   rtmb_code() の setup, parameters, transform, model, generate を使って、独自モデルを書く方法を学べます。

4. RTMB の仕組みと推定アルゴリズム
   MAP 推定、Laplace 近似、MCMC、変分推論などの内部処理を確認できます。