統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

GLMをもう少し理解したい④

前回の記事では、結局GLMというのは以下の方程式:


\mathbf{X^ {T}WXb}^ {(m)} = \mathbf{X^ {T}Wz}

を用いて、 \mathbf{b}を反復的に求めることであると説明しました(IRLS)。

ushi-goroshi.hatenablog.com

そのために必要なパーツとしては \mathbf{W} \mathbf{z}であり、これらは( \mathbf{Y}を除けば) \mu \etaとそれらのそれぞれに対する微分です。 ではそれらの値をどのように求めるのか実際に試してみましょう。

ポアソン回帰

関数の定義

まずはポアソン回帰で確かめてみましょう。以下のように関数を定義します:

get_eta <- function(b) {
   eta <- X %*% b
   return(eta)
}

get_mu <- function(b) {
   mu <- exp(X %*% b)
   return(mu)
}

get_var <- get_mu

GLMにおいて \etaは常に線形予測子によって説明されるものなので、シンプルに \mathbf{X} \mathbf{b}の積によって得られます。一方 \mu \etaを逆リンク関数で \mathbf{Y}の世界に戻したものなので、ポアソン回帰で一般的な逆リンク関数である exp を使います。 3つ目の関数は \mathbf{Y}の分散を得るためのものですが、ポアソン回帰では Var(\mathbf{Y}) = E(\mathbf{Y})を仮定しますので get_mu を再利用しています。

続いて \mathbf{W} \mathbf{z}を得るための関数を定義します:

get_z <- function(b) {
   z <- get_eta(b) + (d$y - get_mu(b)) / get_mu(b)
   return(z)
}

get_W <- function(b) {
   w <- get_mu(b)^2 / get_var(b)
   return(w)
}

get_z の二項目には分母として get_mu を用いました。これは前回の記事で説明した通り \mathbf{z}


z_{i} = \Sigma_{k=1}^ {p} x_{ik}b_{k}^ {(m-1)} + (y_{i} - \mu_{i})( \frac{\partial{\eta_{i}}}{\partial{\mu_{i}}} )

で得られますが1


\frac{\partial{\eta_{i}}}{\partial{\mu_{i}}} = \frac{\partial{log(\mu_{i})}}{\partial{\mu_{i}}} = \frac{1}{\mu_{i}}

となるためです。また get_W の分子には get_mu の二乗を使いましたが、これは


\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}} = \frac{\partial{exp(\eta_{i})}}{\partial{\eta_{i}}} = exp(\eta_{i}) = \mu_{i}

となるためです。

サンプルデータ作成

続いてサンプルデータを作成します。

set.seed(123)
n <- 100
x <- cbind(rep(1, n), runif(n, -1, 1))
b <- c(1, 0.5)
lam <- exp(x %*% b)

d <- data.frame(
   y = rpois(n, lam),
   x = x[, -1])

ここでは(推定されるべき)真の回帰係数として b1 = 1b2 = 0.5 としました。 特に難しいところはないと思いますので、このデータを使って解を求めてみます。

解の推定

以下のように反復の条件を設定します:

iteration <- 0 # イテレータ
b_old <- c(2, 0.3) # 解の初期値
threshold <- 1e-06 # 反復終了の閾値
diff <- 1 # 解の更新前後の差
X <- cbind(rep(1, n), d[, -1]) # 説明変数X
n <- nrow(d) # サンプルサイズ
W <- matrix(0, n, n) # W

最後の W ですが、IRLSにおいて \mathbf{W}は対角行列であり、対角要素のみが更新されるため予め0行列を用意しておきました。

では、反復を開始します:

while (diff > threshold) {

   z <- get_z(b_old) # zを計算
   w <- get_W(b_old) # Wを計算
   diag(W) <- w # Wの対角要素を更新
   
   xwx <- t(X) %*% W %*% X # 左辺を計算
   xwz <- t(X) %*% W %*% z # 右辺を計算
   
   b_new <- solve(xwx) %*% xwz # solveで逆行列を求める
   diff  <- sum(abs(b_old - b_new)) # 解の更新前後の差分
   b_old <- b_new # 解の更新
   
   iteration <- iteration + 1
   cat(sprintf("Iterations: %i, b_New_1: %1.8f, b_New_2: %1.8f \n", 
               iteration, b_new[1], b_new[2]))
   if (iteration > 100) break

}

上を実行すると、下記のような結果が得られます。

Iterations: 1, b_New_1: 1.37767882, b_New_2: 0.37311981 
Iterations: 2, b_New_1: 1.07883183, b_New_2: 0.45722167 
Iterations: 3, b_New_1: 1.02225611, b_New_2: 0.49049195 
Iterations: 4, b_New_1: 1.02040076, b_New_2: 0.49240696 
Iterations: 5, b_New_1: 1.02039846, b_New_2: 0.49241027 
Iterations: 6, b_New_1: 1.02039846, b_New_2: 0.49241027

最終的に b1 = 1.02039846b2 = 0.49241027 で停止したようですが、 glm の解と一致するでしょうか?

> coef(glm(y ~ x, d, family = poisson("log")))
(Intercept)           x 
  1.0203985   0.4924103 

合っていますね!

ロジスティック回帰

続いてロジスティック回帰を試してみます。先に言っておくとこの記事の執筆時点でロジスティック回帰の方は上手く行っていませんので予めご承知おきください。

関数の定義

やるべきことはポアソン回帰と変わりません。まずは関数を以下のように定義します:

get_eta <- function(b) {
   eta <- X %*% b
   return(eta)
}

get_mu <- function(b) {
   mu <- exp(X %*% b) / (1 + exp(X %*% b))
   return(mu)
}

get_var <- function(b) {
   var <- get_mu(b) * (1 - get_mu(b))
   return(var)
}

get_z <- function(b) {
   p <- get_mu(b)
   z <- get_eta(b) + (d$y - p) / (1/p + 1/(1-p))
   return(z)
}

get_W <- function(b) {
   t <- get_mu(b) * (1 - get_mu(b))
   w <- t^2 / get_var(b)
   return(w)
}

eta についてはポアソン回帰と同じですが、逆リンク関数が異なるため mu は定義が違います。ロジスティック回帰における一般的なリンク関数はロジットなので、逆リンク関数にはロジスティック関数を用います。また get_var にはベルヌーイ分布の分散である pq を用いました。

 \mathbf{z} \mathbf{W}については、


logit(\mu) = log(\mu) - log(1-\mu) \\
\frac{\partial{logit(\mu_{i})}}{\partial{\mu_{i}}} = \frac{1}{\mu} + \frac{1}{(1-\mu)}

および


logistic(\eta) = \frac{exp(\eta)}{1 + exp(\eta)}  = \frac{exp(\eta) + exp(\eta)^ {2} - exp(\eta)^ {2}}{(1 + exp(\eta))^ {2}} \\
= \frac{exp(\eta)}{(1 + exp(\eta))^ {2}} = \frac{exp(\eta)}{1 + exp(\eta)} \frac{1}{1 + exp(\eta)} \\
= logistic(\eta) (1 - logistic(\eta))

を使っています。

サンプルデータ作成

先程と同様にサンプルデータを作成します。

set.seed(789)
n <- 100
x <- cbind(rep(1, n), runif(n, -1, 1))
b <- c(1, 0.5)
eta <- x %*% b
p <- exp(eta)/(1 + exp(eta))

d <- data.frame(
   y = rbinom(n, 1, p),
   x = x[, -1])

回帰係数などは特に変更していません。

解の推定

では反復を開始します。ポアソン回帰の時と違い、随分と収束に時間がかかるため反復回数の上限を500回とし、100回ごとに表示しています。また収束の基準も厳し目にしましたが、それ以外は while の中身に変更はありません。

iteration <- 0
b_old <- c(2, 0.3)
threshold <- 1e-08
diff <- 1
X <- cbind(rep(1, n), d[, -1])
n <- nrow(d)
W <- matrix(0, n, n)

while (diff > threshold) {

   z <- get_z(b_old)
   w <- get_W(b_old)
   diag(W) <- w
   
   xwx <- t(X) %*% W %*% X
   xwz <- t(X) %*% W %*% z
   
   b_new <- solve(xwx) %*% xwz
   diff  <- sum(abs(b_old - b_new))
   b_old <- b_new
   
   iteration <- iteration + 1
   if (iteration %% 100 == 0) {
      cat(sprintf("Iterations: %i, b_New_1: %1.8f, b_New_2: %1.8f \n", 
                  iteration, b_new[1], b_new[2]))
   }
   if (iteration > 500) break
}

上記を実行すると、以下のような結果が得られます:

Iterations: 100, b_New_1: 0.91099754, b_New_2: 0.35764440 
Iterations: 200, b_New_1: 0.86993338, b_New_2: 0.33293052 
Iterations: 300, b_New_1: 0.86937860, b_New_2: 0.33250983

一応、500回までは行かずに収束したと判断されたようなので、 b_new を確認してみましょう。

> b_new
          [,1]
[1,] 0.8693712
[2,] 0.3325039

> coef(glm(y ~ x, d, family = binomial("logit")))
(Intercept)           x 
  0.8825275   0.3975239

…うーん、合っていませんね。収束までにも随分と反復していますし、何かおかしいようです。 残念ながらこの理由についてはまだ良くわかっていません。途中の微分の計算が間違っているのか、分散の定義がダメなのか。。。

終わりに

というわけで、GLMではどのように計算が行われているのかを4回に渡って追いかけてきました。 最後はパッとしない結果になってしまいましたが、IRLSがどのように計算を行っているのかを理解できた気がします。これで今までよりももう少し自信を持って glm を使うことができそうですね。


  1. 一部修正あり