GLMをもう少し理解したい④
前回の記事では、結局GLMというのは以下の方程式:
を用いて、を反復的に求めることであると説明しました(IRLS)。
そのために必要なパーツとしてはとであり、これらは(を除けば)、とそれらのそれぞれに対する微分です。 ではそれらの値をどのように求めるのか実際に試してみましょう。
ポアソン回帰
関数の定義
まずはポアソン回帰で確かめてみましょう。以下のように関数を定義します:
get_eta <- function(b) { eta <- X %*% b return(eta) } get_mu <- function(b) { mu <- exp(X %*% b) return(mu) } get_var <- get_mu
GLMにおいては常に線形予測子によって説明されるものなので、シンプルにとの積によって得られます。一方はを逆リンク関数での世界に戻したものなので、ポアソン回帰で一般的な逆リンク関数である exp
を使います。
3つ目の関数はの分散を得るためのものですが、ポアソン回帰ではを仮定しますので get_mu
を再利用しています。
続いてとを得るための関数を定義します:
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
を用いました。これは前回の記事で説明した通りは
で得られますが1、
となるためです。また get_W
の分子には get_mu
の二乗を使いましたが、これは
となるためです。
サンプルデータ作成
続いてサンプルデータを作成します。
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 = 1
、 b2 = 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においては対角行列であり、対角要素のみが更新されるため予め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.02039846
、 b2 = 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
を用いました。
とについては、
および
を使っています。
サンプルデータ作成
先程と同様にサンプルデータを作成します。
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
を使うことができそうですね。
-
一部修正あり↩