Rと手作業で覚えるGLM - 餡子付゛録゛

glmコマンドに頼り切って二項分布（と言うかロジスティック回帰）やポアソン回帰をしている日々ですが、速度狂は最初から何でもC++で･･･と言う話がある*1ので一般化線形回帰（GLM: Generalized Linear Model）の手順を少し詳しく見てみましょう。最尤法を使っても良いのですが、それではGLMにならないので。

1. GLMとは

GLMの概要については『統計モデル入門』の第3章と第4章と関連した付録AとBを読めば把握できますが、一応、概要を説明します。

GLMは、指数分布族の確率分布のモデルに対して使うことができる推定方法です。パラメーター $\theta$ のときに観測値 $y$ が観察される確率 $f(y;\theta)$ が、 $\exp(A(y)B(θ) + C(θ) + D(y))$ の形式で書けるとき、 $f(y;\theta)$ は指数分布族の関数になります。ポアソン分布、正規分布、二項分布などが含まれ、これらは $A(y) = y$ となるので、特に正準形と呼ばれます。

さて、対数尤度関数 $\mathcal{l}(y; \beta)$ が与えられたときの推定を考えましょう。最適化アルゴリズムは色々とありますが、もっとも基本的なものはニュートン・ラフソン法になります。真のパラメーターを $\beta$ 、推定量を $b^{(m)}$ （ $m$ は反復回数で指数ではない）とすると、

$b^{(m)} = b^{(m-1)} - \frac{\partial \mathcal{l}^{(m-1)}}{\partial{b^{(m-1)}}\partial{b^{(m-1)}}} \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}$

となります。ニュートン・ラフソン法は収束が速いのですが、いたるところで2階導関数になるヘッシアンが凸を示す値を返さないと、局所最適解にはまりがちです。準ニュートン法と呼ばれる発展的な最適化アルゴリズムは、この2階導関数の計算を省く工夫していますが、GLMも同様です。

GLMの場合は、2階導関数を、ヘッシアンにマイナスをかけた値になるフィッシャーの情報行列（information matrix） $\mathcal{I}$ で置き換えます（スコア推定）。

$\mathcal{I}^{(m-1)}b^{(m)} = \mathcal{I}^{(m-1)} b^{(m-1)} + \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}$

$\mathcal{I}$ は、分散共分散行列の符号をマイナスにしたものになります*2。情報行列をどうやって計算すればいいのかが問題になりますが、上述の指数分布族の特性と対数尤度関数の一般的特性を前提に式を整理していくと、 $\mathcal{I}$ のi行j列の要素を、

$\mathcal{I}_{ij} = \sum^N_{k=1} \frac{x_{ki}x_{kj}}{var(y_k)} \bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)^2$

と書くことができます。 $N$ はサンプルサイズです。ここで $\mu$ は観測値 $y$ の期待値、 $\eta$ は説明変数 $X$ と係数 $\beta$ の線形結合 $X\beta$ で、 $\mu$ と $\eta$ はリンク関数 $g$ で対応が定義されます。

$g(\mu) = X\beta = \eta$

具体的な導出は、上述の参考文献にざっくり説明*3があるので、そちらを参照してください。推定は以上と具体的なリンク関数を整理することでできますが、もう少し式を整理すると、一般化線形回帰と言う名称の意味が分かります。

以下の要素の対角行列 $W$ と、

$w_{ii} = \frac{1}{var(y_k)} \bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)^2$

以下を要素に持つベクトル $z$ を定義し、

$z_i = \sum_{k=1} x_{ik} b_k^{m-1} + (y_i - \mu_i)\bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)$

スコア法の反復式を整理すると、

$X^t W X b^m = X^t W z$

となり、ウェイト付線形回帰の形式に持っていくことができます。 $W$ も $z$ も $b$ に依存するので有り難味は微妙なのですが。

2. 実践

さて、本題に入り、べたべたとRのコードを貼っていきましょう。

2.1. ポアソン回帰

ポアソン分布の平均と分散になるλの推定は、ウェイト付線形回帰の形式で実行するとそれっぽいです。

# データセットを作成
set.seed(2022)
n <- 50 # サンプルサイズ
X <- matrix(c(rep(1, n), (1:n) %% 3), n, 2)
beta <- matrix(c(1, 2), 2, 1) # 観察できない真のパラメーター
lambda <- X %*% beta
y <- sapply(1:n, function(i){
  rpois(1, lambda=lambda[i])
})
# 推定してみる
b <- matrix(c(0, 0), 2, 1) # 初期値は(0, 0)から
for(i in 1:30){
  W <- diag(c(1/(X %*% beta)))
  z <- (X %*% beta) + (y - X %*% beta)
  L <- t(X) %*% W %*% X
  R <- t(X) %*% W %*% z
  b <- solve(L, R)
}
# 推定結果
b
# 真のパラメーターβと比較
round(b - beta, 6)
# ビルトイン関数のglmでも推定
r_glm = glm(y ~ 0 + X, family = poisson(link = "identity")) # デフォルトはlink = "log"になっているのだが、ここではlink = "identity"
# glmの結果との結果と比較
round(coef(r_glm) - b, 6)

だいたい同じ結果になります。

2.2. ロジスティック回帰

ロジスティック回帰のためだけにGLMを使っている人も多いとか何とか･･･

# データセットを作成
set.seed(2022)
n <- 100 # サンプルサイズ
# 観察できるX
X <- matrix(c(rep(1, n), runif(n, min=0, max=10)), n, 2)
# 観察できない推定パラメーター
beta <- matrix(c(-2, 1), 2, 1)
log_odds <- X %*% beta
pi <- exp(log_odds)/(1 + exp(log_odds))
# 観察できるy
y <- 0 + (runif(n) < pi)*1
# 推定してみる
b <- matrix(0, 2, 1) # 初期値
for(i in 1:30){
  t_lodds <- X %*% b
  t_pi <- exp(t_lodds)/(1 + exp(t_lodds))
  # グラディエント∂l/∂bを計算
  U <- sapply(1:ncol(X), function(i){
    sum(X[,i]*(y - t_pi))
  })

  I <- matrix(0, ncol(X), ncol(X))
  for(i in 1:nrow(I)){
    for(j in 1:ncol(I)){
    # log(π/(1-π)) = Xβ = ηから、せっせと微分して整理すると、∂μ/∂η = π(1-π)になる
    # var(y_k)は、平均と分散が等しいポアソン分布のような場合では無い場合、1と置いて数値演算してよい模様。結局、スコア推定の収束条件は同じになる。
      I[i,j] <- sum(X[,i]*X[,j]*t_pi*(1-t_pi))
    }
  }

  R <- I %*% b + U
  b <- solve(I, R)
}
# ビルトイン関数のglmでも推定
r_glm = glm(y ~ 0 + X, family=binomial(link="logit"))
# glmの結果との結果と比較
round(coef(r_glm) - b, 6)

*1: 無い。

*2:最尤法で係数の標準誤差として、ヘッシアンの逆行列の対角成分を見る理由を思い出しましょう。

*3: ルベーグ積分を未学習の人は、微分と有限ではない区間の積分の交換をしているのが気になるかもです。累積分布関数なので、ルベーグの収束定理を使ってよいわけですが。