Bayesian GLMで過適合を避ける
長らく知人に貸していた「みんなのR」が手元に戻ってきたのでパラパラと見ていたら、見逃していた面白いトピックを見つけたので紹介します。内容としてはこちらの記事とほぼ同じで、基本的には書籍の該当部分を再現しているだけですが、{purrr}のmapを使った書き方にしてみたり、少しだけ変更を加えています。
書籍はこちら。P320の「Bayesian 縮小」からです。
- 作者: Jared P. Lander,Tokyo.R（協力）,高柳慎一,牧山幸史,簑田高志
- 出版社/メーカー: マイナビ
- 発売日: 2015/06/30
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (7件) を見る
分析環境の準備
ライブラリの読み込み
まずは必要なライブラリの準備から。以下のライブラリを使用します。
library(MASS) library(coefplot) library(tidyverse)
今回使うBayesian GLMのための関数bayesglmは{arm}というライブラリにあるのですが、{coefplot}とnamespaceが衝突するためarm::bayesglmとして実行することにします。
データのロード
続いて「みんなのR」の著者が公開しているサンプルデータをこちらのリンク先から取得し、適当な場所に保存します（urlを使ってloadしたのですが、うまくいきませんでした）。
保存したデータを以下のようにloadします。
load("適当なフォルダ/ideo.rdata") # ideoというデータフレームが読み込まれる
フィッティング
まずは、（ベイズではない）通常のモデリングをした場合の問題点を明らかにするためにglmでフィッティングしてみます。ここでideo$Yearごとに同じモデルをあてはめるのですが、折角なのでサブグループごとにフィッティングを行うための方法をいくつか試してみましょう。
1. glmのsubsetで分割する
まずはそのまま、glmのオプションsubsetでYearを指定してデータを分割する方法です。
theYears <- unique(ideo$Year) results_01 <- vector(mode = "list", length = length(theYears)) names(results_01) <- theYears for (i in theYears) { results_01[[as.character(i)]] <- glm(Vote ~ Race + Income + Gender + Education, data = ideo, subset = Year == i, family = binomial(link = "logit")) }
2. purrr::mapを用いる
続いて、同じことを{purrr}ライブラリのmap関数を使って実行します。
results_02 <- ideo %>% split(.$Year) %>% map(~glm(Vote ~ Race + Income + Gender + Education, data = ., family = binomial(link = "logit")))
mapを使えばかなりスッキリ書けますね。なおsplitは{base}の関数です。
3. tidyr::nestと組み合わせる
さらに続いて、{tidyr}のnestと組合わせてみましょう。まずはmapに渡す関数を定義しておきます。
追記：この書き方は書籍「Rではじめるデータサイエンス」の20章「purrrとbroomによる多数のモデル」を参考にしています。
year_model <- function(df) { glm(Vote ~ Race + Income + Gender + Education, data = df, family = binomial(link = "logit")) }
次にideo$Yearでgroup_byし、nest()で入れ子にします。
year_data <- ideo %>% group_by(Year) %>% nest()
nestすると新しいデータフレームのdata列にYearごとのデータフレームが格納されることになり、慣れないとかなり戸惑いがあります。
> year_data # A tibble: 14 x 2 Year data <dbl> <list> 1 1948 <tibble [374 x 7]> 2 1952 <tibble [1,217 x 7]> 3 1956 <tibble [1,245 x 7]> 4 1960 <tibble [882 x 7]> 5 1964 <tibble [1,100 x 7]> 6 1968 <tibble [869 x 7]> 7 1972 <tibble [1,560 x 7]> 8 1976 <tibble [1,292 x 7]> 9 1980 <tibble [829 x 7]> 10 1984 <tibble [1,343 x 7]> 11 1988 <tibble [1,172 x 7]> 12 1992 <tibble [1,304 x 7]> 13 1996 <tibble [1,004 x 7]> 14 2000 <tibble [1,049 x 7]>
このデータに対して先ほど定義した関数を適用します。
results_03 <- map(year_data$data, year_model) names(results_03) <- as.character(unique(ideo$Year))
この方法では結果のlistにYearがnamesとして渡らないため後からつけているのですが、データと名前があっているのか不安になってしまいますね。
結果の確認
では結果を見てみましょう。multiplotはlistの各要素にlmやglmの結果を持つオブジェクトを渡し、関心のある変数の回帰係数を比較することができます。
いずれも同じ結果となっています。
multiplot(results_01, coefficients = "Raceblack", secret.weapon = TRUE) + coord_flip(xlim = c(-20, 10)) multiplot(results_02, coefficients = "Raceblack", secret.weapon = TRUE) + coord_flip(xlim = c(-20, 10)) multiplot(results_03, coefficients = "Raceblack", secret.weapon = TRUE) + coord_flip(xlim = c(-20, 10))
（同じグラフなので二枚は省略）
さてここで結果を見てみると、1964年の係数とその標準誤差がおかしなことになっています。multiplotでplot = FALSEとし、数値を確認してみます。
> (multi_01 <- multiplot(results_01, coefficients = "Raceblack", plot = FALSE)) Value Coefficient HighInner LowInner HighOuter LowOuter Model 1 0.07119541 Raceblack 0.6297813 -0.4873905 1.1883673 -1.045976 1948 2 -1.68490828 Raceblack -1.3175506 -2.0522659 -0.9501930 -2.419624 1952 3 -0.89178359 Raceblack -0.5857195 -1.1978476 -0.2796555 -1.503912 1956 4 -1.07674848 Raceblack -0.7099648 -1.4435322 -0.3431811 -1.810316 1960 5 -16.85751152 Raceblack 382.1171424 -415.8321655 781.0917963 -814.806819 1964 6 -3.65505395 Raceblack -3.0580572 -4.2520507 -2.4610605 -4.849047 1968 7 -2.68154861 Raceblack -2.4175364 -2.9455608 -2.1535242 -3.209573 1972 8 -2.94158722 Raceblack -2.4761518 -3.4070226 -2.0107164 -3.872458 1976 9 -3.03095296 Raceblack -2.6276580 -3.4342479 -2.2243631 -3.837543 1980 10 -2.47703741 Raceblack -2.1907106 -2.7633642 -1.9043839 -3.049691 1984 11 -2.79340230 Raceblack -2.4509285 -3.1358761 -2.1084548 -3.478350 1988 12 -2.82977980 Raceblack -2.4609290 -3.1986306 -2.0920782 -3.567481 1992 13 -4.45297040 Raceblack -3.4433048 -5.4626360 -2.4336392 -6.472302 1996 14 -2.67827784 Raceblack -2.2777557 -3.0788000 -1.8772336 -3.479322 2000
何と1964年の回帰係数は文字通り桁が違っており、誤差も100倍といったスケールで異なっています。
では、このような結果が得られたときにどうしたら良いでしょうか？その答えの一つが、他の年の結果などから「係数はこの辺りにあるだろう」といった事前知識をモデルに取り込むことです。弱情報事前分布をモデルに取り入れることは、ベイズの観点からは縮小と見做すことができ、Bayesian Shrinkage（ベイズ縮小）と呼ばれます。
Bayesian GLMによるフィッティング
では実際にやってみましょう。関数は{arm}ライブラリのbayesglmで、glmと同じように使えますが、事前分布を指定することが可能です。まずは先ほどのモデルに、回帰係数の事前分布として尺度パラメータ2.5のコーシー分布（＝自由度1のt分布）を用いて当てはめてみます。
results_04 <- ideo %>% split(.$Year) %>% # {arm}はlibrary(arm)とせず、arm::bayesglmとして直接呼び出す map(~arm::bayesglm(Vote ~ Race + Income + Gender + Education, data = ., family = binomial(link = "logit"), prior.scale = 2.5, prior.df = 1))
結果の確認
1964年の異常値が解消されたことを確認してみましょう。
multiplot(results_04, coefficients = "Raceblack", secret.weapon = TRUE) + coord_flip()
係数とその誤差がかなり改善されていますね！
事前分布の仮定が他の年の係数に影響を与えてはいないでしょうか？それも確認してみましょう。
multi_04 <- multiplot(results_04, coefficients = "Raceblack", plot = FALSE) > cbind(sprintf("%.3f", multi_01$Value), sprintf("%.3f", multi_04$Value)) [,1] [,2] [1,] "0.071" "0.073" [2,] "-1.685" "-1.635" [3,] "-0.892" "-0.867" [4,] "-1.077" "-1.050" [5,] "-16.858" "-5.100" [6,] "-3.655" "-3.524" [7,] "-2.682" "-2.654" [8,] "-2.942" "-2.864" [9,] "-3.031" "-2.971" [10,] "-2.477" "-2.446" [11,] "-2.793" "-2.744" [12,] "-2.830" "-2.779" [13,] "-4.453" "-4.151" [14,] "-2.678" "-2.620"
1964年は5行目ですが、その他の年の係数に大きな変動はないようです。また異常値が解消されたことで係数のバラつきはかなり小さくなりました。
plot(cbind(multi_01$Value, multi_04$Value))
> apply(cbind(multi_01$Value, multi_04$Value), 2, sd) [1] 4.037096 1.335844
また異常値となっていた５行目を除外すると、相関係数はほぼ1となります。
> correlate(cbind(multi_01$Value[-5], multi_04$Value[-5])) %>% fashion() rowname V1 V2 1 V1 1.00 2 V2 1.00
plot(cbind(multi_01$Value[-5], multi_04$Value[-5]))
おまけ
bayesglmのオプションにあるprior.scaleとprior.dfを変更することで様々な事前分布を仮定することができます。ここをどちらもInfとすると正規分布となり、glmの結果と一致するそうです。
ということで試してみたのですが、Yearごとのデータで当てはめるとエラーとなってしまったので、全体のデータに対して当てはめた結果を比較してみます。
results_05 <- glm(Vote ~ Race + Income + Gender + Education, data = ideo, family = binomial(link = "logit")) results_06 <- arm::bayesglm(Vote ~ Race + Income + Gender + Education, data = ideo, family = binomial(link = "logit"), prior.scale = Inf, prior.df = Inf) # 小数点以下７桁ぐらいまでは一致する > cbind(sprintf("%.7f", coef(results_05)), sprintf("%.7f", coef(results_06))) [,1] [,2] [1,] "0.1592985" "0.1592984" [2,] "-0.2833290" "-0.2833290" [3,] "-2.4806024" "-2.4806025" [4,] "-0.8534029" "-0.8534029" [5,] "-0.3919568" "-0.3919568" [6,] "-0.5816898" "-0.5816898" [7,] "-0.3790876" "-0.3790876" [8,] "-0.0601803" "-0.0601803" [9,] "0.0230734" "0.0230734" [10,] "0.1551473" "0.1551473" [11,] "0.7373030" "0.7373030" [12,] "0.1753482" "0.1753482" [13,] "0.0954255" "0.0954255" [14,] "0.3491240" "0.3491240" [15,] "-0.3277760" "-0.3277760" [16,] "-0.1064615" "-0.1064615" [17,] "0.1802366" "0.1802366" [18,] "-0.1227574" "-0.1227574"
> correlate(cbind(coef(results_05), coef(results_06))) %>% fashion() rowname V1 V2 1 V1 1.00 2 V2 1.00
plot(cbind(coef(results_05), coef(results_05)))
ほとんど一致していますね！
終わりに
この記事では、bayesglm関数を使ったBayesian GLMの当てはめについて紹介しました。また事前分布を変更することでbayesglmの結果がglmと一致することを確認しました。
ベイズでは事前分布のバラつきの程度によってパラメータに対する分析者の「確信度」が表されます。この確信度とデータを観測した時の影響について、追ってブログで紹介したいと思っています。
