1. 全体の流れ

BayesRTMB の内部処理は、おおまかに次の流れで進みます。

rtmb_code()
  -> rtmb_model()
    -> RTMB_Model
      -> build_ad_obj()
        -> RTMB::MakeADFun()
          -> sample() / optimize() / variational() / classic()
            -> MCMC_Fit / MAP_Fit / VB_Fit / Classic_Fit

ユーザーが書くのは、主に rtmb_code() と rtmb_model() までです。

rtmb_model() は、モデルコードとデータを検査し、RTMB_Model という R6 オブジェクトを作ります。この RTMB_Model が、モデル定義、データ、パラメータ定義、変換ブロック、生成量ブロック、ラッパー関数由来のメタデータを保持します。

その後、推定メソッドに応じて RTMB::MakeADFun() のための自動微分オブジェクトが作られます。得られた結果は、推定方法ごとに次の fit object として返されます。

BayesRTMB は、機能的には MCMC、MAP、VB、classic の順に重心を置いています。MCMC は事後分布を直接扱う中心的な推定法、MAP は高速な点推定と近似的な確認、VB は大きなモデルの近似推論、classic は wrapper 由来モデルを頻度主義的に確認するための補助的な入口です。

2. rtmb_code の各ブロック

rtmb_code() は、モデルをいくつかのブロックに分けて書くための関数です。

code <- rtmb_code(
  data = {
    # 入力データの宣言
  },
  setup = {
    # データ前処理
  },
  parameters = {
    # 推定パラメータの宣言
  },
  transform = {
    # パラメータやデータから作る派生量
  },
  model = {
    # 尤度と事前分布
  },
  generate = {
    # 推定後に計算したい生成量
  }
)

すべてのブロックが必須ではありません。最小限必要なのは parameters と model です。

parameters = {
  alpha <- Dim()
  beta  <- Dim(P)
  sigma <- Dim(lower = 0)
  r     <- Dim(G, random = TRUE)
}

transform = {
  mu <- alpha + X %*% beta
}

model = {
  beta ~ normal(0, 10)
  sigma ~ exponential(1)
  Y ~ normal(mu, sigma)
}

Y ~ normal(mu, sigma)

lp <- lp + normal_lpdf(Y, mu, sigma)

generate = {
  y_rep <- rnorm(length(Y), mu, sigma)
}

3. パラメータ表現

ユーザーは、通常、パラメータを自然な尺度で考えます。

• 標準偏差 sigma は正の値。
• 確率 theta は 0 から 1 の範囲。
• ordered cutpoints は昇順。
• simplex は合計が 1。

一方で、勾配ベースの最適化や HMC は、制約のない実数空間で動かすほうが扱いやすくなります。

そのため BayesRTMB は、内部で次の2つの表現を使い分けます。

たとえば、sigma <- Dim(lower = 0) と宣言した場合、ユーザー向けには正の sigma として表示されますが、内部では実数全体を動く値に変換されます。

この変換は、主に次の helper によって扱われます。

• to_unconstrained()
• to_constrained()
• constrained_vector_to_list()
• unconstrained_vector_to_list()
• generate_flat_names()

この層は、fixed、map、ランダム効果、Laplace 近似、summary 表示の整合性に深く関わります。

4. ヤコビアン補正

制約付きのパラメータを非制約空間へ変換すると、確率密度の尺度も変わります。そのため、事後分布を正しく扱うにはヤコビアン補正が必要になります。

概念的には、自然尺度のパラメータを y、非制約尺度のパラメータを x とすると、次の関係を考えます。

\[ \log p_X(x) = \log p_Y(y) + \log \left| \frac{dy}{dx} \right| \]

BayesRTMB では、ユーザーは自然尺度でモデルを書きます。内部では to_unconstrained()、to_constrained()、calc_log_jacobian() を使い、推定方法に応じて必要なヤコビアン補正を目的関数に反映します。

補正の対象は jacobian_target によって制御されます。

用途としては、MCMC では非制約空間で事後分布全体をサンプリングするため、基本的に jacobian_target = "all" を使います。

一方、optimize() で Laplace 近似を使う場合は、Laplace で積分する random 側の密度変換が必要になるため、基本的に jacobian_target = "random" を使います。Laplace 近似を使わない通常の最適化では、内部目的関数上では jacobian_target = "none" が使われます。

通常のユーザーがこの設定を直接操作する必要はほとんどありません。ただし、profile likelihood や Laplace 近似を含む内部処理では、この区別が重要になります。

5. RTMB の自動微分オブジェクト

RTMB は、TMB の自動微分エンジンを R から使うためのパッケージです。

BayesRTMB では、RTMB_Model$build_ad_obj() が RTMB::MakeADFun() を呼び出し、関数値、勾配、ヘッセ行列などを評価できるオブジェクトを作ります。

このとき、モデル関数は通常の数値ではなく AD 型の値を使って評価されます。四則演算、log()、exp()、行列演算、分布関数などの計算過程が、自動微分のための計算グラフとして記録されます。

Stan のようにモデルコードを C++ に書き換えて明示的にコンパイルするのではなく、R で書いたモデル関数を AD 型で評価し、MakeADFun() object として高速に関数値や導関数を評価できるようにする、というイメージです。

rtmb_model() はモデル作成時に pre-check を行います。

• setup や parameters に未定義変数がないか。
• model や transform が初期値で実行できるか。
• log_prob がスカラーの数値を返すか。
• MakeADFun() で勾配が計算できるか。

この段階で構造的なエラーを早めに検出することで、RTMB/TMB 側の分かりにくいエラーをなるべく避ける設計になっています。

6. 推定メソッドごとの違い

fit <- mdl$sample()

fit <- mdl$optimize()

fit <- mdl$variational()

fit <- mdl$classic()

7. optimize と classic の推定結果の違い

BayesRTMB の設計では、optimize() と classic() の責務を明確に分けています。

optimize() は、MAP / ML / marginal MAP のための API です。informative prior を含むモデルも扱います。

classic() は、wrapper 由来かつ flat-prior のモデルだけを対象にする頻度主義的 API です。内部的には固定効果をLaplace近似で周辺化した制約付き最尤法（REML）を行うことで、モーメント法と一致する推定量を得ています。

低水準では、どちらも RTMB による最適化エンジンを共有します。内部の .fit_rtmb() は、MakeADFun() object を作り、最適化を行い、sdreport() や共分散などの raw components を返します。

ただし、.fit_rtmb() 自体は fit object を作りません。

この分離により、同じ最適化エンジンを共有しながら、MAP 推定と頻度主義的分析の表示、自由度、prior correction、後処理を分けています。

8. Laplace 近似と random

階層モデルや潜在変数モデルでは、すべてのパラメータを同時に固定効果のように扱うよりも、一部を積分消去したいことがあります。

BayesRTMB では、Dim(..., random = TRUE) と宣言したパラメータは、Laplace 近似の対象になりえます。

parameters = {
  alpha <- Dim()
  tau   <- Dim(lower = 0)
  r     <- Dim(G, random = TRUE)
}

optimize(laplace = TRUE) では、random = TRUE のパラメータを RTMB::MakeADFun(random = ...) に渡し、Laplace 近似で周辺化します。optimizeではlaplace = TRUEがデフォルトです。

また、optimize(marginal = ...) を使うと、parameters ブロックで宣言された primary parameter を一時的に random として扱い、marginal MAP / empirical Bayes 型の推定を行えます。flat priorにおいて固定効果をmarginalに指定すれば、制約付き最尤法(REML)と一致します。

fit <- mdl$optimize(marginal = c("mean", "delta"))

ここで重要なのが、marginal_vars と laplace_random_vars の区別です。

mixed model では、もともとの random effects と、marginal で指定された primary parameter の両方が Laplace 対象になることがあります(REML) 。そのため、laplace_random_vars は marginal_vars より広い集合になることがあります。

marginal = "auto" は、wrapper が extra$marginal に保存した primary parameter 名を使います。ここに入れられるのは、parameters ブロックで宣言されたパラメータだけです。transform や generate で作られた派生量は指定できません。optimizeでREML推定をしたいときはラッパー関数を使ってmarginal = "auto"とすれば簡単に実行できます。また、そのときは自動的に自由度推定がdf_method = "satterthwaite"となり、t分布で信頼区間を計算します。df_method = Infとすれば正規分布を仮定します。

9. 標準誤差と区間推定

optimize() と classic() では、点推定だけでなく標準誤差や区間推定も扱います。

optimize() の se_method には、主に次の選択肢があります。

se_method = "none" は軽量モードです。標準誤差を 0 とみなすのではなく、標準誤差や信頼区間を計算しない、という意味です。

transform で作られた派生量は、Wald 型ではデルタ法、se_method = "sampling" ではシミュレーションにもとづいて不確実性を伝播します。

generate は基本的に推定後の生成量です。生成量の標準誤差や区間を必要とする場合は、MCMC サンプルや generated_quantities() の結果を使って確認するのが自然です。

sdreport() が信頼できる標準誤差を返せない場合、BayesRTMB は診断メッセージを出します。たとえば、pdHess = FALSE、cov.fixed が返らない、標準誤差が非有限、Hessian fallback、MASS::ginv() fallback などが該当します。

10. AD 型のコードを書くときの注意点

RTMB では、モデルコードは通常の R のように見えますが、内部では自動微分の計算グラフとして記録されます。そのため、普通の R では問題なくても、AD 型の計算では避けたほうがよい書き方があります。

# 避けたい例
if (sigma > 1) {
  lp <- lp + normal_lpdf(y, 0, sigma)
} else {
  lp <- lp + normal_lpdf(y, 0, 1)
}

11. Stan、TMB、RTMB との関係

BayesRTMB は、RTMB をバックエンドとして使います。RTMB は、TMB の自動微分エンジンを R から使えるようにしたパッケージです。

Stan と似ている点は、ユーザーが確率モデルを書き、制約変換や勾配ベースの推定を使えることです。

一方で、BayesRTMB では R の構文でモデルを書き、rtmb_code() のブロックを通じて RTMB の目的関数を作ります。ラッパー関数が生成したモデルも print_code() で確認できるため、標準的な分析から独自モデルへ進みやすい設計になっています。また、テープの記録ではC++へコンパイルをする必要がないので、スピーディーに分析が可能です。

12. まとめ

BayesRTMB の内部構造は、次のように整理できます。

• rtmb_code() は、モデルをブロック構造として記述する入口です。
• rtmb_model() は、データとコードを検査し、RTMB_Model を作ります。
• RTMB_Model は、モデル、パラメータ、変換、生成量、wrapper metadata を保持します。
• build_ad_obj() は、RTMB::MakeADFun() のための自動微分オブジェクトを作ります。
• sample()、optimize()、variational()、classic() は、同じモデル表現から異なる推定を実行します。
• optimize() と classic() は低水準エンジンを共有しますが、責務と返り値は明確に分かれています。
• Laplace 近似では、random = TRUE と marginal の扱い、特に marginal_vars と laplace_random_vars の区別が重要です。

より実践的なモデルの書き方は、モデルコードの書き方 を参照してください。ラッパー関数がどのような rtmb_code() を作るかは、ラッパー関数の使い方 と print_code() が参考になります。階層モデルや GLMM の具体的な使い方は、階層モデル・GLMM で確認できます。