ラッパー関数とは

ラッパー関数は、典型的な統計モデルを簡潔に書くための入り口です。たとえば、重回帰なら次のように書けます。

library(BayesRTMB)
data(debate)

mdl <- rtmb_lm(sat ~ talk * perf, data = debate)

この時点では、まだ推定は実行されていません。モデルを作ったあとに、目的に応じて推定方法を選びます。

fit_mcmc <- mdl$sample()
fit_map  <- mdl$optimize()
fit_cl   <- mdl$classic()

同じ mdl から、ベイズ推定と頻度主義的分析の両方を扱えるところが、BayesRTMB のラッパー関数の大きな特徴です。

ラッパー関数一覧

主なラッパー関数は次のとおりです。ここでは、細かなオプションよりも典型的な使い方を示します。

推定のイメージ

ラッパー関数でモデルを作ったあとは、メソッドを選んで推定します。

MCMC を使う場合は sample() です。

mdl <- rtmb_lm(sat ~ talk * perf, data = debate)
fit <- mdl$sample()
fit$summary()

MAP推定を使う場合は optimize() です。

fit <- mdl$optimize()
fit$summary()

頻度主義的分析を使う場合は classic() です。

fit <- mdl$classic()
fit$summary()

このように、モデルの書き方は同じで、推定方法だけを切り替えられます。

頻度主義的分析

classic() は、BayesRTMB のラッパー関数を頻度主義的な分析として使うためのメソッドです。

rtmb_corr(cbind(sat, perf), data = debate)$classic()

rtmb_corr() のデフォルトは method = "pearson" なので、classic() は cor.test() と同じ Pearson 相関の結果を返します。

## Call:
## Classical estimation via corr 
## 
## Point Estimates and Confidence Intervals:
##           Estimate Std. Error Lower 95% Upper 95%  df t value     Pr    
## corr[rho]  0.29742         NA   0.19059   0.39728 298 5.37755 <.0001 ***

統制変数を入れた場合は、共変量だけを統制した偏相関になります。

rtmb_corr(cbind(sat, perf),
          data = debate,
          covariates = ~ skill)$classic()

##           Estimate Std. Error Lower 95% Upper 95%  df t value     Pr    
## pcorr[rho]  0.29604         NA   0.18896   0.39617 297 5.34134 <.0001 ***

method = "spearman" を指定すると Spearman 相関を使います。method = "reml" を指定すると、RTMB のREML推定に基づく相関モデルとして推定されます。

rtmb_corr(cbind(sat, perf), data = debate, method = "spearman")$classic()
rtmb_corr(cbind(sat, perf), data = debate, method = "reml")$classic()

fit_tab <- rtmb_table(skill, cond, data = debate)$classic()
anova(fit_tab)

## Contingency Table Analysis Tests
## - 
##                                        X-squared df p-value
## Pearson's Chi-squared (from est. prob)   1.20478  2  0.5475
## 
## Fisher's Exact Test for Count Data
## p-value =  0.54310 

data(BigFive)

items <- BigFive[, 1:10]
fit_fa1 <- rtmb_fa(items, nfactors = 1)$classic()
fit_fa2 <- rtmb_fa(items, nfactors = 2)$classic()

anova(fit_fa1, fit_fa2)

## Likelihood Ratio Tests
##    Model Df  logLik    AIC    BIC  Chisq Chi Df Pr(>Chisq)
## fit_fa1 30 -2485.5 5030.9 4970.9     NA     NA         NA
## fit_fa2 40 -2440.4 4960.9 4880.9 90.089     10 5.1421e-15

c(AIC_1 = AIC(fit_fa1),
  AIC_2 = AIC(fit_fa2),
  BIC_1 = BIC(fit_fa1),
  BIC_2 = BIC(fit_fa2))

##    AIC_1    AIC_2    BIC_1    BIC_2 
## 5030.942 4960.853 4970.942 4880.853

print_code()

ラッパー関数は簡単に使えますが、内部ではすべて rtmb_code() に展開されています。どのようなモデルが生成されたかは print_code() で確認できます。

mdl <- rtmb_ttest(sat ~ cond, data = debate)
mdl$print_code()

## === RTMB Model Code ===
## 
## rtmb_code(
##   setup = {
##     Y1 <- as.numeric(na.omit(Y[G == 1]))
##     Y2 <- as.numeric(na.omit(Y[G == 2]))
##   }, 
##   parameters = {
##     mean0 = Dim(1)
##     mean1 = Dim(1)
##     sd = Dim(2, lower = 0)
##   }, 
##   transform = {
##     diff <- mean0 - mean1
##     mean <- (mean0 + mean1)/2
##     sd_pooled <- sqrt((sd[1]^2 + sd[2]^2)/2)
##     delta <- diff/sd_pooled
##   }, 
##   model = {
##     Y1 ~ normal(mean0, sd[1])
##     Y2 ~ normal(mean1, sd[2])
##   }
## )

この例は Welch の t検定に対応するモデルです。通常の検定としては短く呼び出せますが、内部では2群それぞれの平均と標準偏差を推定するモデルとして表現されています。

このように、ラッパー関数は「簡単に使う」ための入り口であると同時に、「どのモデルを推定しているか」を確認しながら、必要に応じて独自の rtmb_code() へ進むための足場にもなります。

欠損値の扱い (Missing Data Handling)

BayesRTMB のラッパー関数は、欠損値（NA）を適切かつ統一的に扱うための missing 引数を備えています。分析手法やデータの形式に応じて、最適な手法を選択できます。

多くのラッパー関数ではデフォルトが "listwise"（リストワイズ削除）となっていますが、多変量解析を扱う一部のラッパー関数では、完全情報最尤法（FIML） や ペアワイズ削除 も選択可能です。

data(BigFive)
# 意図的にランダムな欠損を生成
set.seed(42)
BigFive_miss <- BigFive
for (col in 1:ncol(BigFive_miss)) {
  BigFive_miss[sample(1:nrow(BigFive_miss), 10), col] <- NA
}

# FIMLによる探索的因子分析
mdl_fa <- rtmb_fa(BigFive_miss, nfactors = 5, missing = "fiml")
fit_fa <- mdl_fa$optimize()

# ペアワイズ削除による古典的相関分析
mdl_corr <- rtmb_corr(BigFive_miss[, 1:5], method = "pearson", missing = "pairwise")
fit_corr <- mdl_corr$classic()
fit_corr$summary()

次に読む記事

BayesRTMB の全体像や、より細かな使い方を知りたい場合は、次の記事も参照してください。

• イントロダクション: BayesRTMB の基本的な考え方と、MCMC、MAP推定、変分推論、頻度主義的分析の位置づけ。
• クイックスタート: インストール、最初のモデル、MCMC診断、可視化、t検定のベイズファクター。
• 階層モデル・GLMM・分散分析: rtmb_glmer() を使って mixed model / GLMM を詳しく扱う方法。
• モデルを書く: rtmb_code() を使って独自のモデルを書く方法。
• 内部構造: setup、parameters、transform、model など、BayesRTMB の内部的なモデル構造。