GAMMだけじゃない！Rにおける一般化加法混合モデルの代替手法を比較

2025-05-16

一般化加法混合モデル (GAMMs) とは

GAMMs（ギャムズ）は、統計モデリング手法の一つで、以下の2つの主要な概念を組み合わせたものです。

一般化加法モデル (GAMs: Generalized Additive Models): 線形回帰モデルや一般化線形モデル（GLMs）を拡張したもので、説明変数と目的変数との関係を線形ではなく、非線形な関数（平滑化スプラインなど）の和で表現することを可能にします。これにより、データに内在する複雑な非線形なパターンを柔軟に捉えることができます。例えば、時間経過に伴う反応の変化や、ある環境要因が特定の範囲で非線形に影響する場合などに有効です。
混合モデル (Mixed Models): データに階層構造がある場合（例：同じ被験者から複数回測定されたデータ、複数の学校から生徒のデータが収集された場合など）に用いられます。個体差やグループ間のばらつきを「ランダム効果」としてモデルに組み込むことで、データの独立性の仮定が満たされない場合でも適切な分析が可能になります。固定効果（固定されたパラメータ）とランダム効果（ランダムな変動を持つパラメータ）の両方を考慮する点が特徴です。

GAMMsは、これらの概念を組み合わせることで、階層的な構造を持つデータにおける非線形な関係性を分析するための強力なツールとなります。特に、時系列データや、被験者内デザインの実験データ、地理空間データなど、複雑な構造を持つデータに対して威力を発揮します。

GAMMsのRにおける実装

RでGAMMsを扱うには、主にmgcvパッケージが使われます。このパッケージは、GAMsおよびGAMMsを構築・推定するための非常に強力で柔軟な機能を提供します。

基本的なGAMMの式

GAMMの一般的な式は、以下のように表すことができます。

g(E[Y])=beta_0+f_1(X_1)+f_2(X_2)+dots+Zgamma+epsilon

ここで、

epsilon：誤差項
Zgamma：ランダム効果の部分。Zはランダム効果に対応する計画行列、gammaはランダム効果の係数ベクトル（通常、平均ゼロの正規分布に従うと仮定されます）。
f_i(X_i)：説明変数X_iに対する平滑化関数（非線形な関係を捉える部分）。通常、スプラインなどの基底関数で表現されます。
beta_0：切片
E[Y]：Yの期待値
g()：リンク関数（GLMと同様に、目的変数の種類に応じて適切に選択されます。例：正規分布の場合は恒等関数、二項分布の場合はロジット関数、ポアソン分布の場合は対数関数など）
Y：目的変数

`mgcv`パッケージを使ったGAMMの例

mgcvパッケージでは、gam()関数がGAMMのモデリングに用いられます。ランダム効果を組み込むには、s()関数（平滑化項）の中でby引数を使うか、またはre()関数を使います。

例：被験者ごとの非線形なトレンドと全体的な非線形トレンドを考慮するケース

ある反応（response）が、時間（time）と被験者（subject）によって影響を受けるデータを考えます。このとき、時間の影響は非線形であり、かつ被験者ごとにその非線形なトレンドが異なる可能性があるとします。

# mgcvパッケージをロード
library(mgcv)

# ダミーデータの作成
set.seed(123)
n_subjects <- 20
n_obs_per_subject <- 50
data <- data.frame(
  subject = factor(rep(1:n_subjects, each = n_obs_per_subject)),
  time = rep(seq(0, 10, length.out = n_obs_per_subject), n_subjects)
)
# 真の非線形効果とランダム効果を加える
data$response <- 2 * sin(data$time / 2) + # 全体的な非線形トレンド
                 rnorm(n_subjects, mean = 0, sd = 0.5)[data$subject] * data$time + # 被験者ごとの線形なランダムスロープ
                 rnorm(n_subjects * n_obs_per_subject, sd = 0.8) # 誤差

GAMMのモデル構築

このデータに対して、以下のようなGAMMを構築できます。

# モデルの式:
# s(time): 全体的な非線形効果
# s(time, by = subject, bs = "re"): 被験者ごとのランダムな平滑化効果（ランダムスプライン）
# または s(subject, bs = "re"): ランダム切片 (被験者ごとの平均レベルの差)
# + s(time, subject, bs="fs", m=1): 被験者ごとの非線形なトレンド (ファクタースプライン、より柔軟)

# パターン1：全体的な非線形効果 + 被験者ごとのランダム切片 + 被験者ごとのランダムスロープ
# mgcvでは、ランダム効果は 's(factor_variable, bs="re")' の形式で指定できる
model1 <- gam(response ~ s(time) + s(subject, bs = "re") + s(time, by = subject, bs = "re"),
              data = data,
              method = "REML") # REMLは混合モデルで推奨される推定方法

summary(model1)
plot(model1, pages = 1) # 各平滑化項のプロット

# パターン2：全体的な非線形効果 + 被験者ごとの非線形なランダム効果（より柔軟なモデル）
# この形式では、s(time)が全体的な非線形トレンド、s(time, subject)が各被験者ごとの時間に対する非線形な変動を表す。
# bs="fs" (factor smooth) は、factor_variableによってスムーズが異なることを意味する。
model2 <- gam(response ~ s(time) + s(time, by = subject, bs = "fs"),
              data = data,
              method = "REML")

summary(model2)
plot(model2, pages = 1, select = 1) # 全体的な非線形トレンド
plot(model2, pages = 1, select = 2) # 被験者ごとの非線形な変動 (一部の被験者の例)

# 注意：s(time, by = subject, bs = "re") と s(time, subject, bs = "fs") は
# どちらもランダムな非線形効果を表現できますが、"fs"はより一般的な形で、
# 各レベル（subject）に対して個別のスムーズ関数をフィットします。
# "re"は純粋なランダム効果の平滑化版であり、ランダム切片やランダムスロープの非線形版と解釈できます。

メリット

解釈可能性: 各平滑化項を個別にプロットすることで、各説明変数が目的変数にどのように影響しているかを視覚的に理解しやすいです。
様々なデータタイプに対応: GLMsのフレームワークに基づいているため、正規分布だけでなく、二項分布（ロジスティック回帰）、ポアソン分布（ポアソン回帰）など、様々な目的変数に対応できます。
階層的データの考慮: ランダム効果を導入することで、データ内のグループ構造や繰り返し測定の特性を適切にモデル化し、より正確な推定と推論が可能になります。
非線形性の柔軟なモデリング: 説明変数と目的変数間の複雑な非線形関係を、事前に形式を仮定することなくデータから学習できます。

デメリット

適切な平滑化パラメータの選択: mgcvパッケージは自動的に平滑化パラメータを推定しますが、最適な選択には注意が必要です。過学習や過小学習のリスクがあります。
計算コスト: 大規模なデータセットや多くの平滑化項を含む場合、計算に時間がかかることがあります。
モデルの複雑さ: 線形モデルに比べてモデルが複雑になるため、解釈が難しくなる場合があります。特に、複数の平滑化項や相互作用がある場合。

GAMMsにおけるよくあるエラーとトラブルシューティング

GAMMsのモデリングは、非線形性と階層構造という2つの複雑な要素を同時に扱うため、通常の線形モデルやGLMに比べてエラーや警告が出やすくなります。

収束の問題 (Convergence Issues)

最も頻繁に遭遇する問題の一つです。モデルの推定アルゴリズムが最適解を見つけられない場合に発生します。

考えられる原因と対処法:
1. 過剰なモデルの複雑さ (Overly Complex Model):
  - 原因: データ量に対して平滑化項（s()）の数やk（基底次元）が多すぎる、またはランダム効果の構造が複雑すぎる場合に発生します。特に、カテゴリカル変数のレベル数が少ない場合にby引数を用いたランダム平滑化（s(time, by=subject)）を設定すると、各グループで十分なデータがないために不安定になります。
  - 対処法:
    - kの値を減らす（例: s(x, k=5)）。
    - モデルを単純化する。不要な平滑化項やランダム効果を削除する。
    - ランダム効果の構造を見直す。例えば、ランダムスロープが収束しない場合、まずはランダム切片（s(subject, bs="re")）のみで試す。
    - mgcv::gam()のmethodを"REML"（制限付き最尤推定）から"ML"（最尤推定）に変更する。"REML"の方がバイアスが小さい傾向がありますが、"ML"の方が収束しやすい場合があります（特に小サンプルで）。
2. 共線性/コンカービティ (Collinearity/Concurvity):
  - 原因: 複数の説明変数が互いに強く相関している場合（共線性）や、複数の平滑化項が同じような効果を捉えようとしている場合（コンカービティ）に発生します。GAMsでは、平滑化項同士でもこの問題が起こり得ます。
  - 対処法:
    - gam.check()関数を使ってコンカービティを確認する。特に、concurvity()関数で数値的な指標を確認できます。
    - 相関の高い変数をモデルから除外する、または組み合わせる。
    - 複数の平滑化項が同じようなトレンドを捉えている場合は、それらを統合するか、一方を削除する。
    - by引数を使ったランダム平滑化（例: s(time, by=subject)）と、全体的な平滑化（s(time)）が同時に存在する場合、重複する情報をモデルが推定しようとして問題となることがあります。このような場合は、s(time, subject, bs="fs")のようなファクタースムーズ（factor smooth）を使うことで、全体効果と各グループの偏差を同時にモデル化できます。
3. データの問題 (Data Issues):
  - 原因:
    - データが少なすぎる: 特にランダム効果や平滑化項が多い場合、十分なデータがないと安定した推定ができません。各グループのデータ数や、非線形性を推定するために必要なデータポイントが不足している可能性があります。
    - 変数のスケールが極端に異なる: 数値的な安定性に影響を与えることがあります。
    - 外れ値: 極端な外れ値が存在すると、モデルの推定が困難になることがあります。
    - 目的変数の分布とリンク関数の不一致: 例えば、カウントデータなのに正規分布を仮定しているなど。
    - オーバーディスパーション (Overdispersion): ポアソン分布や二項分布のデータで、分散が平均から期待される以上に大きい場合に発生します。
  - 対処法:
    - データの量を確認し、必要であればデータを増やす。
    - 連続変数をセンタリングまたはスケールする（例: scale()関数）。
    - 外れ値を特定し、その影響を調べる（必要に応じて除去または変換）。
    - 適切な目的変数分布とリンク関数を選択する。カウントデータにはポアソン分布の代わりに負の二項分布（family=nb()）を使う、二項データでオーバーディスパーションがある場合は擬似二項分布（family=quasibinomial()）を使う。
    - オーバーディスパーションの場合、gam.check()の残差プロットを確認するか、dispersionパラメータを推定する。
4. 最適化アルゴリズムの問題 (Optimizer Issues):
  - 原因: gam()関数が使用する最適化アルゴリズムが、デフォルトの設定では最適解に到達できない場合があります。
  - 対処法:
    - control引数を使って最適化設定を変更する。特にnthreads（スレッド数）やmaxit（最大反復回数）を増やす。
    - mgcvの内部で使われている最適化ルーチンを変更する（通常は自動で良いが、特殊なケースでは調整が必要な場合もある）。
エラー/警告メッセージの例:
- Error: system is computationally singular: reciprocal condition number ...
- Warning: Model failed to converge
- Warning: Model is nearly unidentifiable: very large eigenvalue
- Error in solve.default(t(X) %*% X) : system is computationally singular: reciprocal condition number ...

モデルの適合不良 (Poor Model Fit)

モデルが収束しても、データへの適合が悪い場合があります。

問題の症状:
- gam.check()の出力でk-indexが1に近い、または有意なP値を示す。これは、指定した基底次元kが小さすぎて、モデルがデータの非線形性を十分に捉えられていない（under-fitting）ことを示します。
- 残差プロットにパターンが見られる（残差が正規分布に従わない、分散が一定でないなど）。
- 予測が不自然に滑らかすぎたり、逆にギザギザしすぎたりする。

モデルの解釈の問題 (Interpretation Issues)

モデルはフィットしたが、結果の解釈が難しい場合があります。

考えられる原因と対処法:
1. 平滑化項のセンタリング (Centering of Smooths):
  - 原因: mgcvの平滑化項は、通常、平均がゼロになるように制約されます。これにより、平滑化項の切片や線形部分が固定効果の推定値に吸収されることがあります。
  - 対処法: これ自体はエラーではありませんが、解釈に注意が必要です。plot()関数を使って各平滑化項を個別に可視化することで、その変数の非線形な影響を理解できます。もし、全体的な切片や線形効果を明確にしたい場合は、平滑化項に加えて線形項（例: y ~ x + s(x)）もモデルに含めることが検討されますが、これは共線性につながる可能性があるため注意が必要です。
2. 相互作用項の複雑さ (Complexity of Interaction Terms):
  - 原因: GAMMsでは、平滑化項間の相互作用（例: s(x1, x2)）や、線形項と平滑化項の相互作用（例: x1:s(x2)）をモデル化できますが、これらは解釈が複雑になります。
  - 対処法:
    - 2D平滑化の場合は、plot()関数のperspまたはcontour引数で3Dプロットや等高線プロットを作成し、視覚的に把握する。
    - by引数を使った相互作用（例: s(x, by=factor)）の場合、各ファクターレベルごとに平滑化項のプロットを確認する。
    - marginaleffectsパッケージなどのツールを使って、特定の条件下での予測値の変化をプロットすることで、相互作用の効果をより具体的に理解できます。
問題の症状:
- summary()の出力が理解しにくい。
- plot()で描かれる平滑化項のプロットが予想と異なる。
- 特定の変数の効果がどのように現れているか直感的に理解できない。

専門家への相談:
- 上記の手順を試しても解決しない場合は、統計学の専門家や、mgcvパッケージに詳しいコミュニティ（Stack Overflowなど）に相談することを検討しましょう。
Rとmgcvパッケージのアップデート:
- 古いバージョンのRやパッケージでは、バグが存在する場合があります。常に最新の安定版を使用することを推奨します。
エラーメッセージの検索:
- Rのエラーメッセージは非常に具体的です。エラーメッセージをコピーして、Google検索やStack Overflowなどのフォーラムで検索すると、同じ問題に遭遇した他のユーザーの解決策が見つかることが多いです。
gam.check()の活用:
- mgcv::gam()でモデルをフィットした後、必ずgam.check(model_object)を実行し、残差の診断プロットとk-indexを確認する。これは、モデルの適合度と基底次元の適切性を評価する上で非常に重要です。
データとモデルの再確認:
- すべての変数が適切なデータ型（factor, numericなど）であることを確認する。特にカテゴリカル変数はfactor型にする。
- 欠損値がないか確認する。欠損値があると、デフォルトでは該当する行が削除されます。
- 変数の名前や式のスペルミスがないか確認する。
- 最初にシンプルなモデルから始め、徐々に複雑化していく。線形モデル (lm()) やGLM (glm()) で試してみて、問題なく動くか確認する。

GAMMsは、非線形な関係を柔軟にモデル化できる「一般化加法モデル (GAMs)」と、階層的なデータ構造や反復測定データに対応できる「混合モデル (Mixed Models)」を組み合わせたものです。これにより、個体差やグループ間のばらつきを考慮しつつ、説明変数と応答変数間の非線形な関係を捉えることができます。

以下の例では、架空のデータセットを作成し、様々なGAMMsの構築方法と、結果の解釈、診断方法について示します。

事前準備: 必要なパッケージのインストールと読み込み

# 必要なパッケージをインストール（まだインストールしていない場合）
# install.packages("mgcv")
# install.packages("ggplot2") # プロット用
# install.packages("dplyr")   # データ操作用

# パッケージの読み込み
library(mgcv)
library(ggplot2)
library(dplyr)

例1: 基本的なGAMM (ランダム切片)

最もシンプルなGAMMの例として、各被験者（subject）が異なるベースラインを持つ（ランダム切片）と同時に、時間（time）に対する非線形な効果がある場合を考えます。

データの作成

set.seed(123) # 再現性のためのシード設定
n_subjects <- 20 # 被験者数
n_obs_per_subject <- 50 # 各被験者からの観測数

# ダミーデータの生成
data_basic <- expand.grid(time = seq(1, 50, length.out = n_obs_per_subject),
                          subject = factor(1:n_subjects)) %>%
  as.data.frame()

# 被験者ごとのランダム切片
subject_intercepts <- rnorm(n_subjects, mean = 0, sd = 2)
data_basic$random_intercept <- subject_intercepts[as.numeric(data_basic$subject)]

# 真の非線形効果（時間）
data_basic$true_nonlinear_effect <- 3 * sin(data_basic$time / 10) + 0.5 * data_basic$time

# 応答変数 (response) の生成
data_basic$response <- 10 + data_basic$true_nonlinear_effect +
                       data_basic$random_intercept +
                       rnorm(nrow(data_basic), mean = 0, sd = 1.5)

# データの確認
head(data_basic)

GAMMの構築と診断

mgcvパッケージのgam()関数を使用します。

s(subject, bs="re"): subjectに対するランダム切片。bs="re"は「random effect basis」を意味し、各subjectのカテゴリカル変数に対して独立した正規分布のランダム効果を仮定します。
s(time): timeに対する非線形な平滑化項（通常、thin plate splineがデフォルト）。

# GAMMのモデル構築
# method="REML" は、混合モデルで推奨される平滑化パラメータ推定方法
model_basic <- gam(response ~ s(time) + s(subject, bs = "re"),
                   data = data_basic,
                   method = "REML")

# モデルの概要を表示
summary(model_basic)

# モデルの診断プロット
# 残差の診断、k（基底次元）のチェックなど
gam.check(model_basic)

# 平滑化項のプロット
# 全体的な時間の効果と、各被験者のランダム切片のプロット
plot(model_basic, pages = 1, residuals = TRUE, pch = 1) # pages=1で全てを1ページに、residuals=TRUEで残差もプロット

結果の解釈

summary(model_basic)の出力:

edf (estimated degrees of freedom): 平滑化項の「柔軟性」の度合い。s(time)はkに近い値（デフォルトではk=10）になり、非線形性を捉えていることを示します。s(subject)はedfがn_subjects - 1に近い値になり、各被験者ごとのランダムなばらつきを捉えていることを示します。
Approximate significance of smooth terms: 平滑化項の有意性。s(time)は非線形効果が統計的に有意であることを示し、s(subject)はランダム切片の分散が有意であることを示します。
Parametric coefficients: 固定効果（この場合は切片のみ）。

plot(model_basic)の出力:

s(subject)のプロット: 各subjectにおける切片からの偏差（ランダム切片の推定値）が表示されます。
s(time)のプロット: timeがresponseに与える非線形な影響の推定値と、その95%信頼区間が表示されます。

例2: ランダムスロープと非線形なランダム効果

時間経過とともに被験者ごとの傾きが異なる場合（ランダムスロープ）、または被験者ごとに非線形なトレンドが異なる場合を考えます。

set.seed(456)
data_complex <- expand.grid(time = seq(1, 50, length.out = n_obs_per_subject),
                            subject = factor(1:n_subjects)) %>%
  as.data.frame()

# 被験者ごとのランダム切片
subject_intercepts_cplx <- rnorm(n_subjects, mean = 0, sd = 2)
data_complex$random_intercept <- subject_intercepts_cplx[as.numeric(data_complex$subject)]

# 被験者ごとのランダムスロープ（時間の効果が被験者ごとに異なる）
# ここでは、時間に対する線形なランダムスロープを例として含める
subject_slopes_cplx <- rnorm(n_subjects, mean = 0, sd = 0.1)
data_complex$random_slope_effect <- subject_slopes_cplx[as.numeric(data_complex$subject)] * data_complex$time

# 被験者ごとの非線形なランダム効果の基盤（より複雑なパターン）
# これは s(time, by=subject, bs="re") や s(time, subject, bs="fs") で捉える
subject_nonlinear_dev <- sapply(1:n_subjects, function(i) {
  0.5 * sin(data_complex$time[data_complex$subject == i] / sample(5:15, 1))
})
data_complex$random_nonlinear_effect <- as.vector(subject_nonlinear_dev)

# 全体的な非線形効果
data_complex$true_nonlinear_effect <- 5 * cos(data_complex$time / 8)

# 応答変数の生成
data_complex$response <- 20 + data_complex$true_nonlinear_effect +
                         data_complex$random_intercept +
                         data_complex$random_slope_effect +
                         data_complex$random_nonlinear_effect +
                         rnorm(nrow(data_complex), mean = 0, sd = 2)

head(data_complex)

GAMMの構築（ランダムスロープと非線形なランダム効果）

ランダム切片とランダムスロープ（線形） s(time, by = subject, bs = "re") は、timeの効果がsubjectによって異なるランダムスロープを意味しますが、bs="re"を使うと各subjectで線形に変動するランダム効果（より正確には、時間の効果のランダムな「傾き」の平滑化版）を表現します。

model_rs <- gam(response ~ s(time) +                # 全体的な非線形効果
                         s(subject, bs = "re") +    # ランダム切片
                         s(time, by = subject, bs = "re"), # 被験者ごとの時間に対するランダム効果（スロープ）
                   data = data_complex,
                   method = "REML")

summary(model_rs)
gam.check(model_rs)
plot(model_rs, pages = 1)

ファクタースムーズ (Factor Smooth): 被験者ごとの非線形なトレンド s(time, subject, bs="fs") は、subjectの各レベルに対してtimeの非線形効果を推定します。これは、各被験者で異なる非線形なパターンを捉えたい場合に非常に柔軟なアプローチです。m=1は、各グループのスムーズが全体平均から大きく逸脱することにペナルティを与えることで、過学習を防ぐのに役立ちます。

model_fs <- gam(response ~ s(time) +                # 全体的な非線形効果
                         s(time, subject, bs = "fs", m = 1), # 被験者ごとの非線形トレンド（factor smooth）
                   data = data_complex,
                   method = "REML")

summary(model_fs)
gam.check(model_fs)
plot(model_fs, pages = 1) # 各被験者の非線形効果のプロットが含まれる

plot(model_fs)の出力では、s(time)が全体の非線形効果を示し、s(time, subject)が各被験者における全体効果からの非線形な偏差を示します。これにより、各被験者固有の時間トレンドを視覚的に確認できます。

応答変数が二値の場合（例: 成功/失敗、Yes/No）には、二項分布とロジットリンク関数を使用します。

set.seed(789)
data_binary <- expand.grid(time = seq(1, 50, length.out = n_obs_per_subject),
                           subject = factor(1:n_subjects)) %>%
  as.data.frame()

# 被験者ごとのランダム切片（対数オッズスケール）
subject_intercepts_bin <- rnorm(n_subjects, mean = -1, sd = 1)
data_binary$random_intercept_logit <- subject_intercepts_bin[as.numeric(data_binary$subject)]

# 時間の非線形効果（対数オッズスケール）
data_binary$true_nonlinear_effect_logit <- 0.1 * data_binary$time - 0.002 * data_binary$time^2

# 線形予測子 (logit scale)
data_binary$linear_predictor <- -2 + data_binary$true_nonlinear_effect_logit +
                               data_binary$random_intercept_logit

# 確率 (probability) に変換
data_binary$probability <- plogis(data_binary$linear_predictor) # plogisはロジット関数の逆関数

# 応答変数 (binary_response) の生成
data_binary$binary_response <- rbinom(nrow(data_binary), size = 1, prob = data_binary$probability)

head(data_binary)

family = binomial() を指定します。

# 二項分布のGAMM
model_binary <- gam(binary_response ~ s(time) + s(subject, bs = "re"),
                    data = data_binary,
                    family = binomial(),
                    method = "REML")

summary(model_binary)
gam.check(model_binary)
plot(model_binary, pages = 1)

plot(model_binary)のプロットは、線形予測子スケール（対数オッズスケール）での平滑化効果を示します。

これらの例は、Rのmgcvパッケージを用いたGAMMsの基本的な使い方を示しています。

summary()でモデルの要約を、gam.check()でモデルの診断を、plot()で平滑化項を視覚的に確認します。
method="REML"は混合モデルで推奨される推定方法です。
family引数で応答変数の分布を指定します。
by引数やbs="fs"を使って、より複雑なランダム効果（ランダムスロープや非線形なランダム効果）を表現できます。
bs="re"でランダム効果をモデル化します。
s()関数で平滑化項を指定し、kで基底次元を調整できます。
gam()関数が主要な関数です。

lme4 パッケージとの組み合わせ (gamm 関数)

mgcv パッケージには、lme4 パッケージの機能を利用してGAMMsをフィットするための gamm 関数が含まれています。この関数は、特に誤差項の相関構造や、ランダム効果のより複雑な構造を明示的に指定したい場合に有用です。

特徴

計算の安定性: mgcv::gam よりも大規模なデータセットや複雑なランダム効果モデルで収束が安定する場合もあります。
nlme パッケージの機能を利用: gamm 関数は内部で nlme パッケージ（以前は lme4 だったが、現在は nlme を推奨）の lme 関数を呼び出します。これにより、自己相関（corARMA）やヘテロスケダスティシティ（varIdent など）といった、より多様な残差構造をモデルに組み込むことができます。

コード例

# install.packages("nlme") # 必要な場合
library(mgcv)
library(nlme) # gamm関数が内部で利用

# 例1のデータを使用 (data_basic)
# response ~ s(time) + s(subject, bs="re") を gamm で表現

# gamm() の構文:
# fixed = 固定効果の式 (GAMsの平滑化項もここに含める)
# random = ランダム効果の式 (nlme::lme の形式)
# data = データフレーム
# correlation = 誤差項の相関構造 (オプション)
# weights = 残差の分散構造 (オプション)
# family = GLMのファミリー (例: binomial, poisson)

# ランダム切片モデル
model_gamm_basic <- gamm(response ~ s(time),
                         random = list(subject = ~ 1), # subject ごとにランダム切片
                         data = data_basic,
                         method = "REML")

# 結果の確認 (lme と gam の両方の結果が得られる)
summary(model_gamm_basic$gam) # GAMs部分のサマリー
summary(model_gamm_basic$lme) # 混合モデル部分のサマリー

plot(model_gamm_basic$gam) # 平滑化項のプロット

# 時系列データにおける自己相関の考慮 (例: subject 内での時間による自己相関)
# この場合は、各 subject 内で time が順序付けられている必要があります
# data_basic を subject と time でソートしておく
data_basic_sorted <- data_basic %>% arrange(subject, time)

# AR(1) 相関を考慮したモデル
# form = ~ time | subject は、subject ごとに time の順序で相関を考慮することを示す
model_gamm_corr <- gamm(response ~ s(time),
                        random = list(subject = ~ 1),
                        correlation = corARMA(form = ~ time | subject, p = 1, q = 0), # AR(1)
                        data = data_basic_sorted,
                        method = "REML")

summary(model_gamm_corr$gam)
summary(model_gamm_corr$lme)

glmmTMB パッケージ

glmmTMB は、一般化線形混合モデル（GLMMs）を構築するための非常に柔軟なパッケージですが、非線形な共変量の効果を平滑化項として組み込むことも可能です。mgcv のような自動的な平滑化選択機能は少ないですが、多様な分布（ゼロ過剰データなど）やランダム効果構造に対応できる点が強みです。

高速な推定: Template Model Builder (TMB) フレームワークを使用しており、大規模なデータセットや複雑なモデルでの推定が高速です。
柔軟なランダム効果: 多重メンバーシップ（multiple membership）や空間的なランダム効果など、複雑なランダム効果構造を定義できます。
柔軟な分布: ゼロ過剰（zero-inflated）モデル、コンパウンドポアソン、ベータ分布など、mgcv の family オプションでは直接利用できない多くの分布をサポートします。

mgcv の s() 関数に相当する機能を glmmTMB で実現するには、splines パッケージの基底関数（例: ns for natural splines, bs for B-splines）を固定効果の式に含めるか、mgcv からsmooth.constructを利用する必要があります。しかし、より直接的な方法としては、glmmTMB 内でGAMsの平滑化を模倣する「GAM-like」なランダム効果を使用することが可能です。

# install.packages("glmmTMB") # 必要な場合
library(glmmTMB)
library(splines) # スプライン基底用

# 例1のデータを使用 (data_basic)

# glmmTMB でのランダム切片と平滑化効果のモデリング
# smooth.spline() のように、時間に対する平滑化効果をモデル化
# s() とは異なり、平滑化の柔軟性は手動で調整する必要がある
# ここでは、ns (natural spline) を使用して非線形効果を近似
# ランダム効果は (1 | subject) で指定

model_glmmTMB_basic <- glmmTMB(response ~ ns(time, df = 5) + (1 | subject), # df は自由度（柔軟性）
                               data = data_basic,
                               REML = TRUE) # REML推定

summary(model_glmmTMB_basic)

# glmmTMB でのランダムスロープも同様に指定可能
# model_glmmTMB_rs <- glmmTMB(response ~ ns(time, df = 5) + (time | subject), # subject ごとに time のランダムスロープ
#                             data = data_basic,
#                             REML = TRUE)
# summary(model_glmmTMB_rs)

# glmmTMB でのゼロ過剰ポアソンGAMMの例 (ダミーデータ)
# data_zi <- data.frame(
#   x = runif(100),
#   group = factor(rep(1:10, each = 10)),
#   # 真の発生プロセスを定義
#   lambda = exp(1 + 2 * sin(runif(100) * pi * 2) + rnorm(100, 0, 0.5)),
#   # ゼロ過剰の確率
#   p_zero = plogis(-1 + rnorm(100, 0, 0.5)[rep(1:10, each=10)])
# )
# data_zi$y <- rpois(100, data_zi$lambda)
# data_zi$y[runif(100) < data_zi$p_zero] <- 0 # ゼロを挿入

# model_glmmTMB_zi <- glmmTMB(y ~ ns(x, df = 4) + (1 | group),
#                             ziformula = ~ 1, # ゼロ過剰部分のモデル
#                             family = poisson(),
#                             data = data_zi)
# summary(model_glmmTMB_zi)

ベイズ推定 (brms パッケージ)

ベイズ統計モデリングフレームワークである brms は、Stan をバックエンドに用いて、GAMMsを含む非常に広範な統計モデルを構築できます。特に事前分布の組み込みや、推定の不確実性の完全な評価をしたい場合に強力です。

事前分布の指定: ドメイン知識に基づいて事前分布を指定することで、推定の安定性を高めたり、結果に影響を与えることができます。
多様な分布: GLMsでサポートされるものから、ロバストな分布（t分布など）、ゼロ過剰分布、順序ロジスティック回帰など、非常に多様な応答変数に対応します。
柔軟なモデル定義: mgcv の s() 関数に似た構文で平滑化項を定義できます。ランダム効果も柔軟に指定できます。
ベイズアプローチ: パラメータの点推定だけでなく、完全な事後分布が得られます。これにより、より詳細な不確実性の評価や、複数のパラメータ間の関係の理解が可能になります。

brms は mgcv の構文にかなり似ていますが、計算に時間がかかります。

# install.packages("brms") # 必要な場合
# install.packages("rstan") # brms のバックエンド
library(brms)
library(dplyr) # データ操作用

# 例1のデータを使用 (data_basic)

# brms での GAMM 構築
# s() 関数は mgcv と同様に平滑化を指定
# (1 | subject) はランダム切片を指定
# family = gaussian() は正規分布
# cores = 4 などで並列処理 (環境に合わせて調整)
# iter, warmup は MCMC のイテレーション数とウォームアップ数 (収束のために増やす必要あり)

# 計算に時間がかかるため、小規模なデータで試すか、iter を減らすことを検討
# model_brms_basic <- brms(response ~ s(time) + (1 | subject),
#                          data = data_basic,
#                          family = gaussian(),
#                          chains = 4, # MCMCチェーンの数
#                          iter = 2000, # イテレーション数
#                          warmup = 1000, # ウォームアップ（バーンイン）数
#                          cores = 4, # 使用するCPUコア数
#                          seed = 123)

# summary(model_brms_basic)
# plot(model_brms_basic) # パラメータの事後分布
# conditional_effects(model_brms_basic, effects = "time") # 平滑化効果のプロット

# 平滑化項のプロット (brms の plot 関数は s() のプロットも含む)
# plot(model_brms_basic, plot = "smooths") # これだと全スムーズが表示される

# 特定のスムーズをプロット
# plot(conditional_effects(model_brms_basic, effects = "time"), points = TRUE)

brms: **ベイズ推論の利点（事前分布の組み込み、事後分布の完全な評価、不確実性の明示的な定量化）**を享受したい場合に適しています。特に、複雑なモデルで推定が不安定な場合や、事前情報を利用したい場合に検討する価値があります。ただし、計算コストは最も高くなります。
glmmTMB: 多様な分布（特にゼロ過剰データ）や、mgcv では直接サポートされていない複雑なランダム効果構造が必要な場合に非常に強力な選択肢です。ただし、平滑化の柔軟性は手動で調整するか、スプライン基底関数を明示的に指定する必要があります。
mgcv::gamm: mgcv に慣れていて、誤差項に特定の相関構造や分散構造を導入したい場合に最適です。特に時系列データや空間データで残差の自己相関をモデル化したい場合に役立ちます。