Rで学ぶ統計モデリング：GAM(統合的平滑性推定)のコード例と実践

統合的平滑性推定を伴う一般化加法モデル (Generalized Additive Models with Integrated Smoothness Estimation) in R

「Generalized Additive Models with Integrated Smoothness Estimation」、略して GAMs with integrated smoothness estimation は、「R」で利用できる強力な統計モデリング手法の一つです。これは、伝統的な線形モデルや一般化線形モデル（GLMs）を拡張したもので、応答変数と予測変数の間の非線形な関係を柔軟に捉えることができます。特に、平滑化関数 (smooth functions) を用いることで、データに内在する複雑なパターンを学習します。

GAMs の基本的な考え方

従来の線形モデルでは、予測変数と応答変数の関係は以下のように線形に表現されます。

g(E(Y))=β0​+β1​X1​+β2​X2​+⋯+βp​Xp​

ここで、g(cdot) はリンク関数、E(Y) は応答変数 Y の期待値、X_i は予測変数、beta_i は回帰係数です。

一方、GAMs では、線形項の代わりに、予測変数の非線形な関数（平滑化関数）を用います。

g(E(Y))=β0​+f1​(X1​)+f2​(X2​)+⋯+fp​(Xp​)

ここで、f_i(X_i) が予測変数 X_i の滑らかな非線形関数を表します。これらの平滑化関数は、スプライン関数（例：thin plate regression splines, cubic regression splines）や局所重み付け散布図平滑化（LOESS）など、さまざまな種類があります。

統合的平滑性推定 (Integrated Smoothness Estimation)

GAMs の重要な側面の一つが、平滑性の推定 です。平滑化関数は、データに過剰に適合（過学習）しないように、ある程度の「滑らかさ」を持つ必要があります。この滑らかさを制御するために、通常は 正則化 (regularization) の手法が用いられます。

統合的平滑性推定では、この平滑性の度合い（つまり、どの程度非線形な関係を許容するか）をモデルの推定プロセスの中で自動的に決定します。これは、罰則付き最小二乗法 (penalized least squares) や 罰則付き尤度最大化 (penalized likelihood maximization) といった最適化手法を通じて行われます。

具体的には、平滑化関数にペナルティ項（罰則項）を導入し、モデルの当てはまりの良さと関数の複雑さ（滑らかさの度合い）のバランスを取ります。このペナルティ項は、関数の曲率や変動の大きさを測るもので、これが大きいほど関数はより「ギザギザ」になります。モデルの推定では、このペナルティ項を最小化するように平滑化関数の形状が調整されます。

R での GAMs の利用

R では、主に mgcv パッケージの gam() 関数を用いて GAMs を実装します。gam() 関数では、さまざまな種類の平滑化関数を指定したり、統合的平滑性推定を自動的に行うことができます。

基本的な gam() 関数の使用例は以下の通りです。

library(mgcv)

# サンプルデータの作成
set.seed(123)
n <- 100
x <- sort(runif(n))
true_y <- sin(2 * pi * x) + 0.5 * x
y <- true_y + rnorm(n, sd = 0.2)
data <- data.frame(x = x, y = y)

# GAM のモデル fitting
model <- gam(y ~ s(x), data = data)

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

# 予測
new_x <- seq(0, 1, length.out = 100)
predictions <- predict(model, newdata = data.frame(x = new_x))

# 結果のプロット
plot(data$x, data$y, main = "GAM with Integrated Smoothness Estimation")
lines(new_x, predictions, col = "blue", lwd = 2)
lines(new_x, sin(2 * pi * new_x) + 0.5 * new_x, col = "red", lty = 2, lwd = 2) # 真の関数
legend("topleft", legend = c("Observed Data", "GAM Prediction", "True Function"),
       col = c("black", "blue", "red"), lty = c(1, 1, 2), lwd = 2)

上記の例では、s(x) が x に対する平滑化関数を指定しています。gam() 関数は、データに基づいて最適な平滑性を自動的に推定します。summary(model) を実行すると、推定された平滑化関数の有効自由度（edf: estimated degrees of freedom）などが表示され、モデルがどの程度柔軟な関数を学習したかを知ることができます。

• 自動的な平滑性の調整
  統合的平滑性推定により、過学習のリスクを軽減しつつ、データに適切な複雑さのモデルを構築できます。
• 解釈可能性
  個々の平滑化関数の形状を視覚的に解釈することができます。
• データの探索的分析
  予測変数と応答変数の間の潜在的なパターンを発見するのに役立ちます。
• 非線形性の柔軟なモデリング
  線形モデルでは捉えられない複雑な関係性を捉えることができます。

---

R における GAM (統合的平滑性推定付き) の一般的なエラーとトラブルシューティング

mgcv パッケージの gam() 関数を使用する際に遭遇しやすいエラーと、その解決策について解説します。

モデル収束に関するエラー

• トラブルシューティング
  

  • 反復回数を増やす
    gam() 関数の control 引数で maxit をより大きな値に設定します。

    model <- gam(y ~ s(x), data = data, control = list(maxit = 200))
    

  • 収束判定基準を緩める
    control 引数で epsilon (収束の許容誤差) をより大きな値に設定します。ただし、これは最終的な解の精度に影響を与える可能性があるため、注意が必要です。

    model <- gam(y ~ s(x), data = data, control = list(epsilon = 1e-6)) # デフォルトは 1e-7
    

  • 異なる最適化アルゴリズムを試す
    method 引数を変更してみます ("GCV.Cp", "REML", "P-ML" など)。データやモデルの特性によって、最適なアルゴリズムが異なる場合があります。

    model <- gam(y ~ s(x), data = data, method = "REML")
    

  • より単純なモデルを試す
    平滑化関数の複雑さを減らす（例：k の値を小さくする）か、より少ない予測変数を使用してみます。
  • データのスケーリング
    予測変数のスケールが大きく異なる場合、収束を妨げる可能性があります。変数を標準化または中心化することを検討してください。

• 原因
  モデルの最適化アルゴリズムが、指定された反復回数内で収束しなかったことを示しています。これは、モデルが複雑すぎる、データがスパースである、初期値が不適切である、などが原因で起こり得ます。

• • Iteration limit reached.
  • Model convergence problem. Iterations stopped.
  • step halving failed to reduce deviance in smoothing parameter estimation

平滑化パラメータ推定に関するエラー

• トラブルシューティング
  

  • 平滑化基底の次元 (k) の調整
    s() 関数の k 引数は、平滑化関数の柔軟性を制御します。k が小さすぎるとモデルがデータに十分に適合せず、大きすぎると過学習のリスクが高まります。データの複雑さに合わせて k を調整してみてください。

    model <- gam(y ~ s(x, k = 10), data = data) # k の値を調整
    

  • 罰則の種類の変更
    s() 関数の bs 引数で異なる基底関数（罰則の種類）を試してみます。例えば、"tp" (thin plate regression spline) がデフォルトですが、他のオプション ("cr", "ps", "ts") も検討できます。

    model <- gam(y ~ s(x, bs = "cr"), data = data)
    

  • ロバストな回帰
    外れ値の影響が大きい場合は、family 引数でロバストな分布（例：gaussian(link = "identity", ...) に weights を指定する、または MASS パッケージの関数を使用する）を検討します。
  • データの確認
    欠損値、極端な外れ値、共線性などの問題がないかデータを確認し、必要に応じて前処理を行います。

• 原因
  平滑化パラメータの最適化プロセスで問題が発生しています。これは、モデルの特定化が不適切である、データに極端な値が含まれている、などが原因となることがあります。

• エラーメッセージの例
  

  • A problem occurred in the penalized likelihood calculation.
  • Inner iteration failed.

モデルの特定化に関するエラー

• トラブルシューティング
  

  • formula の確認
    モデルの formula が正しく記述されているか、応答変数と予測変数の名前がデータフレーム内の変数名と一致しているかを確認します。
  • データフレームの確認
    モデルで使用する変数が含まれているデータフレームが正しく指定されているか (data 引数) を確認します。
  • 変数の存在確認
    使用しようとしている変数が R の環境内で定義されているかを確認します。

• 原因
  モデルの formula の指定が誤っているか、必要な変数がデータフレーム内に存在しない場合に発生します。

• エラーメッセージの例
  

  • Error in model.frame.default(formula = y ~ s(x), data = data) : variable 'y' not found (変数がデータフレームに存在しない)
  • Error in eval(predvars, data, env) : object 'x' not found (変数がスコープ内に見つからない)

• リンク関数の選択
  GLM の枠組みで GAM を使用する場合、応答変数の分布に適したリンク関数を選択する必要があります。不適切なリンク関数は、モデルの収束不良や不適切な当てはまりを引き起こす可能性があります。
• 自由度の超過
  平滑化関数の自由度 (k の値) がデータポイント数に近い場合、過学習のリスクが高まります。k の値を適切に設定することが重要です。summary(model) で推定された有効自由度 (edf) を確認し、必要に応じて k を調整します。
• 多重共線性
  複数の予測変数間に強い相関がある場合、モデルの推定が不安定になることがあります。相関の高い変数を削除したり、次元削減の手法を検討したりします。
• NA 値の取り扱い
  データに NA 値が含まれている場合、gam() 関数はデフォルトでそれらを無視します。必要に応じて、na.action 引数を変更したり、事前に NA 値を処理したりする必要があります。

トラブルシューティングの一般的なヒント

• パッケージのバージョンを確認
  古いバージョンのパッケージを使用している場合、既知のバグに遭遇する可能性があります。update.packages("mgcv") でパッケージを最新の状態に更新してみてください。
• 残差分析
  plot(model) を実行して残差プロットを確認し、モデルの仮定からの逸脱がないか評価します。
• 再現可能な最小限のコードを作成する
  問題を特定しやすくするために、エラーを再現できる最小限のデータとコードを作成して試してみてください。
• エラーメッセージをよく読む
  エラーメッセージは、問題の原因を示唆する重要な情報を含んでいます。

---

R における GAM (統合的平滑性推定付き) のプログラミング例

ここでは、mgcv パッケージの gam() 関数を用いた具体的なコード例を通して、GAM の基本的な使い方から、より高度な応用までを解説します。

基本的な GAM の適用 (単一の平滑化項)

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

# サンプルデータの作成
set.seed(123)
n <- 100
x <- sort(runif(n))
true_y <- sin(2 * pi * x) + 0.5 * x
y <- true_y + rnorm(n, sd = 0.2)
data <- data.frame(x = x, y = y)

# 基本的な GAM のモデル fitting
# y を x の平滑化関数 s(x) でモデル化
model1 <- gam(y ~ s(x), data = data)

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

# 予測
new_x <- seq(0, 1, length.out = 100)
predictions1 <- predict(model1, newdata = data.frame(x = new_x))

# 結果のプロット
plot(data$x, data$y, main = "Basic GAM with s(x)")
lines(new_x, predictions1, col = "blue", lwd = 2)
lines(new_x, true_y[order(x)], col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("Observed Data", "GAM Prediction", "True Function"),
       col = c("black", "blue", "red"), lty = c(1, 1, 2), lwd = 2)

この例では、応答変数 y を単一の予測変数 x の滑らかな関数 s(x) でモデル化しています。s() 関数は、thin plate regression spline をデフォルトの平滑化関数として使用し、データのパターンを柔軟に捉えます。

複数の平滑化項の追加

# 複数の予測変数を持つサンプルデータの作成
set.seed(456)
n <- 150
x1 <- sort(runif(n))
x2 <- runif(n)
true_y2 <- sin(pi * x1) + 2 * (x2 - 0.5)^2
y2 <- true_y2 + rnorm(n, sd = 0.3)
data2 <- data.frame(x1 = x1, x2 = x2, y = y2)

# 複数の平滑化項を持つ GAM のモデル fitting
# y を x1 の平滑化関数 s(x1) と x2 の平滑化関数 s(x2) でモデル化
model2 <- gam(y ~ s(x1) + s(x2), data = data2)

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

# 予測 (ここでは x1 と x2 のグリッドを作成して予測)
x1_grid <- seq(0, 1, length.out = 30)
x2_grid <- seq(0, 1, length.out = 30)
new_data2 <- expand.grid(x1 = x1_grid, x2 = x2_grid)
predictions2 <- predict(model2, newdata = new_data2)
prediction_matrix <- matrix(predictions2, nrow = 30, ncol = 30)

# 結果の可視化 (等高線プロット)
contour(x1_grid, x2_grid, prediction_matrix, main = "GAM with s(x1) + s(x2)")
points(data2$x1, data2$x2, pch = 16, cex = 0.5)

この例では、二つの予測変数 x1 と x2 が応答変数 y に与える非線形な影響を、それぞれ s(x1) と s(x2) でモデル化しています。

線形項との組み合わせ

# 線形効果を持つ予測変数を含むサンプルデータの作成
set.seed(789)
n <- 120
x3 <- sort(runif(n))
z <- 2 * x3 + 1
true_y3 <- sin(2 * pi * x3) + z
y3 <- true_y3 + rnorm(n, sd = 0.4)
data3 <- data.frame(x3 = x3, z = z, y = y3)

# 平滑化項と線形項を組み合わせた GAM のモデル fitting
# y を x3 の平滑化関数 s(x3) と z の線形効果でモデル化
model3 <- gam(y ~ s(x3) + z, data = data3)

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

# 予測
new_x3 <- seq(0, 1, length.out = 100)
new_z <- 2 * new_x3 + 1
predictions3 <- predict(model3, newdata = data.frame(x3 = new_x3, z = new_z))

# 結果のプロット
plot(data3$x3, data3$y, main = "GAM with s(x3) + linear(z)")
lines(new_x3, predictions3, col = "blue", lwd = 2)
lines(new_x3, true_y3[order(x3)], col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("Observed Data", "GAM Prediction", "True Function"),
       col = c("black", "blue", "red"), lty = c(1, 1, 2), lwd = 2)

ここでは、非線形な影響を持つ x3 と、線形な影響を持つ z の両方をモデルに含めています。z は通常の線形項として扱われます。

交互作用項のモデリング

# 二つの予測変数の交互作用を含むサンプルデータの作成
set.seed(101)
n <- 100
u <- runif(n)
v <- runif(n)
true_y4 <- (2 * u - 1) * (2 * v - 1)
y4 <- true_y4 + rnorm(n, sd = 0.1)
data4 <- data.frame(u = u, v = v, y = y4)

# テンソル積平滑化 (te) を用いた交互作用のモデリング
model4 <- gam(y ~ te(u, v), data = data4)

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

# 予測 (グリッドデータで予測)
u_grid <- seq(0, 1, length.out = 20)
v_grid <- seq(0, 1, length.out = 20)
new_data4 <- expand.grid(u = u_grid, v = v_grid)
predictions4 <- predict(model4, newdata = new_data4)
prediction_matrix4 <- matrix(predictions4, nrow = 20, ncol = 20)

# 結果の可視化 (3D プロット)
library(rgl)
persp3d(u_grid, v_grid, prediction_matrix4, col = "lightblue",
        xlab = "u", ylab = "v", zlab = "y", main = "GAM with te(u, v)")
points3d(data4$u, data4$v, data4$y, col = "red", size = 3)

te() 関数は、複数の予測変数間の非線形な交互作用をモデル化するために使用されます。ここでは、u と v の交互作用が y に与える影響を捉えています。

グループ化されたデータのモデリング (ランダム効果)

# グループ化されたサンプルデータの作成
set.seed(112)
n_group <- 10
n_obs_per_group <- 20
group <- factor(rep(1:n_group, each = n_obs_per_group))
x5 <- runif(n_group * n_obs_per_group)
random_effects <- rnorm(n_group, sd = 0.5)
true_y5 <- sin(pi * x5) + random_effects[as.numeric(group)]
y5 <- true_y5 + rnorm(n_group * n_obs_per_group, sd = 0.2)
data5 <- data.frame(x = x5, y = y5, group = group)

# ランダム効果を含む GAM のモデル fitting
# s(group, bs = "re") はグループごとのランダム切片を表す
model5 <- gam(y ~ s(x) + s(group, bs = "re"), data = data5)

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

# グループごとのランダム効果の推定値
ranef(model5)

s(group, bs = "re") は、グループ group ごとのランダムな切片をモデルに組み込むための指定です。これにより、グループ間のベースラインの違いを考慮したモデリングが可能になります。

一般化加法モデル (GLM の枠組みでの GAM)

# 二値応答変数のサンプルデータの作成 (ロジスティック回帰の例)
set.seed(131)
n <- 100
x6 <- sort(runif(n))
linear_predictor <- 2 + 5 * x6 + 3 * sin(2 * pi * x6)
probabilities <- 1 / (1 + exp(-linear_predictor))
y6 <- rbinom(n, 1, probabilities)
data6 <- data.frame(x = x6, y = y6)

# ロジスティック GAM のモデル fitting
# family = binomial は二項分布 (ロジスティック回帰) を指定
model6 <- gam(y ~ s(x), family = binomial, data = data6)

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

# 予測 (確率)
new_x6 <- seq(0, 1, length.out = 100)
probabilities_predicted <- predict(model6, newdata = data.frame(x = new_x6), type = "response")

# 結果のプロット
plot(data6$x, data6$y, main = "Logistic GAM")
lines(new_x6, probabilities_predicted, col = "blue", lwd = 2)
lines(new_x6, probabilities, col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("Observed (0/1)", "Predicted Probability", "True Probability"),
       col = c("black", "blue", "red"), lty = c(1, 1, 2), lwd = 2)

family 引数に適切な分布とリンク関数を指定することで、二値データ（ロジスティック回帰）、カウントデータ（ポアソン回帰など）、ガンマ分布に従うデータなど、さまざまな種類の応答変数を GAM でモデル化できます。

平滑化関数の基底と次元の指定 (k と bs)

# 平滑化関数の基底 (bs) と次元 (k) を指定する例
set.seed(141)
n <- 100
x7 <- sort(runif(n))
true_y7 <- exp(-5 * (x7 - 0.5)^2)
y7 <- true_y7 + rnorm(n, sd = 0.05)
data7 <- data.frame(x = x7, y = y7)

# cubic regression spline (bs = "cr") を使用し、基底次元を k = 8 に設定
model7 <- gam(y ~ s(x, bs = "cr", k = 8), data = data7)

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

# 予測とプロット (省略)

bs 引数で平滑化関数の種類（例："tp" (thin plate), "cr" (cubic regression), "ps" (P-spline) など）を、k 引数で基底関数の次元（自由度の上限）を指定できます。k は、モデルの柔軟性に影響を与える重要なパラメータです。

---

以下に、いくつかの代替となり得る方法や関連するパッケージを紹介しますが、それぞれに特徴や適用範囲の違いがあることに注意してください。

gam パッケージ (オリジナルの GAM 実装)

• 注意点
  gam パッケージは、mgcv ほど活発に開発されておらず、最新の GAM の手法や機能が実装されていない場合があります。
• 柔軟性
  mgcv ほど多様な平滑化関数やモデルの構造をサポートしていません。
• 平滑性推定
  平滑化パラメータは、通常、交差検定 (cross-validation) などの外部的な方法によって選択されます。mgcv のようにモデル推定プロセスに統合されていません。
• 特徴
  Trevor Hastie 氏によるオリジナルの gam パッケージは、GAM の基本的な枠組みを提供します。mgcv パッケージが罰則付き尤度最大化と自動的な平滑性推定に重点を置いているのに対し、gam パッケージは主に backfitting アルゴリズムを用いてモデルを推定します。

gamlss パッケージ (一般化加法位置・尺度・形状モデル)

• 注意点
  gamlss はより複雑なモデルを扱うことができる反面、学習コストが高く、mgcv とは異なるモデリング哲学に基づいています。
• 柔軟性
  非常に柔軟な分布とリンク関数の選択肢を提供し、さまざまな種類の応答変数に対応できます。
• 平滑性推定
  gamlss でも平滑化関数 (pb() 関数など) を使用できますが、平滑性推定の方法や統合の度合いは mgcv と異なります。
• 特徴
  gamlss パッケージは、応答変数の平均だけでなく、分散や形状などの分布パラメータも予測変数の関数としてモデル化できる、より広範なフレームワークを提供します。GAM もその特殊なケースとして扱うことができます。

機械学習の手法 (柔軟な非線形モデリング)

• 注意点
  機械学習モデルは、統計的な推論や平滑化された関数の解釈を主な目的としない場合に使用されます。
• 例 (Random Forest)
  

  # randomForest パッケージのインストール
  install.packages("randomForest")
  library(randomForest)
  
  # サンプルデータの作成 (mgcv の例と同じ)
  set.seed(123)
  n <- 100
  x <- sort(runif(n))
  true_y <- sin(2 * pi * x) + 0.5 * x
  y <- true_y + rnorm(n, sd = 0.2)
  data <- data.frame(x = x, y = y)
  
  # Random Forest モデルの fitting
  model_rf <- randomForest(y ~ x, data = data)
  
  # 予測
  new_x <- seq(0, 1, length.out = 100)
  predictions_rf <- predict(model_rf, newdata = data.frame(x = new_x))
  
  # 結果のプロット (省略)
  

• 柔軟性
  非常に高い柔軟性を持ち、複雑なデータ構造に対応できますが、モデルの解釈が難しい場合があります (ブラックボックスモデル)。
• 平滑性推定
  これらの手法は、明示的な平滑化関数の概念を持たず、データのパターンを学習する過程で暗黙的に複雑さを制御します (正則化など)。統合的平滑性推定とは異なるアプローチです。
• 特徴
  Random Forest、Gradient Boosting Machines (GBM)、ニューラルネットワークなどの機械学習アルゴリズムは、予測変数と応答変数の間の複雑な非線形関係を捉えることができます。

スプライン回帰のための他のパッケージ

• earth パッケージ
  Multivariate Adaptive Regression Splines (MARS) を実装しており、データ駆動型で区分的な線形関数を組み合わせることで非線形性をモデル化します。自動的な変数選択や非線形性の検出に役立ちますが、伝統的な GAM の平滑化関数とは異なります。
• splines パッケージ
  自然スプライン (ns()) や B スプライン (bs()) などの基本的なスプライン関数を提供します。これらの関数を線形回帰モデルや GLM の予測変数として使用することで、非線形な関係をモデル化できますが、平滑化パラメータの選択は手動で行うか、交差検定などの方法を用いる必要があります。統合的平滑性推定は行われません。

統合的平滑性推定を伴う GAM を行うための最も強力で柔軟なツールは、依然として mgcv パッケージの gam() 関数です。他のパッケージや手法は、それぞれ異なる特徴や目的を持っており、特定の状況下では代替となり得ますが、GAM の核となる機能 (柔軟な非線形モデリングと自動的な平滑性推定) を完全に代替するわけではありません。