一般化加法モデル (GAM)

概要¶

一般化加法モデル（Generalized Additive Models, GAM） は、Hastie と Tibshirani が 1986 年の論文

Hastie, T., & Tibshirani, R. (1986). Generalized additive models. Statistical Science, 1(3), 297–310.

で提案した統計モデル。

一般化線形モデル（GLM）は応答変数の平均 $\mu$ と線形予測子 $\eta = \mathbf{x}^\top \boldsymbol{\beta}$ をリンク関数 $g$ で結ぶ。しかし線形予測子は各特徴量の効果が線形であると仮定しており、非線形な関係を柔軟に表現できない。

GAM は線形予測子を 滑らかな非線形関数の和（加法的構造） に置き換えることで、解釈可能性を保ちながら非線形性を扱う。

各説明変数の効果を個別に可視化できる（部分依存プロットが自然に得られる）
非線形効果を柔軟にモデリングできる
ブラックボックスな機械学習手法とは異なり、各変数の貢献が明示的
- → 交互作用項を入れたGA $^2$ M、Explainable Boosting Machineといった説明性の高い機械学習手法へ発展していった

Hastie & Tibshirani (1986) の論文は GAM の枠組みを整備し、バックフィッティング（backfitting）アルゴリズムという実用的な推定法を提案した点で重要である。

モデル¶

加法モデル¶

応答変数 $Y$ と説明変数 $X_1, \dots, X_p$ について、加法モデル（additive model）は

Y = \alpha + \sum_{j=1}^p f_j(X_j) + \varepsilon

(1)

と表される。ここで

$\alpha$ は切片
$f_j : \mathbb{R} \to \mathbb{R}$ は各説明変数に対する 非線形の滑らかな関数（スムーザー）
$\varepsilon$ は平均ゼロの誤差項

である。識別可能性（identifiability）のため、各 $f_j$ には $\mathbb{E}[f_j(X_j)] = 0$ という中心化の制約を課す。

線形モデルは $f_j(X_j) = \beta_j X_j$ の特殊ケースである。

一般化加法モデル¶

Hastie & Tibshirani (1986) はこれを GLM の枠組みへ拡張した。応答変数の条件付き平均 $\mu = \mathbb{E}[Y \mid X_1, \dots, X_p]$ に対して、リンク関数 $g$ を用いて

g(\mu) = \alpha + f_1(X_1) + f_2(X_2) + \cdots + f_p(X_p)

(2)

とモデル化する。

応答変数の分布	リンク関数 $g$	モデル
正規分布	$g(\mu) = \mu$ （恒等写像）	加法モデル
ベルヌーイ分布	$g(\mu) = \log\frac{\mu}{1-\mu}$ （ロジット）	加法ロジスティック回帰
ポアソン分布	$g(\mu) = \log(\mu)$ （対数）	加法ポアソン回帰

滑らか関数 $f_j$ の表現¶

滑らか関数 $f_j$ の具体的な形は事前に固定せず、データから推定する。代表的な選択肢は：

スプライン（cubic spline）: 結節点（knot）を設けた区分多項式
LOESS / 局所多項式回帰: 局所加重回帰（→ LOWESSのノートも参照）
カーネル回帰: 重み関数を用いた局所平均
自然スプライン（natural spline）

各 $f_j$ の滑らかさは、平滑化パラメータ（スパン、バンド幅、自由度）で制御する。

アルゴリズム：バックフィッティング¶

問題設定¶

$n$ 個のサンプル $(\mathbf{x}_i, y_i)$ 、 $i = 1, \dots, n$ において、 $\mathbf{x}_i = (x_{i1}, \dots, x_{ip})$ とする。

加法モデルの目的関数（正規誤差の場合）は

\min_{\alpha, f_1, \dots, f_p} \sum_{i=1}^n \left( y_i - \alpha - \sum_{j=1}^p f_j(x_{ij}) \right)^2 + \sum_{j=1}^p \lambda_j \int \{f_j''(t)\}^2 \, dt

(3)

第2項は $f_j$ の曲率に対するペナルティで、 $\lambda_j \geq 0$ は滑らかさを制御するパラメータである。

バックフィッティングアルゴリズム¶

Hastie & Tibshirani (1986) が提案したバックフィッティングは、各 $f_j$ を交互に更新する座標降下法的アプローチである。

アイデア: $j$ 番目の関数 $f_j$ を推定する際、他の関数の寄与を除いた 部分残差（partial residual） に対してスムーザーを適用する。

ここで $\mathcal{S}_j[\cdot]$ は $j$ 番目の変数に対するスムーザーである。すなわち、入力 $\{x_{ij}\}$ と目標値 $\{r_{ij}\}$ （部分残差）のペアに対してスムージングを行い、各点での推定値 $\hat{f}_j(x_{ij})$ を返す写像である。

収束¶

Buja, Hastie & Tibshirani (1989) により、スムーザーが線形かつ対称半正定値である（例：スプラインスムーザー）場合、バックフィッティングは一意な解に収束することが示されている。

収束の判定は、前のイテレーションからの関数の変化量

\frac{1}{n} \sum_{i=1}^n \left( \hat{f}_j^{\text{new}}(x_{ij}) - \hat{f}_j^{\text{old}}(x_{ij}) \right)^2

(8)

が十分小さくなったときに行う。

GAM における IRLS との組み合わせ¶

非正規分布（ロジスティック、ポアソンなど）の場合、GLM の 反復再重み付き最小二乗法（IRLS: Iteratively Reweighted Least Squares） とバックフィッティングを組み合わせる。

具体的には、各 IRLS ステップで調整済み応答変数（adjusted dependent variable）

z_i = \hat{\eta}_i + (y_i - \hat{\mu}_i) \cdot g'(\hat{\mu}_i)

(9)

を計算し、これに対して重み $w_i = 1 / \{V(\hat{\mu}_i) \cdot [g'(\hat{\mu}_i)]^2\}$ 付きのバックフィッティングを適用する。ここで $V(\mu)$ は分散関数である。

実装¶

以下では正規分布・恒等リンクの加法モデルをバックフィッティングで実装する。スムーザーには LOESS（局所多項式回帰）を用いる。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib_fontja
from statsmodels.nonparametric.smoothers_lowess import lowess

np.random.seed(42)

データ生成¶

真のモデルを

y = \sin(x_1) + x_2^2 / 5 + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, 0.3^2)

(10)

として、 $n=300$ のサンプルを生成する。

n = 300
x1 = np.random.uniform(-np.pi, np.pi, n)
x2 = np.random.uniform(-3, 3, n)

f1_true = np.sin(x1)
f2_true = x2**2 / 5
y = f1_true + f2_true + np.random.normal(0, 0.3, n)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(x1, y, alpha=0.3, s=10)
axes[0].set(xlabel="$x_1$", ylabel="$y$", title="$y$ vs $x_1$")
axes[1].scatter(x2, y, alpha=0.3, s=10)
axes[1].set(xlabel="$x_2$", ylabel="$y$", title="$y$ vs $x_2$")
plt.tight_layout()
plt.show()

スムーザーの定義¶

バックフィッティングで使用するスムーザーを定義する。ここでは LOESS を用いる。frac パラメータがバンド幅に相当し、滑らかさを制御する。

def smoother(x: np.ndarray, y: np.ndarray, frac: float = 0.4) -> np.ndarray:
    """LOESS スムーザー: x の順序に沿った推定値を返す"""
    result = lowess(endog=y, exog=x, frac=frac, return_sorted=True)
    x_sorted, y_fitted = result[:, 0], result[:, 1]
    return np.interp(x, x_sorted, y_fitted)

バックフィッティングの実装¶

Hastie & Tibshirani (1986) のアルゴリズムを実装する。

def backfitting(
    X: np.ndarray,
    y: np.ndarray,
    frac: float = 0.4,
    max_iter: int = 100,
    tol: float = 1e-5,
):
    """
    加法モデルのバックフィッティング推定

    Parameters
    ----------
    X : (n, p) 説明変数行列
    y : (n,) 応答変数
    frac : LOESSのバンド幅パラメータ
    max_iter : 最大イテレーション数
    tol : 収束判定閾値

    Returns
    -------
    alpha : 切片の推定値
    F : (n, p) 各 f_j の推定値
    """
    n, p = X.shape

    # 初期化
    alpha = y.mean()
    F = np.zeros((n, p))  # F[:, j] = f_j(x_{1j}), ..., f_j(x_{nj})

    for iteration in range(max_iter):
        F_old = F.copy()

        for j in range(p):
            # 部分残差: j 番目以外の関数の寄与を除く
            partial_residual = y - alpha - F.sum(axis=1) + F[:, j]

            # スムーザーを適用
            f_j_new = smoother(X[:, j], partial_residual, frac=frac)

            # 中心化（識別可能性）
            f_j_new -= f_j_new.mean()
            F[:, j] = f_j_new

        # 収束判定: 全関数の変化量の平均
        delta = np.mean((F - F_old) ** 2)
        if delta < tol:
            print(f"収束: {iteration + 1} イテレーション")
            break
    else:
        print(f"警告: {max_iter} イテレーションで収束しませんでした")

    return alpha, F

X = np.column_stack([x1, x2])
alpha_hat, F_hat = backfitting(X, y, frac=0.3)

収束: 3 イテレーション

推定結果の可視化¶

各 $\hat{f}_j$ を真の関数と重ねて比較する。

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# f1 の比較
order1 = np.argsort(x1)
axes[0].scatter(x1, y - alpha_hat - F_hat[:, 1], alpha=0.2, s=10, label="部分残差")
axes[0].plot(x1[order1], F_hat[order1, 0], "r-", lw=2, label="$\\hat{f}_1$ (GAM推定)")
axes[0].plot(
    x1[order1],
    (f1_true - f1_true.mean())[order1],
    "k--",
    lw=2,
    label="真の $f_1$",
)
axes[0].set(xlabel="$x_1$", ylabel="$f_1(x_1)$", title="$f_1$ の推定")
axes[0].legend()

# f2 の比較
order2 = np.argsort(x2)
axes[1].scatter(x2, y - alpha_hat - F_hat[:, 0], alpha=0.2, s=10, label="部分残差")
axes[1].plot(x2[order2], F_hat[order2, 1], "r-", lw=2, label="$\\hat{f}_2$ (GAM推定)")
axes[1].plot(
    x2[order2],
    (f2_true - f2_true.mean())[order2],
    "k--",
    lw=2,
    label="真の $f_2$",
)
axes[1].set(xlabel="$x_2$", ylabel="$f_2(x_2)$", title="$f_2$ の推定")
axes[1].legend()

plt.tight_layout()
plt.show()

予測精度の確認¶

y_pred = alpha_hat + F_hat.sum(axis=1)
residuals = y - y_pred

rmse = np.sqrt(np.mean(residuals**2))
ss_res = np.sum(residuals**2)
ss_tot = np.sum((y - y.mean())**2)
r2 = 1 - ss_res / ss_tot

print(f"RMSE : {rmse:.4f}")
print(f"R²   : {r2:.4f}")

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(y, y_pred, alpha=0.3, s=10)
lims = [min(y.min(), y_pred.min()), max(y.max(), y_pred.max())]
axes[0].plot(lims, lims, "k--")
axes[0].set(xlabel="実測値 $y$", ylabel="予測値 $\\hat{y}$", title="実測値 vs 予測値")

axes[1].scatter(y_pred, residuals, alpha=0.3, s=10)
axes[1].axhline(0, color="k", linestyle="--")
axes[1].set(xlabel="予測値 $\\hat{y}$", ylabel="残差", title="残差プロット")

plt.tight_layout()
plt.show()

RMSE : 0.3027
R²   : 0.8951

pygam を使った実装¶

実用上は pygam ライブラリを使うと便利である。スプラインベースの GAM を手軽に実装できる。

try:
    from pygam import LinearGAM, s

    gam = LinearGAM(s(0) + s(1)).fit(X, y)
    print(gam.summary())

    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    titles = ["$f_1(x_1)$", "$f_2(x_2)$"]
    for i, ax in enumerate(axes):
        XX = gam.generate_X_grid(term=i)
        pdep, confi = gam.partial_dependence(term=i, X=XX, width=0.95)
        ax.plot(XX[:, i], pdep, label="pygam 推定")
        ax.fill_between(XX[:, i], confi[:, 0], confi[:, 1], alpha=0.3, label="95% CI")
        ax.set(xlabel=f"$x_{i+1}$", ylabel=titles[i], title=titles[i])
        ax.legend()
    plt.tight_layout()
    plt.show()

except ImportError:
    print("pygam がインストールされていません。`pip install pygam` でインストールできます。")

LinearGAM                                                                                                 
=============================================== ==========================================================
Distribution:                        NormalDist Effective DoF:                                      22.834
Link Function:                     IdentityLink Log Likelihood:                                 -1079.0852
Number of Samples:                          300 AIC:                                             2205.8385
                                                AICc:                                            2210.1406
                                                GCV:                                                0.1054
                                                Scale:                                               0.091
                                                Pseudo R-Squared:                                   0.9037
==========================================================================================================
Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
================================= ==================== ============ ============ ============ ============
s(0)                              [0.6]                20           12.3         1.11e-16     ***         
s(1)                              [0.6]                20           10.6         1.11e-16     ***         
intercept                                              1            0.0          1.11e-16     ***         
==========================================================================================================
Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
         which can cause p-values to appear significant when they are not.

WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
         known smoothing parameters, but when smoothing parameters have been estimated, the p-values
         are typically lower than they should be, meaning that the tests reject the null too readily.
None

/tmp/ipykernel_39153/2112680051.py:5: UserWarning: KNOWN BUG: p-values computed in this summary are likely much smaller than they should be. 
 
Please do not make inferences based on these values! 

Collaborate on a solution, and stay up to date at: 
github.com/dswah/pyGAM/issues/163 

  print(gam.summary())

参考文献¶

Hastie, T., & Tibshirani, R. (1986). Generalized additive models. Statistical Science, 1(3), 297–310.
Buja, A., Hastie, T., & Tibshirani, R. (1989). Linear smoothers and additive models. The Annals of Statistics, 17(2), 453–510.
Hastie, T., & Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall.
Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). CRC Press.