--- title: "BayesRTMB の概要" description: "BayesRTMB の考え方、モデルの作り方、推定方法の全体像を紹介します。" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{BayesRTMB の概要} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, eval = FALSE) ``` # BayesRTMB とは **BayesRTMB** は、RTMB を計算エンジンとして用いながら、R の中で統計モデルを記述し、推定するためのパッケージです。 BayesRTMB では、標準的な分析は `rtmb_lm()` や `rtmb_glmer()` などのラッパー関数から始められます。一方で、必要に応じて `rtmb_code()` を使い、Stan に近い感覚で独自のモデルを書くこともできます。 1つのモデルオブジェクトから、MCMC、MAP 推定、変分推論、そして頻度主義的分析を呼び出せる点が BayesRTMB の基本的な考え方です。 ```{r eval=FALSE} fit_mcmc <- mdl$sample() fit_map <- mdl$optimize() fit_vb <- mdl$variational() fit_freq <- mdl$classic() ``` このページでは、BayesRTMB の位置づけと、パッケージを使ううえで最初に知っておくとよい全体像を紹介します。詳しいコードの書き方や個別の分析例は、他の vignette で扱います。 ## 記事へのリンク 以下の記事があります。 - **[日本語版の概要](ja-introduction.html)** - **[クイックスタート](ja-quick_start.html)** - **[ラッパー関数の使い方](ja-wrapper_functions.html)** - **[階層モデル・GLMM・分散分析](ja-rtmb_glmer.html)** - **[モデルコードの書き方](ja-writing_models.html)** - **[RTMB の仕組みと推定アルゴリズム](ja-rtmb_internals.html)** # 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つ目は、標準的な分析を行うための **ラッパー関数** です。通常はこちらから始めるのが簡単です。 ```{r eval=FALSE} data(debate) mdl <- rtmb_lm(sat ~ talk + perf, data = debate) fit <- mdl$optimize() fit ``` 2つ目は、モデルを自分で書くための **`rtmb_code()`** です。ラッパー関数では表現しにくいモデルや、新しいモデルを試したい場合に使います。 ```{r eval=FALSE} 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 を使います。 ```{r eval=FALSE} Y ~ normal(mu, sigma) beta ~ normal(0, 10) sigma ~ exponential(1) ``` これは内部的には対数密度を足し合わせるコードに変換されます。 # ラッパー関数から始める まずはラッパー関数を使った例を見てみます。ここでは `debate` データを使い、満足度 `sat` を `talk` と `perf` で説明する線形回帰モデルを考えます。 ```{r eval=FALSE} data(debate) mdl_lm <- rtmb_lm( sat ~ talk + perf, data = debate ) ``` この時点では、まだ推定は実行されていません。`mdl_lm` は、モデル定義とデータを保持した `RTMB_Model` オブジェクトです。 MAP 推定を行うには、`optimize()` を使います。 ```{r eval=FALSE} fit_map <- mdl_lm$optimize() fit_map ``` ```text ## 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 のモデルに対して利用できます。 ```{r eval=FALSE} fit_freq <- mdl_lm$classic() fit_freq ``` ```text ## 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()` で確認できます。 ```{r eval=FALSE} mdl_lm$print_code() ``` ```text ## === 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 を持つ回帰モデルは次のように書けます。 ```{r eval=FALSE} 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 推定だけでは不確実性を十分に表現しにくい場合に使います。 ```{r eval=FALSE} set.seed(1) fit_mcmc <- mdl_lm$sample() fit_mcmc$summary() ``` ```text ## 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()` は、事後分布または尤度を最大化する点推定を行います。モデルの確認や、複雑なモデルの初期検討に向いています。 ```{r eval=FALSE} fit_map <- mdl_lm$optimize() fit_map$summary() ``` ```text ## 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"` を使うと、近似正規サンプルによる不確実性の伝播を使って、派生量の標準誤差や区間を計算できます。 ```{r eval=FALSE} fit_map <- mdl_lm$optimize(se_method = "sampling") ``` ## 変分推論 `variational()` は、ADVI によって事後分布を近似します。MCMC より高速におおまかな結果を得たい場合に使います。ただし、事後分布の不確実性を過小評価することがあるため、最終的な不確実性評価には注意が必要です。 ```{r eval=FALSE} fit_vb <- mdl_lm$variational( method = "meanfield", iter = 3000, num_estimate = 4 ) ``` ## 頻度主義的分析 `classic()` は、ラッパー関数で作られた flat prior のモデルに対して、頻度主義的な推定結果を返します。 ```{r eval=FALSE} fit_freq <- mdl_lm$classic() fit_freq$summary() ``` ```text ## 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` として宣言できます。 ```{r eval=FALSE} 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 近似によって周辺化されます。 ラッパー関数を使う場合は、次のように混合モデルを書けます。 ```{r eval=FALSE} 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 による周辺尤度を計算できます。 ```{r eval=FALSE} set.seed(1) fit_mcmc$bridgesampling() ``` 表示は、たとえば次のように出力されます。 ```text ## Bridge Sampling Converged: LogML = -404.591 (Error = 0.0032, ESS = 604.0) ``` 戻り値は数値の log marginal likelihood で、`error` と `ess` の属性を持ちます。 Bayes factor を計算したい場合は、`bayes_factor()` を使います。 ```{r eval=FALSE} bf <- fit_mcmc$bayes_factor(fixed = list("b[talk]" = 0)) bf ``` ```text ## --- 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` を定義します。 ```{r eval=FALSE} generate = { # 観測ごとの対数尤度を計算して log_lik に代入する（sum = FALSE を指定） log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE) } ``` ※ ラッパー関数（`rtmb_lm()` など）を使う場合は、`WAIC = TRUE` を引数に指定することで、自動的に内部で `log_lik` が生成されます。 ```{r eval=FALSE} mdl_lm <- rtmb_lm(sat ~ talk + perf, data = debate, WAIC = TRUE) fit_mcmc <- mdl_lm$sample() # WAIC の計算 fit_mcmc$WAIC() ``` 計算結果は以下のように出力されます。 ```text ## 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 の全体像がつかめたら、次は目的に応じて以下の記事を読むのがおすすめです。 1. **[クイックスタート](ja-quick_start.html)** 基本的なモデルを実際に動かしながら、`rtmb_code()` と `rtmb_model()` の流れを確認できます。 2. **[ラッパー関数の使い方](ja-wrapper_functions.html)** 回帰、一般化線形モデル、混合モデル、t 検定など、標準的な分析をラッパー関数で実行する方法を確認できます。 3. **[階層モデル・GLMM・分散分析](ja-rtmb_glmer.html)** `rtmb_glmer()` を使った階層モデル、GLMM、順序モデル、残差相関、可視化の流れを確認できます。 4. **[モデルコードの書き方](ja-writing_models.html)** `setup`, `parameters`, `transform`, `model`, `generate` の各ブロックを使って、自分でモデルを書く方法を詳しく説明します。 5. **[RTMB の仕組みと推定アルゴリズム](ja-rtmb_internals.html)** Laplace 近似、制約付きパラメータ、MCMC、変分推論などの内部処理を確認できます。