---
title: "ラッパー関数の使い方"
description: "BayesRTMB のラッパー関数で、MCMC、MAP推定、変分推論、頻度主義的分析を手早く実行する方法を紹介します。"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{ラッパー関数の使い方}
  %\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 では、`rtmb_code()` を使ってモデルを直接書くこともできますが、よく使う分析についてはラッパー関数が用意されています。

ラッパー関数は、`lm()` や `glm()` のような短い書き方から、内部的に一貫した `rtmb_code()` を生成します。そのため、同じモデルに対して MCMC、MAP推定、変分推論、そして頻度主義的分析を切り替えて使うことができます。

このページでは、各ラッパー関数の細かい理論やオプションではなく、「どのような雰囲気で使うか」を中心に紹介します。

# ラッパー関数とは

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

```{r}
library(BayesRTMB)
data(debate)

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

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

```{r}
fit_mcmc <- mdl$sample()
fit_map  <- mdl$optimize()
fit_cl   <- mdl$classic()
```

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

# ラッパー関数一覧

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

| 関数 | 分析 | 典型的なコード |
| :--- | :--- | :--- |
| `rtmb_ttest()` | t検定 | `rtmb_ttest(sat ~ cond, data = debate)` |
| `rtmb_lm()` | 線形回帰 | `rtmb_lm(sat ~ talk * perf, data = debate)` |
| `rtmb_glm()` | 一般化線形モデル | `rtmb_glm(y ~ x, data = df, family = "bernoulli")` |
| `rtmb_lmer()` | 線形混合モデル | `rtmb_lmer(y ~ x + (1 | group), data = df)` |
| `rtmb_glmer()` | 一般化線形混合モデル | `rtmb_glmer(y ~ x + (1 | group), data = df, family = "poisson")` |
| `rtmb_corr()` | 相関・偏相関 | `rtmb_corr(cbind(sat, perf), data = debate)` |
| `rtmb_table()` | クロス表分析 | `rtmb_table(skill, cond, data = debate)` |
| `rtmb_loglinear()` | 対数線形モデル | `rtmb_loglinear(~ A + B + A:B, data = df)` |
| `rtmb_fa()` | 因子分析 | `rtmb_fa(items, nfactors = 3)` |
| `rtmb_irt()` | 項目反応理論 | `rtmb_irt(items, model = "2PL")` |
| `rtmb_lrt()` | 潜在ランク理論 | `rtmb_lrt(y, k = 3)` |
| `rtmb_mdu()` | 多次元展開法 | `rtmb_mdu(Y, ndim = 2)` |
| `rtmb_mediation()` | 媒介分析 | `rtmb_mediation(list(m ~ x, y ~ x + m), data = df)` |
| `rtmb_mixture()` | 混合分布モデル | `rtmb_mixture(y ~ x, k = 3, data = df)` |

`rtmb_lm()`、`rtmb_glm()`、`rtmb_lmer()`、`rtmb_glmer()`では、formula
の環境に変数がある場合は `data` を省略できます。省略時はformula内で変数名を
そのまま使い、`$`、`[[`、`.` を含む式では明示的に `data` を指定します。

IRTで `prior_normal()` を指定した場合、自由に推定される識別力パラメータには
既定で `a ~ exponential(1 / 2)`、すなわち平均2の指数事前分布が使われます。
rateを変更する場合は `prior_normal(a_rate = ...)` と指定します。

回帰系ラッパーでは、`std = TRUE` を指定すると、固定効果の標準化係数が `b_std` として追加されます。Gaussianモデルではデザイン行列の各列と応答変数の標準偏差を使い、それ以外のfamilyでは各列の標準偏差だけを使うため、係数はlink尺度に残ります。factorや交互作用から作られた列も対象です。

```{r}
mdl_std <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  std = TRUE
)
fit_std <- mdl_std$sample()
fit_std$summary()
```

# 推定のイメージ

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

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

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

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

```{r}
fit <- mdl$optimize()
fit$summary()
```

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

```{r}
fit <- mdl$classic()
fit$summary()
```

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

# 事後予測チェック

回帰系ラッパーのMAP、MCMC、変分推論の結果では、事後予測データを直接生成できます。
`type = "auto"` は、観測モデルが `_lpdf` なら密度、`_lpmf` なら棒グラフを使います。
MCMCと変分推論では事後標本を使います。サンプリング型の標準誤差を保存していない
MAPでは、MAP点に条件づけて複製データを生成します。

```{r}
yrep <- fit$posterior_predict(draws = 100, seed = 123)
fit$pp_check(type = "auto", draws = 100, seed = 123)
fit$pp_check(stat = mean, draws = 100, seed = 123)
```

`x` にモデルデータの列名、または観測数と同じ長さのベクトルを指定すると、
散布図形式のチェックになります。観測値と、事後予測平均および95%事後予測区間を
同じ図に表示します。予測にはモデル内のすべての項が使われ、`x` は横軸だけを選びます。

```{r}
fit$pp_check(x = "talk", draws = 100, seed = 123)
fit$pp_check(x = ".fitted", draws = 100, seed = 123)
```

予約値 `x = ".fitted"` を指定すると、各観測の事後予測平均を横軸に置きます。
モデル全体のキャリブレーションを確認するための表示になり、45度線も追加されます。
散布図の凡例は自動配置され、`legend_position`、`legend_cex`、または
`show_legend = FALSE` で調整できます。
予測区間はデフォルトで95%であり、`interval` で変更できます。

階層モデルでは `random` で予測水準を選びます。`"conditional"` は推定された
集団効果を使い、`"population"` は集団効果を除き、`"simulate"` は各事後標本から
新しい集団効果を生成します。

```{r}
fit_lmer$pp_check(random = "conditional")
fit_lmer$pp_check(random = "population")
fit_lmer$pp_check(random = "simulate")
```

# 頻度主義的分析

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

```{r}
rtmb_corr(cbind(sat, perf), data = debate)$classic()
```

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

```text
## 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 ***
```

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

```{r}
rtmb_corr(cbind(sat, perf),
          data = debate,
          covariates = ~ skill)$classic()
```

```text
##           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推定に基づく相関モデルとして推定されます。

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

## S3メソッド

`classic()` の結果は `Classic_Fit` オブジェクトなので、`anova()`、`AIC()`、`BIC()` などのS3メソッドを使えます。

クロス表分析では、1つの対象に対して `anova()` を使うと、独立性の検定を表示できます。

```{r}
fit_tab <- rtmb_table(skill, cond, data = debate)$classic()
anova(fit_tab)
```

```text
## 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 
```

複数のモデルを比較する場合は、`anova(fit1, fit2)` と書けます。たとえば、因子分析で因子数の異なるモデルを比較できます。

```{r eval=FALSE}
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)
```

```text
## 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
```

`AIC()` や `BIC()` も同じように使えます。

```{r eval=FALSE}
c(AIC_1 = AIC(fit_fa1),
  AIC_2 = AIC(fit_fa2),
  BIC_1 = BIC(fit_fa1),
  BIC_2 = BIC(fit_fa2))
```

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

# print_code()

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

```{r}
mdl <- rtmb_ttest(sat ~ cond, data = debate)
mdl$print_code()
```

```text
## === RTMB Model Code ===
## 
## rtmb_code(
##   setup = {
##     Y <- as.numeric(.data$response)
##     G <- as.integer(.data$group)
##     Y1 <- as.numeric(na.omit(Y[G == 1]))
##     Y2 <- as.numeric(na.omit(Y[G == 2]))
##   }, 
##   parameters = {
##     total_mean = Dim(1)
##     sd = Dim(1, lower = 0)
##     diff = Dim(1)
##   }, 
##   transform = {
##     delta <- diff/sd
##     mean0 <- total_mean + diff/2
##     mean1 <- total_mean - diff/2
##   }, 
##   model = {
##     Y1 ~ normal(mean0, sd)
##     Y2 ~ normal(mean1, sd)
##   }
## )
```

この既定例は等分散を仮定する2標本 t検定に対応し、2群で共通の標準偏差を推定します。`var.equal = FALSE` を指定すると、2群の標準偏差を別々に推定する Welch 型のモデルになります。

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

# 事後分析

回帰、GLM、GLMM のラッパーで交互作用を含むモデルを推定した場合は、`conditional_effects()` や `simple_effects()` を使って、予測値や単純効果を確認できます。

```{r}
mdl_int <- rtmb_lm(sat ~ talk * perf, data = debate)
fit_int <- mdl_int$sample()

ce <- conditional_effects(fit_int, effect = "talk:perf")
plot(ce)

simple_effects(fit_int, effect = "talk:perf")
```

MCMC fit では事後サンプルに基づいて区間を計算します。MAP fit や classic fit でも、標準誤差やシミュレーションに基づく不確実性が利用できる場合には、同じ関数で条件付き効果や単純効果を確認できます。

# 欠損値の扱い (Missing Data Handling)

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

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

| ラッパー関数 | 指定可能な `missing` オプション | デフォルト | 内容説明 |
| :--- | :--- | :--- | :--- |
| `rtmb_lm`, `rtmb_glm`, `rtmb_lmer`, `rtmb_glmer`, `rtmb_ttest` | `"listwise"` | `"listwise"` | モデル式（formula）に含まれる変数に `NA` が1つでもある行（観測値）をすべて除外します。 |
| `rtmb_fa`, `rtmb_corr` | `"listwise"`, `"fiml"` | `"listwise"` | `"listwise"` は完全データのみから共分散行列（十分統計量）を計算し高速に推定します。`"fiml"` は行ごとの欠損パターンに基づき、観測されている変数のみの周辺対数尤度を動的に足し合わせて推定します（完全情報最尤法）。 |
| `rtmb_corr` （追加オプション） | `"pairwise"` | | `method = "pearson"` または `"spearman"` のときのみ使用可能。R標準の `cor(..., use = "pairwise.complete.obs")` と完全に一致するペアワイズ削除結果を `classic()` メソッドなどで出力します。 |
| `rtmb_mdu` | `"listwise"`, `"fiml"` | `"listwise"` | `"fiml"` は評定データの欠損値（`NA`）を尤度計算ループ内で安全にスキップし、有効なデータのみで推定します。`"listwise"` は欠損を含む対象者を丸ごと除外します。 |
| `rtmb_irt` | `"listwise"`, `"fiml"` | `"fiml"` | `"fiml"` は項目ごとの欠損値を自動的にスキップし、手元にある回答のみを用いて能力や項目パラメータを推定します（デフォルト）。`"listwise"` を指定すると、1項目でも欠損がある回答者をすべて除外します。 |

### 分析例：因子分析における完全情報最尤法（FIML）

データに一部の欠損が含まれている場合、`missing = "fiml"` を指定することで、一部しか回答していない参加者も含めた全データを用いて因子負荷量を頑健に推定できます。

```r
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()
```

### 分析例：相関分析におけるペアワイズ削除（Pairwise）

R標準の `cor()` や `cor.test()` のペアワイズ削除の挙動と一致させたい場合は、`missing = "pairwise"` を指定して `classic()` を実行します。

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

---

# 次に読む記事

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

- [イントロダクション](ja-introduction.html): BayesRTMB の基本的な考え方と、MCMC、MAP推定、変分推論、頻度主義的分析の位置づけ。
- [クイックスタート](ja-quick_start.html): インストール、最初のモデル、MCMC診断、可視化、t検定のベイズファクター。
- [階層モデル・GLMM・分散分析](ja-rtmb_glmer.html): `rtmb_glmer()` を使って mixed model / GLMM を詳しく扱う方法。
- [モデルを書く](ja-writing_models.html): `rtmb_code()` を使って独自のモデルを書く方法。
- [内部構造](ja-rtmb_internals.html): `setup`、`parameters`、`transform`、`model` など、BayesRTMB の内部的なモデル構造。
- [分析リファレンス](ja-analysis_reference.html): fit object のメソッド、モデル比較、固定パラメータ、分布、AD テープ化の参照情報。
