BayesRTMB とは

BayesRTMB は、RTMB を計算エンジンとして用いながら、R の中で統計モデルを記述し、推定するためのパッケージです。

BayesRTMB では、標準的な分析は rtmb_lm() や rtmb_glmer() などのラッパー関数から始められます。一方で、必要に応じて rtmb_code() を使い、Stan に近い感覚で独自のモデルを書くこともできます。

1つのモデルオブジェクトから、MCMC、MAP 推定、変分推論、そして頻度主義的分析を呼び出せる点が BayesRTMB の基本的な考え方です。

fit_mcmc <- mdl$sample()
fit_map <- mdl$optimize()
fit_vb <- mdl$variational()
fit_freq <- mdl$classic()

このページでは、BayesRTMB の位置づけと、パッケージを使ううえで最初に知っておくとよい全体像を紹介します。詳しいコードの書き方や個別の分析例は、他の vignette で扱います。

記事へのリンク

以下の記事があります。

BayesRTMB でできること

BayesRTMB の主な特徴は次の通りです。

R の中でモデルを書ける
rtmb_code() を使い、parameters, transform, model, generate などのブロックに分けてモデルを記述できます。
ラッパー関数からすぐに分析できる
回帰、一般化線形モデル、混合モデル、t 検定、相関、因子分析、IRT、潜在ランク理論、多次元展開法などを、専用のラッパー関数から実行できます。
同じモデルから複数の推定方法を使える
MCMC、MAP 推定、変分推論、頻度主義的分析を、同じ RTMB_Model オブジェクトから呼び出せます。
ランダム効果を効率的に扱える
random = TRUE と宣言したパラメータは、MAP 推定では Laplace 近似によって周辺化できます。
wrapper で作られたモデルを確認できる
print_code() を使うと、ラッパー関数が内部で作っている rtmb_code() を確認できます。
モデル比較にも対応している
MCMC の結果から bridge sampling による周辺尤度や Bayes factor を計算できます。

2つの入口

BayesRTMB には、大きく分けて2つの入口があります。

1つ目は、標準的な分析を行うための ラッパー関数 です。通常はこちらから始めるのが簡単です。

data(debate)

mdl <- rtmb_lm(sat ~ talk + perf, data = debate)
fit <- mdl$optimize()
fit

2つ目は、モデルを自分で書くための rtmb_code() です。ラッパー関数では表現しにくいモデルや、新しいモデルを試したい場合に使います。

dat <- list(
  Y = debate$sat,
  X = data.matrix(debate[c("talk", "perf")])
)

code <- rtmb_code(
  setup = {
    N <- length(Y)
    P <- ncol(X)
  },
  parameters = {
    alpha = Dim()
    beta = Dim(P)
    sigma = Dim(lower = 0)
  },
  transform = {
    mu <- alpha + X %*% beta
  },
  model = {
    Y ~ normal(mu, sigma)
    alpha ~ normal(0, 10)
    beta ~ normal(0, 10)
    sigma ~ exponential(1)
  }
)

mdl <- rtmb_model(dat, code)

model ブロックでは、次のような sampling syntax を使います。

Y ~ normal(mu, sigma)
beta ~ normal(0, 10)
sigma ~ exponential(1)

これは内部的には対数密度を足し合わせるコードに変換されます。

ラッパー関数から始める

まずはラッパー関数を使った例を見てみます。ここでは debate データを使い、満足度 sat を talk と perf で説明する線形回帰モデルを考えます。

data(debate)

mdl_lm <- rtmb_lm(
  sat ~ talk + perf,
  data = debate
)

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

MAP 推定を行うには、optimize() を使います。

fit_map <- mdl_lm$optimize()
fit_map

## Starting RTMB optimization...
## 
## 
## Call:
## MAP Estimation via RTMB
## 
## Negative Log-Posterior: 395.58
## Approx. Log Marginal Likelihood (Laplace): -404.61
## 
## Point Estimates and 95% Wald CI:
##  variable  Estimate  Std. Error  Lower 95%  Upper 95% 
## Intercept   1.83366     0.21105    1.42000    2.24732 
## b[talk]     0.28694     0.05291    0.18323    0.39064 
## b[perf]     0.15632     0.02987    0.09777    0.21486 
## sigma       0.90453     0.03693    0.83497    0.97988

頻度主義的な線形回帰として推定したい場合は、classic() を使います。classic() は、ラッパー関数で作られた flat prior のモデルに対して利用できます。

fit_freq <- mdl_lm$classic()
fit_freq

## Starting RTMB optimization...
## 
## 
## Call:
## Classical estimation via lm 
## 
## Log-Likelihood: -402.220, AIC: 812.440, BIC: 827.255
## 
## Point Estimates and Confidence Intervals:
##           Estimate Std. Error Lower 95% Upper 95%  df  t value     Pr    
## Intercept  1.83366    0.21212   1.41621   2.25110 297  8.64454 <.0001 ***
## b[talk]    0.28694    0.05318   0.18229   0.39159 297  5.39580 <.0001 ***
## b[perf]    0.15632    0.03002   0.09724   0.21539 297  5.20703 <.0001 ***
## sigma      0.90909    0.03730   0.83856   0.98554 297     NA

ラッパー関数がどのようなモデルコードを作っているかは、print_code() で確認できます。

mdl_lm$print_code()

## === RTMB Model Code ===
## 
## rtmb_code(
##   setup = {
##     mf <- model.frame(formula, df)
##     Y <- model.response(mf)
##     X_full <- model.matrix(formula, mf)
##     X <- X_full[, colnames(X_full) != "(Intercept)", drop = FALSE]
##     N <- length(Y)
##     K <- ncol(X)
##   }, 
##   parameters = {
##     Intercept <- Dim(1)
##     b <- Dim(K)
##     sigma <- Dim(num_sigma_groups, lower = 0)
##   }, 
##   model = {
##     # Transform
##     eta <- Intercept + X %*% b
##     # Likelihood (Data)
##     Y ~ normal(eta, sigma)
##     # Priors
##   }
## )

print_code() は、ユーザーが rtmb_model() で同じモデルを再現できるように、setup, parameters, transform, model などのブロックを表示します。

自分でモデルを書く

より柔軟なモデルを書きたい場合は、rtmb_code() を使います。

rtmb_code() の主なブロックは次の通りです。

setup: データから必要な量を前処理する
parameters: 推定するパラメータを宣言する
transform: パラメータから派生量や線形予測子を作る
model: 尤度と事前分布を書く
generate: 予測値や生成量を計算する

たとえば、prior_normal() に近い素直な normal / exponential prior を持つ回帰モデルは次のように書けます。

dat <- list(
  Y = as.numeric(debate$sat),
  X = data.matrix(debate[c("talk", "perf", "skill")])
)

code_reg <- rtmb_code(
  setup = {
    N <- length(Y)
    P <- ncol(X)
  },
  parameters = {
    alpha = Dim()
    beta = Dim(P)
    sigma = Dim(lower = 0)
  },
  transform = {
    mu <- alpha + X %*% beta
  },
  model = {
    Y ~ normal(mu, sigma)
    alpha ~ normal(0, 10)
    beta ~ normal(0, 2.5)
    sigma ~ exponential(1)
  },
  generate = {
    Y_rep <- rnorm(length(Y), mu, sigma)
  }
)

mdl_reg <- rtmb_model(dat, code_reg)

setup に前処理を書くことで、parameters や model ブロックを読みやすくできます。また、print_code() で表示したときにも、モデルの再現に必要な処理が見えるようになります。

推定方法の使い分け

BayesRTMB は Bayes 推定を基本にしています。そのため、ここでは MCMC、MAP 推定、変分推論、そして頻度主義的分析の順に整理します。

メソッド	主な目的
`$sample()`	NUTS による MCMC サンプリングを行う
`$optimize()`	MAP 推定または最尤推定を高速に行う
`$variational()`	ADVI による近似事後分布を求める
`$classic()`	flat prior のラッパーモデルを頻度主義的分析として扱う

MCMC

sample() は、NUTS による MCMC サンプリングを行います。事後分布全体を見たい場合や、MAP 推定だけでは不確実性を十分に表現しにくい場合に使います。

set.seed(1)

fit_mcmc <- mdl_lm$sample()

fit_mcmc$summary()

##  variable     mean    sd      map     q2.5    q97.5  ess_bulk  ess_tail  rhat 
## lp         -397.79  1.54  -396.75  -401.63  -395.94       859      1236  1.01
## Intercept     1.82  0.22     1.78     1.38     2.26       629      1023  1.01
## b[talk]       0.29  0.05     0.28     0.18     0.40       665      1205  1.00
## b[perf]       0.16  0.03     0.15     0.10     0.22       671       990  1.01
## sigma         0.91  0.04     0.91     0.84     1.00       934      1254  1.01

MAP 推定

optimize() は、事後分布または尤度を最大化する点推定を行います。モデルの確認や、複雑なモデルの初期検討に向いています。

fit_map <- mdl_lm$optimize()
fit_map$summary()

## Starting RTMB optimization...
## 
## 
## Call:
## MAP Estimation via RTMB
## 
## Negative Log-Posterior: 395.58
## Approx. Log Marginal Likelihood (Laplace): -404.61
## 
## Point Estimates and 95% Wald CI:
##  variable  Estimate  Std. Error  Lower 95%  Upper 95% 
## Intercept   1.83366     0.21105    1.42000    2.24732 
## b[talk]     0.28694     0.05291    0.18323    0.39064 
## b[perf]     0.15632     0.02987    0.09777    0.21486 
## sigma       0.90453     0.03693    0.83497    0.97988

se_method = "sampling" を使うと、近似正規サンプルによる不確実性の伝播を使って、派生量の標準誤差や区間を計算できます。

fit_map <- mdl_lm$optimize(se_method = "sampling")

変分推論

variational() は、ADVI によって事後分布を近似します。MCMC より高速におおまかな結果を得たい場合に使います。ただし、事後分布の不確実性を過小評価することがあるため、最終的な不確実性評価には注意が必要です。

fit_vb <- mdl_lm$variational(
  method = "meanfield",
  iter = 3000,
  num_estimate = 4
)

頻度主義的分析

classic() は、ラッパー関数で作られた flat prior のモデルに対して、頻度主義的な推定結果を返します。

fit_freq <- mdl_lm$classic()
fit_freq$summary()

## Starting RTMB optimization...
## 
## 
## Call:
## Classical estimation via lm 
## 
## Log-Likelihood: -402.220, AIC: 812.440, BIC: 827.255
## 
## Point Estimates and Confidence Intervals:
##           Estimate Std. Error Lower 95% Upper 95%  df  t value     Pr    
## Intercept  1.83366    0.21212   1.41621   2.25110 297  8.64454 <.0001 ***
## b[talk]    0.28694    0.05318   0.18229   0.39159 297  5.39580 <.0001 ***
## b[perf]    0.15632    0.03002   0.09724   0.21539 297  5.20703 <.0001 ***
## sigma      0.90909    0.03730   0.83856   0.98554 297     NA

lm や glm 型のモデルでは、検定統計量は通常の表記に近い形で表示されます。混合モデルでは、必要に応じて自由度補正や Laplace 近似が使われます。

ランダム効果と Laplace 近似

階層モデルや混合モデルでは、パラメータを random = TRUE として宣言できます。

Y <- debate$sat
group <- as.integer(factor(debate$group))
G <- length(unique(group))

dat_icc <- list(Y = Y, group = group, G = G)

code_icc <- rtmb_code(
  parameters = {
    mu = Dim()
    sigma = Dim(lower = 0)
    tau = Dim(lower = 0)
    r = Dim(G, random = TRUE)
  },
  model = {
    Y ~ normal(mu + r[group] * tau, sigma)
    r ~ normal(0, 1)
    tau ~ exponential(1)
    sigma ~ exponential(1)
  },
  generate = {
    icc <- tau^2 / (tau^2 + sigma^2)
  }
)

mdl_icc <- rtmb_model(dat_icc, code_icc)
fit_icc <- mdl_icc$optimize(laplace = TRUE)

MAP 推定では、laplace = TRUE が既定値です。random = TRUE とされたパラメータは、Laplace 近似によって周辺化されます。

ラッパー関数を使う場合は、次のように混合モデルを書けます。

mdl_hlm <- rtmb_glmer(
  sat ~ talk + perf + (1 | group),
  data = debate,
  family = "gaussian"
)

fit_hlm <- mdl_hlm$optimize()

モデル比較（bridge sampling / WAIC）

bridge sampling による周辺尤度と Bayes factor

MCMC の結果に対して、bridge sampling による周辺尤度を計算できます。

set.seed(1)

fit_mcmc$bridgesampling()

表示は、たとえば次のように出力されます。

## Bridge Sampling Converged: LogML = -404.591 (Error = 0.0032, ESS = 604.0)

戻り値は数値の log marginal likelihood で、error と ess の属性を持ちます。

Bayes factor を計算したい場合は、bayes_factor() を使います。

bf <- fit_mcmc$bayes_factor(fixed = list("b[talk]" = 0))
bf

## --- Bayes Factor Analysis (Bridge Sampling) ---
## Bayes Factor (BF12) : 142963.4
## Log Bayes Factor    : 11.8703 (Approx. Error = 0.0040)
## Evidence            : Decisive evidence for Model 1 
## Comparison model    : Parameters fixed at list("b[talk]" = 0)

WAIC によるモデル比較

BayesRTMB では、rtmb_code() の generate ブロックで、観測ごとの対数尤度を log_lik という変数に格納しておくと、推定後に $WAIC() メソッドを使って WAIC（広く使える情報量基準）を計算できます。

たとえば、手書きモデルの generate ブロックに以下のように log_lik を定義します。

  generate = {
    # 観測ごとの対数尤度を計算して log_lik に代入する（sum = FALSE を指定）
    log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE)
  }

※ ラッパー関数（rtmb_lm() など）を使う場合は、WAIC = TRUE を引数に指定することで、自動的に内部で log_lik が生成されます。

mdl_lm <- rtmb_lm(sat ~ talk + perf, data = debate, WAIC = TRUE)
fit_mcmc <- mdl_lm$sample()

# WAIC の計算
fit_mcmc$WAIC()

計算結果は以下のように出力されます。

## WAIC
##   WAIC (elpd scale): -400.07398
##   -2 WAIC: 800.14796
##   p_waic: 4.70921
##   lppd: -395.36477
##   nobs: 300, draws: 4000

WAIC() メソッドは、MCMC (MCMC_Fit) や変分推論 (VB_Fit)、および se_method = "sampling" を指定した MAP 推定 (MAP_Fit) の結果に対して利用できます。

次のステップ

BayesRTMB の全体像がつかめたら、次は目的に応じて以下の記事を読むのがおすすめです。

クイックスタート
基本的なモデルを実際に動かしながら、rtmb_code() と rtmb_model() の流れを確認できます。
ラッパー関数の使い方
回帰、一般化線形モデル、混合モデル、t 検定など、標準的な分析をラッパー関数で実行する方法を確認できます。
階層モデル・GLMM・分散分析
rtmb_glmer() を使った階層モデル、GLMM、順序モデル、残差相関、可視化の流れを確認できます。
モデルコードの書き方
setup, parameters, transform, model, generate の各ブロックを使って、自分でモデルを書く方法を詳しく説明します。
RTMB の仕組みと推定アルゴリズム
Laplace 近似、制約付きパラメータ、MCMC、変分推論などの内部処理を確認できます。

BayesRTMB の概要