The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

BayesRTMB の概要

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、変分推論などの内部処理を確認できます。

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.