「ガウス過程と機械学習」第3章のグラフ（一部）を作図する
※ 3/13 一部修正
いま「ガウス過程と機械学習」を読んでいるのですが、第3章のグラフが非常に面白かったので自分でも作図してみたところちょっと腹落ちするものがありました。せっかくなので記事にしておきます。今回グラフを作成したのは図3.7と3.8（P69、70）の2つです。
ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)
- 作者: 持橋大地,大羽成征
- 出版社/メーカー: 講談社
- 発売日: 2019/03/09
- メディア: 単行本（ソフトカバー）
- この商品を含むブログを見る
順番は前後してしまうのですが、図3.8から作図します。まずは必要な関数を定義しておきましょう。本にしたがい、共分散行列を計算するRBFカーネル（ガウスカーネル）のtheta1とtheta2はともに1にしてあります。
## 1~4までの数列を適当な刻み幅で生成する generate_vector <- function(by, from = 1, to = 4) seq(by, from = from, to = to) ## xからRBFカーネルを使って共分散行列を作成する get_cov_mat <- function(x, theta1 = 1, theta2 = 1) { vec <- x theta1 * exp(-(apply(as.matrix(x), 1, function(x) x-vec))^2/theta2) }
※ 3/13 apply中のfunctionの引数を修正しました
上記の関数を使って1から4までの数列xを生成し、その共分散行列Kから生成されるyとのプロットを描く関数を以下のように定義します。なるべく本のグラフと近づけるために、共分散行列のヒートマップはcorrplotパッケージのcorrplot関数を使い、method = "shade"を指定しました。
## xからKとyを生成し、そのプロットを描く gen_y_plot <- function(by) { x <- generate_vector(by) K <- get_cov_mat(x) mu <- rep(0, length(x)) y <- MASS::mvrnorm(n = 1, mu = mu, Sigma = K) par(mfcol = c(1, 2)) plot(x, y, xlim = c(0, 5)) corrplot::corrplot(K, method = "shade") }
数列の刻み幅を変えながらいくつかプロットしてみます。
set.seed(123) gen_y_plot(1) gen_y_plot(0.5) gen_y_plot(0.2) gen_y_plot(0.1)
本にある通り、刻み幅を小さくすることで順に滑らかな曲線になっていくのがよくわかります。
続いて図3.7を作成します。先ほどとほぼ同じです。
gen_x_plot <- function(by) { x <- generate_vector(by, -4, 4) K <- get_cov_mat(x) mu <- rep(0, length(x)) y <- MASS::mvrnorm(n = 1, mu = mu, Sigma = K) plot(x, y, type = "l", ylim = c(-4, 4)) }
plotを重ねるためにpar(new = T)を使っているのですが、一回目だけは作図を0から始める必要があるので処理を分けています。もっと良い方法がある気がしていますが、ひとまず欲しいグラフは得られたので良しとします。
set.seed(123) for (i in 1:4) { if(i == 1) { gen_x_plot(0.05) next } par(new = T) gen_x_plot(0.05) }
自分で書いてみると理解できているかの確認になりますし、パラメータを変えたときの挙動を見ることで理解が深まりますので本当にオススメです。この後もなるべく自分で作図してみようと思います。
データを小集団に分割しながら線形回帰の解を推定する
背景
突然ですが、一般に線形回帰と言えば以下の正規方程式：
※ 7/3 上の式を修正
をbについて解くことで得られると教わり、そのまま理解していることが多いのではないでしょうか。
これ自体は決して間違っていないのですが、実装においては計算精度の問題から、逆行列ではなくQR分解を活用して解いている場合があります。例えばRでも、以前の記事においてlmのソースコードをたどっていった結果、ハウスホルダー変換によってQR分解が行われていることを確認しました。
過去の記事はこちら。 ushi-goroshi.hatenablog.com
ここでlmと逆行列およびQR分解による解の推定値をちょっと見てみましょう。適当にデータを作成します。
set.seed(123) n <- 100 b <- c(1, 1.5) # 切片と回帰係数 x <- cbind(1, rnorm(n)) y <- x %*% b + rnorm(n)
また、それぞれによる解の推定方法を以下のように定義します。
## lm(.fit)を使う my_lm <- function() { coef(lm.fit(x, y)) } ## 逆行列で解く my_solve <- function() { solve(crossprod(x, x)) %*% crossprod(x, y) } ## QR分解で解く my_qr <- function() { solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y }
上で定義した関数は、いずれも同じ解を返します：
> cbind(my_solve(), my_qr(), my_lm()) [,1] [,2] [,3] x1 0.8971969 0.8971969 0.8971969 x2 1.4475284 1.4475284 1.4475284
一緒の値になっていますね。少し脱線しますが、ついでに計算時間も見てみましょう：
time_1000 <- data.frame(microbenchmark::microbenchmark(my_solve(), my_qr(), my_lm(), times = 1000))
> library(ggplot2) > ggplot(time_1000, aes(x = expr, y = log(time), group = expr)) + geom_violin() + coord_flip() + labs(x = "functions") + NULL
逆行列を用いた場合が一番早く、QR分解を用いたものが最も遅いようでした。なおこのグラフは横軸が対数となっていることに注意してください。
さて、このようにして線形回帰の解はQR分解を使って求めることができますが、実は計算を工夫することで、Xを小集団に分割した上でそれぞれのデータからX全体の解を得ることができます。これが何を意味するかというと、メモリに全部載せきれないような大きいデータであっても解を推定したり、あるいは線形回帰であっても並列に計算を回すことができる、ということです1。
もともと今回の記事を書こうと思ったのは、以前に「線形回帰はデータを分割して並列計算できる」という話を知人から聞いたことをふと思い出したのがきっかけです。当時は何を言っているのか今いち理解できなかったのですが、大変わかりやすい下記の記事を見つけたため、写経した内容をメモしておきます。
freakonometrics.hypotheses.org
手順
実装に取り掛かる前に手順について簡単に理解しておきましょう。まずXをQR分解すると、冒頭に示した正規方程式から得られるは以下のようになります：
QR分解によって得られる行列Qは直交行列であるため、となります。またここで積の逆行列は
という性質があることから、
となります。すなわちQR分解によって得られた行列Rの逆行列と、行列Qの転置があれば良いことになります。先ほどmy_qrを定義したときは説明なく示しましたが、これは下のように書けます：
## my_qrの定義（再掲） solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y
問題は、このおよび
をどのようにして小集団から再構成するか、ということになりますが、これは以下の手順で計算できるようです：
- 共通処理
- X、yをそれぞれ小集団に分割する
- 各小集団のXをQR分解する
を計算する
を計算する
- 1-2で得られたQを2-2で得たQに乗じる（
- 1-2で得られたQを2-2で得たQに乗じる（
- 2と3の結果およびyにより解を得る
- 3-1で得たQ'にyを乗じる
- 両者を乗じる
なおこの手順で確かにや
が再構成できることは確認できたのですが、これがなぜ上手くいくのかについては残念ながら調べても分かりませんでした。もしご存知でしたら誰か教えてください。
実装
それでは実装に入りますが、先にデータをすべて使った時の回帰係数を確認しておきましょう。サンプルデータにはcarsを使い、目的変数をdist、説明変数をspeedとした単回帰を回してみます。
> lm(dist ~ speed, data = cars)$coefficients (Intercept) speed -17.579095 3.932409
切片とspeedの回帰係数がそれぞれ-17.579、3.932と推定されました。冒頭でも確認した通り、lmの結果は下記の方法と一致します。
y <- cars$dist x <- cbind(1, cars$speed)
> cbind( + solve(crossprod(x, x)) %*% crossprod(x, y), + solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y + ) [,1] [,2] [1,] -17.579095 -17.579095 [2,] 3.932409 3.932409
この数値を、分割した小集団に対する計算結果から再び得ることが目標となります。
1. 共通処理
1. X、yをそれぞれ小集団に分割する
それではxを小集団に分割した上で解を推定していきます。今回はデータを5つに分割しましょう。xは50行のデータなので各データセットには10行ずつ割り当てられます。各データをlist形式で保存しておきます。
# 分割するデータの数 m <- 5 n_per_d <- nrow(x) / m # 割り切れなかった場合用 if (nrow(x) %% m != 0) m <- m + 1 xlist <- list() # 各xの保存用リスト ylist <- list() # 各yの保存用リスト for (i in 1:m) { if(i == m) { xlist[[i]] = x[((i-1) * n_per_d + 1):nrow(x), ] ylist[[i]] = y[((i-1) * n_per_d + 1):nrow(x)] } xlist[[i]] = x[(i-1) * n_per_d + 1:n_per_d, ] ylist[[i]] = y[(i-1) * n_per_d + 1:n_per_d] }
このような形でデータが保存されます：
> head(xlist[[1]]) [,1] [,2] [1,] 1 4 [2,] 1 4 [3,] 1 7 [4,] 1 7 [5,] 1 8 [6,] 1 9
2. 各小集団のXをQR分解する
次に各小集団をQR分解し、その結果として得られる行列QおよびRをそれぞれ保存しておきましょう。リストの各要素は、更にそれぞれQとRを要素に持つリストとなります。
QR1 <- list() # 各データセットに対するQR分解の結果を保存するリスト for (i in 1:m) { QR1[[i]] <- list(Q = qr.Q(qr(xlist[[i]])), R = qr.R(qr(xlist[[i]]))) }
この時点でQR1は、10行2列の行列Qと2行2列の上三角行列Rを要素に持つリストになっています。
> str(QR1) List of 5 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -25.3 7.48 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -39.53 2.55 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -48.38 3.48 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -58.82 2.9 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -71.47 5.87
2. R^{-1}を計算する
1. 各小集団からのRを統合する
続いてQR1に格納された行列Rを、rbindで一つにまとめます。
R1 <- c() for(i in 1:m) { R1 <- rbind(R1, QR1[[i]]$R) }
2. 再度QR分解してRを得る
このR1を再度QR分解し、 その行列Rを得ます。
R2 <- qr.R(qr(R1))
3. Rの逆行列を求める（R^{-1}）
この逆行列が、当初求めようとしていたものの1つになります。
R_inv <- solve(R2)
では、このR_invがデータ全体を使って求めたを同じ値になっているかを確認してみましょう。
> R_inv [,1] [,2] [1,] 0.1414214 0.41606428 [2,] 0.0000000 -0.02701716 > solve(qr.R(qr(x))) [,1] [,2] [1,] -0.1414214 -0.41606428 [2,] 0.0000000 0.02701716
あれ？符号が反転していますね。。
3. Q'を計算する
ひとまず置いておいて、先に進みましょう。
1. 1-2で得られたQを2-2で得たQに乗じる（Q'）
先ほどR2を計算したときと同じQR分解で、今度は行列Qを得ます。
Q1 <- qr.Q(qr(R1))
さらに説明変数の数（今回は2）ごとにデータを分割します。
## 説明変数の数 p <- ncol(x) Q2list <- list() for(i in 1:m) { Q2list[[i]] <- Q1[(i-1) * p + 1:p, ] }
このQ2listに、最初にQR分解した結果の行列Q（QR1$Q）を掛け合わせます。
Q3list <- list() for(i in 1:m) { Q3list[[i]] <- QR1[[i]]$Q %*% Q2list[[i]] }
ここで得られたQ3listはデータ全体を使ってQR分解したときのになっているはずです。確認してみましょう：
> head(cbind( + do.call("rbind", Q3list), + qr.Q(qr(x)))) [,1] [,2] [,3] [,4] [1,] 0.1414214 0.3079956 -0.1414214 -0.3079956 [2,] 0.1414214 0.3079956 -0.1414214 -0.3079956 [3,] 0.1414214 0.2269442 -0.1414214 -0.2269442 [4,] 0.1414214 0.2269442 -0.1414214 -0.2269442 [5,] 0.1414214 0.1999270 -0.1414214 -0.1999270 [6,] 0.1414214 0.1729098 -0.1414214 -0.1729098
また符号が反転してますね。。。
原因はわかりませんが、も符号が反転していたので、結果的には元に戻るはずです。気にしないで進めましょう。
4. 2と3の結果およびyにより解を得る
1. 3-1で得たQ'にyを乗じる
上で計算された行列をylistと乗じ、結果を要素ごとに足し合わせます。
Vlist <- list() for(i in 1:m) { Vlist[[i]] <- t(Q3list[[i]]) %*% ylist[[i]] } sumV <- Vlist[[1]] for(i in 2:m) { sumV <- sumV + Vlist[[i]] }
2. 両者を乗じる
最後に、2-3で得たと
sumVを掛け合わせれば解が得られるはずです。どうでしょうか？
> cbind( + R_inv %*% sumV, + solve(crossprod(x, x)) %*% crossprod(x, y), + solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y + ) [,1] [,2] [,3] [1,] -17.579095 -17.579095 -17.579095 [2,] 3.932409 3.932409 3.932409
同じですね！
終わりに
今回はデータを小集団に分割しながら線形回帰の解を推定する、ということを紹介しました。今の時代にどうしても必要な知識かと言えばそんなこともありませんが、知っておくとと良いこともあるよ、ということで。
なおこの記事の冒頭で紹介したこちらのページでは、さらに「データソースが複数に分かれている」条件でも線形回帰の解が推定できることを示しています（例えばデータを格納しているサーバーが複数に分かれており、しかもデータのコピーが難しい状況を想定しているようです）。こちらはなかなか実用的なのではないでしょうか？
freakonometrics.hypotheses.org
おまけ
上記の工程をtidyに実行しようとすると、以下のようになるようです（こちらから）
library(tidyverse) X <- data_frame(intercept = 1, speed = cars$speed) %>% as.matrix() y <- cars$dist mats <- X %>% as_data_frame() %>% mutate( id = rep(1:5, each = 10) , y = y ) %>% ## this is where partitioning happens nest(-id) %>% mutate( X = map(data, ~ .x %>% select(-y) %>% as.matrix()), y = map(data, ~ .x %>% pull(y)) ) %>% ## We calculate QR decomposition for each partition independently mutate( Q2 = map(X, ~ .x %>% qr() %>% qr.Q()), R1 = map(X, ~ .x %>% qr() %>% qr.R()) ) df_collect <- mats$R1 %>% do.call(what = 'rbind', args = .) data.frame(dimension = c('rows', 'columns'), cbind(X %>% dim(), df_collect %>% dim())) ## Number of groups for nesting can be automatically inferred m2 <- dim(mats$R1[[1]])[2] ## The map-stage QR-decomposition Q1 = df_collect %>% qr %>% qr.Q R2 = df_collect %>% qr %>% qr.R ## For some reason this did not work with a `mutate` command... mats$Q1 = Q1 %>% as_data_frame() %>% mutate(id = ceiling(row_number() / m2)) %>% nest(-id) %>% mutate(data = map(data, ~ as.matrix(.x))) %>% pull(data) v_sum = mats %>% mutate(Q3_t = map2(.x = Q2, .y = Q1, .f = ~ t(.x %*% .y))) %>% mutate(V = map2(.x = Q3_t, .y = y, .f = ~ .x %*% .y)) %>% pull(V) %>% reduce(`+`) t(solve(R2) %*% v_sum)
-
果たして今の時代にどれほどのニーズがあるのかわかりませんが。。。↩
過小分散なカウントデータを扱いたい
背景
カウントデータをモデリングしようと思ったとき、まず思い浮かべる分布といえばポアソン分布だと思います。しかしポアソン分布は期待値と分散が等しいという性質があるため、実際のデータに当てはめようとすると、（ポアソン分布から期待されるものよりも）分散が大きい（過分散）または小さい（過小分散）という問題に直面することがあります。
そのようなとき、過分散ならともかく1、過小分散の場合にどんな分布を当てはめれば良いのか知らなかったのですが、先日某所でダブルポアソン分布というものを教えてもらったので試してみます。
doublepoissonを触ってみる
いきなりですが、ダブルポアソン分布で調べてみたところRではrmutilパッケージでdoublepoisson関数が使えるようなので早速インストールして使ってみます。
install.packages("rmutil") library(rmutil)
rmutilでは、runifやdnormなどと同じようにr/d/p/qを頭に付けたdoublepois関数が使えます。ダブルポアソン分布に従う乱数を生成してみましょう。Overdispersion parameterとしてsが指定できるのでひとまず小さめの値にしてみます。
set.seed(123) n <- 500 m <- 5 # 平均 s <- 2 # Overdispersion parameter d <- data.frame(y = rdoublepois(n, m, s))
> head(d) y 1 4 2 6 3 5 4 7 5 8 6 3
ぱっと見た感じでは通常のポアソンと変わりありません。このときの平均と分散は：
> mean(d$y) [1] 5.016 > var(d$y) [1] 2.452649
でした。うーん、分散がsと結構違ってますが、これ合っているんですかね。。。
ひとまず次に進めます。rpoisと比べ分散が小さいデータを生成できているか確認してみましょう。カウントごとの頻度を描き2、その上に（通常の）ポアソン分布の確率関数を引いてみます。
## プロット用に最大値を得る max_count <- max(d$y) ## 平均がmのポアソン分布における確率密度を得る nrm_poi <- dpois(0:max_count, m) ## データの分布を得る dens <- table(d$y) / nrow(d) ## プロット用に最大密度を得る。まるめの単位はよしなに。 max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 ## 欠損に対応するために0からmax_countまでのベクトルを用意して値を代入する dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens ## ヒストグラムではなく棒グラフで作図するので、座標を指定する xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) ## 棒グラフに確率密度関数を載せる plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "") rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue")
通常のポアソン分布で期待される密度（青線）と比較するとダブルポアソン分布では平均値周りの密度が高くなっており、全体としてバラツキが小さくなっているようです。
ではsの値を変えるとどうなるでしょうか？簡単に変更できるように関数化しておきます。
plot_double_poisson <- function(n, m, s) { d <- data.frame(y = rdoublepois(n, m, s)) max_count <- max(d$y) nrm_poi <- dpois(0:max_count, m) dens <- table(d$y) / nrow(d) max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "", main = paste0("s = ", s), cex.main = 0.7) rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue") }
sを1, 3, 5, 10とした場合のグラフを描いてみます。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 5, 1) plot_double_poisson(500, 5, 3) plot_double_poisson(500, 5, 5) plot_double_poisson(500, 5, 10)
おや？ sが小さい方がデータのバラツキが大きくなっていますね。sが10のときにはデータの半数以上が5となっているようです。またsが1のときにはほぼポアソン分布と一致しているように見えます3。この傾向はmを変更しても変わらないようです。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 10, 1) plot_double_poisson(500, 10, 3) plot_double_poisson(500, 10, 5) plot_double_poisson(500, 10, 10)
想像していた結果とは異なるのですが、このOverdispersion parameterを変更することでポアソン分布では扱えなかった過小分散なカウントデータに対応できそうです。
切片を推定してみる
glmで推定してみる
それでは続いて、この過小分散なデータに対してGLMを当てはめたときの結果を確認してみましょう。改めてサンプルデータを作成しますが、十分に過小分散となるようsは10としておきましょう。
set.seed(123) n <- 500 m <- 5 s <- 10 d <- data.frame(y = rdoublepois(n, m, s))
glmを使ってフィッティングしますが、ここで切片として5が推定されれば問題ありません。またリンク関数はlogとします。
res_glm <- glm(y ~ 1, d, family = poisson("log"))
> coef(res_glm) (Intercept) 1.609838
切片は1.609と表示されていますが、逆リンク関数expで戻してみると、
> exp(coef(res_glm)) (Intercept) 5.002 > mean(d$y) [1] 5.002
一致しています。過小分散であっても切片は正しく推定されるようです。しかしポアソン分布であると仮定した場合、その期待値と分散は等しい（E(Y) = Var(Y)）はずなのですが、実際の分散を見てみると：
> var(d$y) [1] 0.4709379
大きく異なります。
dglmで推定してみる
それではこのようなデータに対してどのようにモデリングすれば良いかということなのですが、調べてみたところdglmというパッケージ・関数があるようなので使ってみましょう。
install.packages("dglm") library(dglm)
> res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log")) stats::glm.fit(x = Z, y = d, weights = rep(1/2, N), offset = doffset, でエラー: 'x' の中に NA/NaN/Inf があります
あらら、何らかのエラーで止まってしまいました。原因を探りたいところですが、ひとまずsを小さくして先に進みます。
set.seed(123) n <- 500 m <- 5 s <- 3 d <- data.frame(y = rdoublepois(n, m, s))
res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log"))
今度は大丈夫なようですので結果を見てみましょう。
> exp(res_dglm$coefficients[1]) (Intercept) 5.01 > mean(d$y) [1] 5.01
今度もmは推定できているようです。分散はどうでしょうか。
> exp(res_dglm$dispersion.fit$coefficients) (Intercept) 0.3372774 > var(d$y) [1] 1.685271
なんか、やっぱり分散は違いますね。。。
実はここまで触れてこなかったのですが、ダブルポアソン分布についてはこちらのブログを参考にさせて頂きました：
この記事ではダブルポアソン分布の分散として
と紹介されていますので、もしかしたら分散はmとの積を求める必要があるかもしれません。
> sqrt(exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients)) (Intercept) 1.299908 > var(d$y) [1] 1.685271
さっきよりもだいぶ近い値が得られるようになりましたが、右辺が平方根になっているのが気になります（分散なので）。取ってみましょう。
> exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients) (Intercept) 1.68976 > var(d$y) [1] 1.685271
だいぶ近い！…のですが、分散の定義…。
回帰係数を推定してみる
glm、dglmで推定してみる
なかなか良さそうな感触が得られたので、切片だけでなく回帰係数も推定してみましょう。以下のようにサンプルデータを作成します。sは3としました。
set.seed(123) n <- 500 b <- 1.5 x <- runif(n, 3, 5) m <- b * x s <- 3 y <- rdoublepois(n, m, s) d <- data.frame(y, x)
平均mを与えたときのダブルポアソン分布を見てみましょう。
hist(d$y, col = "gray", ylim = c(0, 0.5), probability = T, main = "", xlab = "", breaks = "scott")
また平均と分散を確認すると：
> mean(d$y) [1] 5.982 > var(d$y) [1] 2.751178
このようになりました。
それぞれのデータを、glmとdglmそれぞれでフィッティングしてみます。まずはglmから。
res_glm <- glm(y ~ x, d, family = poisson("log"))
> exp(coef(res_glm)) (Intercept) x 2.057136 1.302974
回帰係数は真の値(1.5)にまぁまぁ近いものが推定されています。続いてdglmでフィッティングすると：
res_dglm <- dglm(y ~ x, ~1, data = d, family = poisson("log"))
> exp(res_dglm$coefficients) (Intercept) x 2.057136 1.302974
glmと同じ結果が得られました。
では分散を見てみましょう。mにはd$xの平均値を用いることにします。
m_hat <- exp(res_dglm$coefficients %*% c(1, mean(d$x)))
> m_hat [,1] [1,] 5.914556 > mean(m) [1] 5.985851
> exp(res_dglm$dispersion.fit$coefficients) * m_hat [,1] [1,] 1.947131
うーん、合いません（真の値は3）。。。
ちなみにこのデータをlmでフィッティングするとどうなるかと言うと：
res_lm <- lm(y ~ x, d)
> var(res_lm$residuals) [1] 1.935673
あまり変わりませんね。
optimで推定してみる
もう少しチャレンジしてみましょう。dglmの代わりにoptimで推定してみます。目的関数を以下のように定義します：
my_obj_fun <- function(b) { b1 <- b[1] b2 <- b[2] b3 <- b[3] ret <- sum(ddoublepois(d$y, exp(b1 + d$x * b2), b3, log = TRUE)) ret }
パラメータのベクトルを適当に用意して、optimにかけてみましょう。
b <- c(0, 1, 3) res_opt <- optim(b, my_obj_fun, control = list(fnscale = -1))
結果を確認します：
> exp(res_opt$par)[-3] # 回帰係数 [1] 2.050437 1.303756 > res_opt$par[3] # Overdispersion [1] 3.064928
おぉ！回帰係数はglm、dglmの結果と大体一致しています。またsの推定値は3.06となっており、これも真の値と近しいものになっています。optimいいですね。
終わりに
今回は過小分散なカウントデータを扱うためにrmutilとdglmというパッケージを使ってみました。optimによる推定は思いの外うまくいったようですが、dglmの結果は少しパッとしませんでした。ただしsの求め方については本当にこれで合っているのか自信がないため、推定できているかについての判断は保留です。ご存知の方がいらっしゃれば教えてください。
FindBestSplitを書いてみる
背景
前回、前々回の記事でrandomForestを使ってみたのですが、ソースコードを読んでいるとノードの分割においてfindbestsplitというサブルーチンが使われていることに気が付きました。このサブルーチン自体はこちらのL191に定義されているのでそれを読めばわかる（はずな）のですが、もう少しわかりやすい説明はないかなーと探してみたところ、こんな解説記事を見つけました。
http://dni-institute.in/blogs/cart-algorithm-for-decision-tree/
これによると、どうやらfindbestsplitは
というステップによって最良の閾値を探しているようです。それほど難しくなさそうなので、これをRで書いてみましょう。
実装
findbestsplitを実装するためには、以下のような関数が必要となりそうです。
- データ、説明変数を与えると閾値の候補を返す（return_threshold_values）
- データ、目的変数、説明変数、閾値を与えるとGini係数を返す（return_gini_index）
- 現在のGini係数との差分が最大となる（最良な）閾値を返す（return_best_value）
まずは１つ目から書いてみましょう。
1. データ、説明変数を与えると閾値の候補を返す関数
先ほど紹介したページでは閾値の候補を生成するための方法として
One of the common approach is to find splits /cut off point is to take middle values
との説明がありましたので、これに倣います。以下のように書いてみました：
- 説明変数列を
uniqueする sortで並び替える- 各要素について1つ前の値との差分（
diff）を取る - 差分を2で割る
sort後の列（最初の要素は除く）に加える
return_threshold_values <- function(dat, col) { uniq_vals <- sort(unique(dat[, col])) diffs <- diff(uniq_vals) / 2 thre_vals <- uniq_vals[-1] - diffs return(thre_vals) }
2. データ、目的変数、説明変数、閾値を与えるとGini係数を返す関数
続いてGini係数を求める関数を定義します。ここでは引数として閾値も与え、後ほどapplyで閾値候補をまとめて並列に処理することを考えました。Gini係数の求め方はここを参考に、以下のように書きました：
- 説明変数を所与の閾値で1/0のカテゴリに振り分ける
- 全体および各カテゴリのサンプルサイズを求める
- 目的変数 × 説明変数による混同行列の各要素の割合（の二乗）を計算する
- Gini係数を求める
return_gini_index <- function(dat, target, col, val) { d <- dat[, c(target, col)] d$cat <- ifelse(d[, col] > val, 1, 0) n0 <- nrow(d) n1 <- sum(d$cat) n2 <- n0 - n1 p11 <- (nrow(d[d$cat == 1 & d$target == 1, ]) / n1)^2 p12 <- (nrow(d[d$cat == 1 & d$target == 0, ]) / n1)^2 p21 <- (nrow(d[d$cat == 0 & d$target == 1, ]) / n2)^2 p22 <- (nrow(d[d$cat == 0 & d$target == 0, ]) / n2)^2 gini_val <- (n1/n0) * (1 - (p11 + p12)) + (n2/n0) * (1 - (p21 + p22)) return(gini_val) }
3. 現在のGini係数との差分が最大となる（最良な）閾値を返す関数
上記の処理によってある閾値におけるGini係数を求めることが出来ましたので、これをapplyで並列化します。また実際のところ必要な値はGini係数そのものではなく、現時点におけるGini係数との差分なので、それも計算しましょう。これまでに定義した関数を使って以下のように書きます：
- 閾値候補のベクトルを
listにする - 閾値候補リストを
return_gini_indexにsapplyで渡す - 現時点のGini係数を計算する
- 差分の最大値、Gini係数の最小値を得る
- 差分が最大となる閾値を得る
- 変数名と合わせて返す
差分の最大値やGini係数の最小値は別に必要ないのですが、参考のために取っておきます。
return_best_val <- function(dat, target, col) { thre_vals <- as.list(return_threshold_values(dat, col)) gini_vals <- sapply(thre_vals, return_gini_index, dat = dat, target = target, col = col) p1 <- sum(dat$target) / nrow(dat) p2 <- 1 - p1 current_gini <- 1 - ((p1)^2 + (p2)^2) max_gap <- max(current_gini - gini_vals) min_gini <- min(gini_vals) max_thre_val <- thre_vals[[which(max_gap == current_gini - gini_vals)]] return(c(col, max_thre_val, round(max_gap, 2), round(min_gini, 2))) }
では、このようにして定義した関数を実際に当てはめまめてみます。irisを使って適当に以下のようなデータを用意します。
my_iris <- iris my_iris$target <- ifelse(my_iris$Species == "setosa", 1, 0) dat <- my_iris[, -5]
まずはSepal.Lengthの最良な閾値を取得してみます。
> return_best_val(dat, "target", "Sepal.Length") [1] "Sepal.Length" "5.45" "0.3" "0.14"
これが合っているのかはわかりませんが、続いて全部の変数に同時に当てはめてみましょう。lapplyとdo.callを使います。
> do.call("rbind", + lapply(as.list(colnames(dat)[1:4]), return_best_val, dat = dat, target = "target")) [,1] [,2] [,3] [,4] [1,] "Sepal.Length" "5.45" "0.3" "0.14" [2,] "Sepal.Width" "3.35" "0.17" "0.28" [3,] "Petal.Length" "2.45" "0.44" "0" [4,] "Petal.Width" "0.8" "0.44" "0"
各変数について最良な閾値を取得することが出来たようです。ところで4列目（Gini係数の最小値）を確認すると、Petal.LengthとPetal.Widthで0になっていますが、Gini係数が0ということは完全に分離されていることを意味します。確かめてみましょう。
> plot(dat$Petal.Length, col = dat$target + 1) > abline(h = 2.45) > plot(dat$Petal.Width, col = dat$target + 1) > abline(h = 0.8)
確かに完全に分離しているようです。
終わりに
今回はrandomForestの中でも使われているfindbestsplitというアルゴリズムをRで書いてみました。実際には決定木やRandom Forestは一連の処理を再帰的に繰り返しているのですが、最も重要なポイントはこちらになるのだと思います。なおPythonのDecisionTreeClassifierでも同じようなアルゴリズムとなっているのか確認したかったのですが、Pythonのソースコードの表示方法が良くわかりませんでした。
おしまい。
randomForestで有効な交互作用を発見したい
背景
GLMは発想がわかりやすく解釈性も高くて良いアルゴリズム1なのですが、線形の仮定があるため変数間の交互作用を見るのが苦手です。実際のプロジェクトでGLMを使った結果を見せ、
- 変数の組み合わせ効果みたいなものは見れないの？
- この変数は条件によって効き方が違うんだよね〜
みたいな指摘を受けて困った経験があったりしないでしょうか。そんな時に使えるテクニックを同僚から教えてもらったので、備忘がてらメモしておきます。勝手に公開して怒られる可能性もありますが。。。
概要
手法の概要ですが、話としてはシンプルで「もしも有効な変数の組み合わせ（交互作用）が存在しているのであれば、Random Forestの各決定木において、ノードの分岐に使われる変数の順番として出現しやすいペアがあるのではないか」ということです。例えば変数X1とX2の間に交互作用があれば、決定木においてX1が選択された場合、続く分岐ではX2が選択されやすくなるのではないでしょうか。
実装
上記のアイディアを実現するために、以下のように実装してみます：
- Random Forestでモデルを作る
- 各決定木から分岐に用いられた変数ペアを得る
- 出現回数のカウントを取る
- 交互作用効果を確かめてみる
1. Random Forestでモデルを作る
まずはRandom Forestでモデルを作ります。randomForestパッケージを使ってサクッと作りましょう。
### libraryの読み込み library(randomForest) library(tidyverse)
データには前回記事と同じTelco Customer Churnを使いますが、前回の反省を踏まえてread.csvを使います。
前回の記事はこちら。
d <- read.csv("./Data/WA_Fn-UseC_-Telco-Customer-Churn.csv") %>% as_data_frame()
本来ならここから一つ一つの変数を観察するところですが、今回はそれが目的ではないので欠損だけ埋めておきます。
> colSums(apply(d, c(1, 2), is.na)) customerID gender SeniorCitizen Partner Dependents tenure 0 0 0 0 0 0 PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection 0 0 0 0 0 0 TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod 0 0 0 0 0 0 MonthlyCharges TotalCharges Churn 0 11 0
TotalChargesに欠損があるようですね。
> summary(d$TotalCharges) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 18.8 401.4 1397.5 2283.3 3794.7 8684.8 11
MedianとMeanに差があるので分布が偏っていそうです。ひとまずNAはMedianで埋めておきましょう。
d2 <- d %>% mutate(TotalCharges = if_else(is.na(.$TotalCharges), median(.$TotalCharges, na.rm = T), .$TotalCharges))
customerIDは変数として使えないので除外しましょう。またrandomForestはカテゴリ数が53より多い変数を扱えないので、カテゴリ数をチェックしておきます。
cat_vars <- sapply(d2, is.factor)
> apply(d2[, cat_vars], 2, function(x) length(unique(x))) customerID gender Partner Dependents PhoneService MultipleLines 7043 2 2 2 2 3 InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV 3 3 3 3 3 3 StreamingMovies Contract PaperlessBilling PaymentMethod Churn 3 3 2 4 2
大丈夫そうですね。customerIDだけ落としておきます。
d3 <- d2 %>% select(-customerID)
randomForestを当てはめます。目的変数はChurnです。
set.seed(123) result <- randomForest(Churn ~ ., d3, ntree = 500)
> result Call: randomForest(formula = Churn ~ ., data = d3, ntree = 500) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 4 OOB estimate of error rate: 20.3% Confusion matrix: No Yes class.error No 4643 531 0.1026285 Yes 899 970 0.4810059
精度とかは気にしません。
2. 各決定木から分岐に用いられた変数ペアを得る
次にRandom Forestから変数ペアを取得します。そのためにはRandom Forestの各決定木について、どの変数がどの順番で分岐に用いられたかを知る必要があります。
まずはRandom Forestから各決定木の結果を取ってきましょう。 getTree 関数を使います。過去記事も参考にしてください。
> getTree(result, 1, labelVar = TRUE) left daughter right daughter split var split point status prediction 1 2 3 Contract 1.000 1 <NA> 2 4 5 OnlineBackup 1.000 1 <NA> 3 6 7 TechSupport 2.000 1 <NA> 4 8 9 tenure 2.500 1 <NA> 5 10 11 InternetService 2.000 1 <NA> 6 12 13 Dependents 1.000 1 <NA> # 中略
getTree は分岐に用いられた変数とその分岐先ノードなどを返します。ここで必要なのは1~3列目なのですが、それぞれ「左の子ノード」「右の子ノード」「分岐に用いられた変数」を意味しています。例えば1行目を見るとノードの分岐に Contract が用いられたことがわかります。
このテーブルに行インデックスを列として追加しましょう。なおこれ以降のコードはこちらの記事を参考にさせて頂きました。
tree_tbl <- getTree(result, 1, labelVar = TRUE) %>% # labelVar = Fだとエラー rownames_to_column() %>% mutate(rowname = as.integer(rowname))
作成した tree_tbl では分岐に用いた変数（ split var ）はわかりますが、変数のペアはわかりません。例えば1行目の Contract で分岐された子ノード（2と3）は、次にそれぞれが OnlineBackup と TechSupport で分岐されています（2行目、3行目）ので、 Contract - OnlineBackup と Contract - TechSupport という変数ペアが出現したことがわかるような形に整形したいですね。
各行には「分岐に用いた変数」と「分岐先の子ノードの番号」がありますので、「分岐先の子ノード（左右両方）」に「分岐元の変数」を追加すれば欲しいものが得られそうです。
まずはノードと変数のマスタを用意しましょう。
var_name <- tree_tbl %>% select(rowname, "split var") %>% rename(split_var =`split var`) %>% # スペースを`_`に修正 unique() %>% filter(!is.na(.$split_var))
続けて left daughter と right daughter それぞれに rowname でJOINします。
> tree_tbl %>% + left_join(var_name, by = c("left daughter" = "rowname")) %>% + left_join(var_name, by = c("right daughter" = "rowname")) %>% + select(rowname, `split var`, `split_var.x`, `split_var.y`) %>% + na.omit() rowname split var split_var.x split_var.y 1 1 Contract OnlineBackup TechSupport 2 2 OnlineBackup tenure InternetService 3 3 TechSupport Dependents StreamingMovies 4 4 tenure InternetService InternetService 5 5 InternetService Dependents OnlineSecurity 6 6 Dependents TotalCharges tenure # 中略
良さそうですね。ただしこのままでは後の工程で使いにくいのでもう少し加工します。
> tree_tbl %>% + left_join(var_name, by = c("left daughter" = "rowname")) %>% + left_join(var_name, by = c("right daughter" = "rowname")) %>% + select(`split var`, `split_var.x`, `split_var.y`) %>% + na.omit() %>% + rename(from_var = `split var`, + left = `split_var.x`, + right = `split_var.y`) %>% + gather(key = node, value = to_var, -from_var) %>% + select(-node) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges # 中略
これで必要なアウトプットが得られました。あとは上記の加工を関数化しておき、各決定木に当てはめれば良さそうです。
get_var_pairs <- function(tree_num, rf = result) { # 決定木の結果を得る tree_tbl <- getTree(rf, tree_num, labelVar = TRUE) %>% rownames_to_column() %>% mutate(rowname = as.integer(rowname)) var_name <- tree_tbl %>% select(rowname, "split var") %>% rename(split_var =`split var`) %>% # スペースを`_`に修正 unique() %>% filter(!is.na(.$split_var)) out <- tree_tbl %>% left_join(var_name, by = c("left daughter" = "rowname")) %>% left_join(var_name, by = c("right daughter" = "rowname")) %>% select(`split var`, `split_var.x`, `split_var.y`) %>% na.omit() %>% rename(from_var = `split var`, left = `split_var.x`, right = `split_var.y`) %>% gather(key = node, value = to_var, -from_var) %>% select(-node) return(out) }
試してみましょう。
> get_var_pairs(5, result) from_var to_var 1 PaymentMethod Contract 2 Contract TechSupport 3 OnlineSecurity tenure 4 TechSupport PaperlessBilling 5 tenure TechSupport 6 TotalCharges StreamingTV # 中略
これを全ての決定木に当てはめます。 purrr の map_dfr を使ってみましょう。
var_pairs <- map_dfr(as.matrix(1:5), get_var_pairs, result) # 少しだけ実行
> dim(var_pairs) [1] 2950 2 > head(var_pairs) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges
素直に for で書いた時と同じ結果になっていますでしょうか？
tmp <- c() for (i in 1:5) { tmp <- bind_rows(tmp, get_var_pairs(i, result)) }
> dim(tmp) [1] 2950 2 > head(tmp) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges
合っているようなので全ての結果を取得します。 randomForest では作成する決定木の数を ntree で指定するので 1:ntree で全ての決定木に適用できるのですが、直接 ntree を取ってくることは出来ないようなので length(treesize()) を使います。エラーは気にしないことにします。
var_pairs <- map_dfr(as.matrix(1:length(treesize(result))), get_var_pairs, result)
30万弱の変数ペアが得られました。
3. 出現回数のカウントを取る
さっそく変数ペアのカウントを取ってみましょう。
> var_pairs %>% + group_by(from_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 19 x 2 from_var cnt <chr> <int> 1 TotalCharges 36908 2 MonthlyCharges 36592 3 tenure 34132 4 PaymentMethod 21244 5 gender 17626 6 MultipleLines 14496 # 中略
分岐の始点となった変数（分岐元）の数は19ですが、分析対象のデータセットには目的変数を含んで20列だったので、分岐元とならない変数はなかったようです。一方で分岐元としての出現頻度には大きなばらつきがあり、 TotalCharges 、 MonthlyCharges 、 tenure が選ばれやすいようですね。
ちなみに varImpPlot で変数重要度を見てみると、これらはいずれも上位に付けており、4位以下と大きな隔たりがあるようです。
varImpPlot(result)
続いて分岐の終点となった変数（分岐先）についても見てみましょう。
> var_pairs %>% + group_by(to_var) %>% + # group_by(from_var, to_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 19 x 2 to_var cnt <chr> <int> 1 TotalCharges 47553 2 MonthlyCharges 47137 3 tenure 37964 4 PaymentMethod 20138 5 gender 14242 6 Partner 12567 # 中略
同じく19変数ありますので、全ての変数は分岐元・分岐先ともに出現しています。分岐元と同じく出現頻度はばらつきがあり、出現しやすい変数としては、 TotalCharges 、 MonthlyCharges 、 tenure となっています。これは少し意外ですね。てっきり分岐元に選ばれやすい変数と分岐先に選ばれやすい変数は違うものになると思っていましたが。
せっかくなので分岐元と分岐先で選ばれやすさが異なるか、可視化してみましょう。
plt <- var_pairs %>% group_by(from_var) %>% summarise(cnt_f = n()) %>% left_join(var_pairs %>% group_by(to_var) %>% summarise(cnt_t = n()), by = c("from_var" = "to_var")) %>% gather(var_type, cnt, -from_var) %>% rename(var = from_var) ggplot(plt, aes(x = reorder(var, -cnt), y = cnt, fill = var_type)) + geom_bar(stat = "identity", position = "dodge") + # scale_color_brewer(palette = "Set2") + # facet_wrap(~var_type, nrow = 2) + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + NULL
上位3変数は分岐元として選ばれやすいですが、分岐先としては更に頻度が多くなっています。また変数間の順位にはほとんど変動はないようですね。
組み合わせでも見てみましょう。
> var_pairs %>% + group_by(from_var, to_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 355 x 3 # Groups: from_var [19] from_var to_var cnt <chr> <chr> <int> 1 MonthlyCharges TotalCharges 6032 2 TotalCharges MonthlyCharges 6018 3 TotalCharges TotalCharges 5858 4 MonthlyCharges MonthlyCharges 5623 5 tenure MonthlyCharges 5602 6 tenure TotalCharges 5407 # 中略
上位10件の組み合わせを見ると、9行目の PaymentMethod を除いていずれも上位3変数の組み合わせになっています。同じ変数のペアも出てきているので、レンジを絞る形で分岐条件として選ばれているようですね、なるほど。なお19 * 19 = 361なので発生していない組み合わせがあるようですが、一部ですね。
4. 交互作用効果を確かめてみる
ひとまず目的としていた分析は以上となります。今回のデータセットおよび分析条件を用いた場合、 TotalCharges 、 MonthlyCharges 、 tenure の3変数が（組み合わせの意味でも、重要度の意味でも）影響の大きい変数であるようです。したがって冒頭のようなクライアントからの指摘があった場合には、特にこの3変数を中心に他の変数との交互作用を確認していくと良いのではないでしょうか。
と、ここまで書いたところで一つ疑問が浮かんできました。このような組み合わせ効果は、GLMでも発見できないでしょうか？例えば二次の交互作用項を準備しておき、Lassoで変数選択させるとこれらの組み合わせが残らないでしょうか？
やってみましょう。yとxを用意します。
library(glmnet) # glmnetを使う y <- as.matrix(ifelse(d3$Churn == "Yes", 1, 0)) tmp <- scale(model.matrix(Churn ~ .^2 , d3)) x_vars <- which(colSums(apply(tmp, c(1,2), is.nan)) == 0) x <- cbind(1, tmp[, x_vars])
> dim(x) [1] 7043 403
Lassoは回帰係数の絶対値に対して罰則が与えられるため、説明変数のスケールを揃えておく必要があります。そのため model.matrix でダミー化したあと scale で正規化しています。またその際に分散が0であるために NaN となってしまう変数は除外しています。
このデータでLassoにかけてみましょう。まずは適切な lambda を得るために cv.glmnet を使いますが、計算時間が少しかかるため nfolds は5にしておきましょう。
res_lasso_cv <- cv.glmnet(x = x, y = y, family = "binomial", alpha = 1, nfolds = 5)
> res_lasso_cv$lambda.min [1] 0.004130735
このときにDevianceが最小となる lambda は0.004130735のようです。プロットも見ておきましょう。
plot(res_lasso_cv)
ここで縦に引かれた二本の破線はそれぞれ lambda.min および lambda.1se を表しています。 lambda.1se は lambda.min の1SD以内で最も罰則を与えたときの lambda を示すようですね。詳しくは以下を参照してください。
https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
ひとまず lambda.min を与えたときの結果を確認しましょう。
res_lasso <- glmnet(x = x, y = y, family = "binomial", alpha = 1, lambda = res_lasso_cv$lambda.min)
# このときの回帰係数の絶対値 > as.data.frame(as.matrix(res_lasso$beta)) %>% + rownames_to_column() %>% + filter(s0 != 0) %>% + mutate(abs_beta = abs(s0)) %>% + arrange(desc(abs_beta)) %>% + select(rowname, abs_beta, s0) rowname abs_beta s0 1 tenure 7.708817e-01 -7.708817e-01 2 ContractTwo year 4.658562e-01 -4.658562e-01 3 ContractOne year 2.723315e-01 -2.723315e-01 4 InternetServiceFiber optic 2.640571e-01 2.640571e-01 5 InternetServiceNo 2.202043e-01 -2.202043e-01 6 tenure:PaymentMethodMailed check 1.714972e-01 -1.714972e-01 # 中略
むむ。。Lassoにおいても tenure は影響の大きい（回帰係数の絶対値が大きい）変数として選ばれましたが、 TotalCharges と MonthlyCharges はいませんね。
っていうか、
> as.data.frame(as.matrix(res_lasso$beta)) %>% + rownames_to_column() %>% + filter(rowname %in% c("TotalCharges", "MonthlyCharges")) rowname s0 1 MonthlyCharges 0 2 TotalCharges 0
Lassoで落とされてる。。。
TotalCharges がどのような影響を示すのか partialPlot で見てみましょう。
partialPlot(result, as.data.frame(d3), "TotalCharges") # tibbleのままだとエラーになる
グニャグニャですね。特定のレンジでは影響が大きいものの、他ではそうでもないということなんでしょうか。だからRandom Forestのような非線形のアルゴリズムだと効果が認められる一方、Lassoのような線形のアルゴリズムでは拾いきれないのかもしれません2。これは素直に、Random Forestの結果から効果のありそうな組み合わせ変数を見つけ、分布を見ながら組み込んだ方が良さそうです。
しかしy軸が1を超えるのはなぜなんでしょうか。。。
終わりに
今回の分析はRandom Forestの結果から交互作用の良い候補を見つけようという趣旨でした。また同様の結果がLassoからも得られるかを検証しましたが、両者の結果は異なるものとなりました。Random Forestは非線形な効果を捉えることができるアルゴリズムなのでこちらの結果から有効な変数ペアを絞り込み、一つずつ検証していくスタイルが良さそうです。
randomForestではCharacterは使わないようにしよう
RのrandomForestを使っていてはまったのでメモしておきます。
①目的変数がcharacterだと分類として扱ってくれない
最初にはまったのがこちらでした。目的変数がcharacterだとカテゴリ変数として扱ってもらえないため、分類ではなく回帰としてプログラムが進んでしまい、エラーが返ります。
まずは以下のようにデータを読み込みます。ちなみにこのデータはTelco Customer Churnで、
KaggleからDLしてきました。
library(tidyverse) library(randomForest) d <- read_csv("./Data/WA_Fn-UseC_-Telco-Customer-Churn.csv")
このデータを使って下記のようにrandomForestを実行します：
> randomForest(Churn ~ gender, d) y - ymean でエラー: 二項演算子の引数が数値ではありません
「何このエラー？？」と思いながらrandomForestの中身を見てみると、以下のような記述が見つかります。
addclass <- is.null(y) classRF <- addclass || is.factor(y)
classRFというのはカテゴリ変数であるかを判定するbooleanなのですが、is.factorで判断しているんですね。ここでFALSEが返ると、後の工程で分類のためのプログラム（Cで書かれたもの、多分コレのL38 ~ L540）ではなく、回帰用のプログラム（多分コレのL22 ~ L340）が呼ばれてしまうようです。
ちなみにrandomForestはformulaでもmatrixでも対応できるような総称関数になっているので、
コンソールでrandomForestと叩いても中身を見ることはできません。そのような場合、まずはmethodsでどのような関数が含まれているかを確認しましょう。
> methods(randomForest) [1] randomForest.default* randomForest.formula* see '?methods' for accessing help and source code
randomForestという関数はrandomForest.defaultとrandomForest.formulaという関数を総称しているようです。前者がメインのようなので以下のコマンドを実行します。
getS3method("randomForest", "default")
そうすると先程の関数定義を確認することができます。さらにちなみに、lookupというパッケージを用いることでそういった総称関数やCで書かれた関数などをRStudio上でシンタックスハイライトさせながら表示することができます。大変便利ですので、是非こちらの記事を参考にして使ってみてください。
②説明変数がcharacterだとダミー化してくれない
次にはまったのがこちらでした。さきほどエラーが返ってきたrandomForestを、目的変数をfactorに直して実行してみましょう：
d2 <- d %>% mutate(Churn = as.factor(Churn))
> randomForest(Churn ~ gender, d2) 強制変換により NA が生成されました randomForest.default(m, y, ...) でエラー: 外部関数の呼び出し (引数 1) 中に NA/NaN/Inf があります
NAがありますというエラーなのですが、このデータではNAはTotalCharges列にしかありませんので、どうやら途中で生成されているようです。ちなみにRStudioでRMarkdownを使っているとき、チャンク内で実行すると上記のエラーが表示されるのですが、コンソールで実行すると以下のようにもう少し情報が追加されます：
> randomForest(Churn ~ gender, d2) randomForest.default(m, y, ...) でエラー: 外部関数の呼び出し (引数 1) 中に NA/NaN/Inf があります 追加情報: 警告メッセージ: data.matrix(x) で: 強制変換により NA が生成されました
data.matrixが悪さしているようですね。私は普段チャンク内で実行させることが多いため、
この表示に気付かなくて時間を無駄にしました。
randomForest.defaultを見ると以下の記述があります。
if (is.data.frame(x)) { xlevels <- lapply(x, mylevels) ncat <- sapply(xlevels, length) ncat <- ifelse(sapply(x, is.ordered), 1, ncat) x <- data.matrix(x)
ここでxをdata.matrix(x)でダミー化しようとしたものの、characterであったために失敗しているようですね。
> d2 %>% select(Churn, gender) %>% str() Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 7043 obs. of 2 variables: $ Churn : Factor w/ 2 levels "No","Yes": 1 1 2 1 2 2 1 1 2 1 ... $ gender: chr "Female" "Male" "Male" "Male" ...
> d2 %>% select(Churn, gender) %>% data.matrix(.) %>% head() 強制変換により NA が生成されました Churn gender [1,] 1 NA [2,] 1 NA [3,] 2 NA [4,] 1 NA [5,] 2 NA [6,] 2 NA
それではfactorに直して実行してみましょう。
d3 <- d2 %>% select(Churn, gender) %>% mutate(gender = as.factor(.$gender))
> str(d3) Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 7043 obs. of 2 variables: $ Churn : Factor w/ 2 levels "No","Yes": 1 1 2 1 2 2 1 1 2 1 ... $ gender: Factor w/ 2 levels "Female","Male": 1 2 2 2 1 1 2 1 1 2 ...
> randomForest(Churn ~ gender, d3) Call: randomForest(formula = Churn ~ gender, data = d3) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 1 OOB estimate of error rate: 26.54% Confusion matrix: No Yes class.error No 5174 0 0 Yes 1869 0 1
ようやく動くようになりました。
このようなエラーにはまらないためには、例えばread_csvの代わりに（StringsAsFactorsをFALSEにしないで）read.csvを使うか、randomForestの代わりにrangerを使うという手があります。
> ranger::ranger(Churn ~ gender, d) Ranger result Call: ranger::ranger(Churn ~ gender, d) Type: Classification Number of trees: 500 Sample size: 7043 Number of independent variables: 1 Mtry: 1 Target node size: 1 Variable importance mode: none Splitrule: gini OOB prediction error: 26.54 %
もとのデータでも動いていますね！
終わりに
今回のエラーはモデリングに入る前にダミー化を行っていれば防げたものでしたが、Rだとカテゴリ変数をそのまま渡しても動くものが多いのでついサボってしまいました。普段からダミー化をする習慣が身についている人や、多分Pythonでははまらないんでしょうね。
ただ今回はエラーを追いかけながら総称関数のソースの便利な確認方法を知ることができたので良かったです。
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 を使うことができそうですね。
-
一部修正あり↩
