スタイン推定量を確かめてみる
いきなりですが、最近購入したベイズモデリングの世界を読んでいて非常に面白い話題を発見しました。
- 作者: 伊庭幸人
- 出版社/メーカー: 岩波書店
- 発売日: 2018/01/18
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (4件) を見る
P118から始まるスタイン推定量についてなのですが、その直感的な意味について本文から引用すると:
それぞれのを、全部の{}の平均値の方向にだけ引っ張ってやる
ことで、そのものを推定量とするよりも、パラメータの真値との二乗誤差の期待値が小さくなるそうです。これはいわゆる縮小推定量の一種で、引っ張りの程度をデータから適応的に求めるのがこの推定量のキモなのだそうです。
スタイン推定量の定義自体はそれほど難しいものでもなく、以下のような式によって表されます:
ここで、
です。この記事では、スタイン推定量が本当によりも二乗誤差が小さくなるのか、以下のように実験してみます。
実験内容
内容としては非常に簡単なもので、平均と分散が共に既知である正規分布からいくつかのサンプルを生成します。そのサンプルをそのまま用いた場合と、スタイン推定量を用いた場合とで、真のパラメータとの二乗誤差の期待値:
]
にどれほど差が生じるかを確認しています。
パラメータ設定
実験では、一回の試行につき平均が{1, 2, ..., n
}、分散が1
である正規分布に従うサンプルをn
個生成することとしました。このn
個のデータから各サンプルのスタイン推定量を求め、真値である{1, 2, ..., n
}との二乗誤差を計算します。同時にデータそのものを用いた場合の二乗誤差も計算しておきます。上記の実験を1000
回繰り返し、最終的にどちらが小さい値となるかを比較しました。
n <- 10 mu <- 1:n sig <- 1 K <- 1000 res <- matrix(0, K, 2)
実験開始
上記のパラメータに基づき、以下のように計算を行います。
set.seed(123) for (i in 1:K) { ## サンプルを発生させる tmp <- rnorm(n, mu, sig) ## スタイン推定値を求める s2 <- (1/(n-3)) * sum((tmp - mean(tmp))^2) a <- sig / s2 s_i <- (1-a)*tmp + (a/n)*sum(tmp) ## 結果を保存する res[i, 1] <- sum((s_i - mu)^2)/n # スタイン推定値を用いた二乗誤差の期待値 res[i, 2] <- sum((tmp - mu)^2)/n # データそのものを用いた二条誤差の期待値 }
結果
それでは結果です。
> mean(res[, 1] / res[, 2]) [1] 0.9618857 > sum(res[, 1] > res[, 2]) [1] 310
スタイン推定量を用いた場合、データそのものの場合と比較して二乗誤差が96%と小さくなっていることが確認できます。
終わりに
今回の実験は以上で終了です。ではこれの何が面白いのかと言うと、スタイン推定量と階層ベイズモデルとの関連性についてなんですね。現在、別に書こうとしている記事で階層ベイズモデルに触れているのですが、その説明として
パラメータを生成する分布を仮定することで、パラメータに緩やかな縛りをかける
といった言い回しを考えています。ここでふと「階層ベイズでは(緩やかな縛りをかけることで)推定値にバイアスを与えているわけだが、それにも関わらずモデルとして良くなるというのは、一体どういう意味なのだろう?」という疑問が浮かんできました。バイアスを与えているわけなのだから、誤差が大きくなってモデルとしては劣化しそうな気がしませんか?
その時にこのスタイン推定量のことを思い出しました。つまり、階層ベイズではパラメータに対する分布を仮定することで推定値にバイアスを与えているが、スタイン推定量の性質を以って全体的には誤差を小さくすることができているということではないかと1。
書籍にはスタイン推定量が誤差を減少させる証明も付いていますが、まだ読み込めていませんので上記の理解が間違っているかもしれません。しかしこの推定量自体も非常に面白く、紹介したいと思ったので実験してみた次第です。
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
また異常値となっていた5行目を除外すると、相関係数はほぼ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) # 小数点以下7桁ぐらいまでは一致する > 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と一致することを確認しました。
ベイズでは事前分布のバラつきの程度によってパラメータに対する分析者の「確信度」が表されます。この確信度とデータを観測した時の影響について、追ってブログで紹介したいと思っています。
ベイジアンネットワークで変数間の関連性を見る
実務でデータ分析を行うときにしばしば変数間の因果関係の有無およびその強さについて問われることがあるのですが、その回答としてベイジアンネットワークを使いたいと思っていたので調べてみました。
この記事ではベイジアンネットワークを扱うためのパッケージ ({bnlearn}および{BNSL})を用いた構造の推定およびそのTipsを示します。
分析環境の準備
ライブラリの読み込み
まずは必要なライブラリを読み込みます。今回はRにおけるベイジアンネットワークのメジャーなライブラリである{bnlearn}に加え、{BNSL}を用います。
library(bnlearn) library(BNSL)
サンプルデータの作成
続いてこちらの記事を参考にサンプルデータを作成します。一部の係数を修正しましたが、記事同様にHeightとSBPがBMIに刺さり、BMIがFBSに刺さる構造となっています。
set.seed(123) norm <- rnorm(4000) Data <- matrix(norm, nrow = 1000, ncol = 4, byrow = T) colnames(Data) <- c("Height", "BMI", "SBP", "FBS") Data2 <- Data <- as.data.frame(Data) Data2$Height <- 170 + Data$Height * 10 Data2$SBP <- 120 + Data$SBP * 10 Data2$BMI <- 22 + Data$Height * 5 + Data$SBP * 5 Data2$FBS <- 90 + (Data2$BMI - 22) * 5 + Data$FBS * 10
以下のようなデータが得られます。
> head(Data2) Height BMI SBP FBS 1 164.3952 26.99116 135.5871 115.66090 2 171.2929 24.95102 124.6092 92.10449 3 163.1315 24.68614 132.2408 107.02886 4 174.0077 21.22465 114.4416 103.99239 5 174.9785 27.99603 127.0136 115.25225 6 159.3218 11.53086 109.7400 30.36538
分析実行その1
構造学習
それではベイジアンネットワークによる関連性の推定に移ります。ベイジアンネットワークはデータから変数間の関連性を推定する構造学習と,推定された構造を元にパラメータを推定するフィッティングの二段階で実施されます。
今回、構造学習のアルゴリズムにはGrow-Shrinkを指定しました。Webで調べた限りだとHill-Climbingが使われる事も多い様子ですが、先ほど参照した記事によるとGrow-Shrinkはグラフ理論による裏付けがあるとのこと。関数はgsです(Hill-Climbingはhc)。
str_01_gs <- gs(Data2) str_01_hc <- hc(Data2) par(mfrow = c(1, 2)) plot(str_01_gs) plot(str_01_hc)
gsでは設定した通りの構造が推定されましたが、hcではBMIとSBPの関係性が逆転しています。これ以降、基本的にgsを使うことにします。
フィッティング
ではこの構造を用いてパラメータを推定しましょう。なおここで得られるパラメータ推定値は個別にlmをかけたものと同じとなり、例えばBMIであればlm(BMI ~ Height + SBP)と等しい値となります。
fit_01 <- bn.fit(str_01_gs, Data2) coef_01 <- coefficients(fit_01) > coef_01 $Height (Intercept) 170.091 $BMI (Intercept) Height SBP -123.0 0.5 0.5 $SBP (Intercept) 120.1408 $FBS (Intercept) BMI -20.240828 5.019595
> coef(lm(BMI ~ Height + SBP, Data2)) (Intercept) Height SBP -123.0 0.5 0.5
同じ値となっていますね。
分析実行その2
構造学習
次にサンプルデータのgaussian.testを用いてパラメータを推定します。元の構造はわかりませんが、gsによって推定した形は以下のようです。
str_02 <- gs(gaussian.test) plot(str_02)
フィッティング
この構造を用いてパラメータの推定を行おうとするとエラーが出ます。どうやらグラフの構造が悪く、フィッティングできないようです。
> fit_02 <- bn.fit(str_02, gaussian.test) bn.fit(str_02, gaussian.test) でエラー: the graph is only partially directed.
そこでstr_02の中身を確認すると、modelが[partially directed graph]となっていることがわかります。またarcsのundirected arcsも1となっています。
> str_02 Bayesian network learned via Constraint-based methods model: [partially directed graph] nodes: 7 arcs: 7 undirected arcs: 1 directed arcs: 6 average markov blanket size: 4.00 average neighbourhood size: 2.00 average branching factor: 0.86 learning algorithm: Grow-Shrink conditional independence test: Pearson's Correlation alpha threshold: 0.05 tests used in the learning procedure: 144 optimized: FALSE
再度str_02の構造を確認するとノードBD間のエッジが無向になっていることがわかります。
そこでこの構造にB -> Dのエッジを設定してみましょう。ベイジアンネットワークではこのように、試行錯誤の結果や分析者の知見を取り込みながら、データ全体の構造を推定していきます。
この修正によりstr_03のundirected arcsが0となり、modelがpartially directed graphではなくノード間の関連性を表したものとなっていることを確認しましょう。
> (str_03 <- set.arc(str_02, "B", "D")) Bayesian network learned via Constraint-based methods model: [A][B][E][G][C|A:B][D|B][F|A:D:E:G] nodes: 7 arcs: 7 undirected arcs: 0 directed arcs: 7 average markov blanket size: 4.00 average neighbourhood size: 2.00 average branching factor: 1.00 learning algorithm: Grow-Shrink conditional independence test: Pearson's Correlation alpha threshold: 0.05 tests used in the learning procedure: 144 optimized: FALSE
ちゃんとB -> D間の矢印が有向になっています。
> plot(str_03)
では、この構造を用いて再度フィッティングしてみましょう。
> fit_03 <- bn.fit(str_03, gaussian.test) > (coef_03 <- coefficients(fit_03)) $A (Intercept) 1.007493 $B (Intercept) 2.039499 $C (Intercept) A B 2.001083 1.995901 1.999108 $D (Intercept) B 5.995036 1.498395 $E (Intercept) 3.493906 $F (Intercept) A D E G -0.006047321 1.994853041 1.005636909 1.002577002 1.494373265 $G (Intercept) 5.028076
今度は上手くいったようですね!
分析実行その3
続いて直接効果と間接効果の分離を試みます。再度gsで構造を推定し、B -> Dのエッジと共にB -> Fのエッジも追加しましょう。これでBはFに対して直接刺さるエッジと、Dを経由してFに刺さるエッジを持つことになります。
str_04 <- set.arc(str_03, "B", "F") par(mfrow = c(1, 3)) plot(str_02) # gsで推定した構造 plot(str_03) # B -> Dに修正 plot(str_04) # B -> Fを追加
フィッティング
続いてフィッティングです。B -> Fのエッジの有無で係数がどのように変化するかを比較してみましょう。必要な部分だけ抜粋します。
fit_04 <- bn.fit(str_04, gaussian.test) > coefficients(fit_04) ~ 省略 ~ $F (Intercept) A B D E G -0.4123157 1.9947919 -0.1025991 1.0737532 1.0023721 1.4943286 ~ 省略 ~
> coefficients(fit_03) ~ 省略 ~ $F (Intercept) A D E G -0.006047321 1.994853041 1.005636909 1.002577002 1.494373265 ~ 省略 ~
切片に多少の変化があっただけで、係数としては大きな変化はないようです。どうやらBの変動による寄与はほとんどDを経由してFに影響しており、BのFへの直接の効果はなさそうです。
分析実行その4
構造学習
今度は構造の推定に{BNSL}を試してみましょう。
str_05 <- bnsl(gaussian.test) plot(str_05)
どうもgsを用いた時と比べて構造がかなり異なるようですね。比較してみましょう。
par(mfrow = c(1, 2)) plot(str_02) plot(str_05)
"A -> C"や"G -> F"など、方向性が逆転する関係が見られる上、Eも独立してしまいました。やはりデータのみから因果関係を推定するのは不安定ということになるのでしょうか。
フィッティング
この構造でフィッティングしてみます。フィッティング自体はbnlearn::bn.fitを用いています(結果は省略)。
fit_05 <- bn.fit(str_05, gaussian.test)
分析実行その5
さらに続いて、データを標準化した場合の結果も確認してみます。scaleでデータを標準化します。
my_gaussian <- as.data.frame(scale(gaussian.test)) str_06 <- bnsl(my_gaussian) par(mfrow = c(1, 2)) plot(str_04) plot(str_06)
何と、標準化すると構造が変わってしまいました。。。フィッティングは省略します。
分析実行その6
最後に、ライブラリ{MASS}に格納されているサンプルデータbirthwtを用いて推定してみます。なおbnlearnではintegerは使えないので、各列をnumericに変換し、簡単のために変数も一部に限定します。
構造学習
最後なので、構造学習のアルゴリズムには色々なものを指定してみます。
vars <- c("low", "race", "smoke", "ht", "ui") my_birthwt <- as.data.frame((do.call("cbind", lapply(birthwt[, vars], as.numeric)))) str_07_01 <- hc(my_birthwt) str_07_02 <- gs(my_birthwt) str_07_03 <- iamb(my_birthwt) str_07_04 <- fast.iamb(my_birthwt) str_07_05 <- inter.iamb(my_birthwt) str_07_06 <- tabu(my_birthwt) str_07_07 <- mmhc(my_birthwt) str_07_08 <- rsmax2(my_birthwt) str_07_09 <- mmpc(my_birthwt) str_07_10 <- si.hiton.pc(my_birthwt) str_07_11 <- chow.liu(my_birthwt) str_07_12 <- aracne(my_birthwt) par(mfrow = c(1, 3)) plot(str_07_01) plot(str_07_02) plot(str_07_03) par(mfrow = c(1, 3)) plot(str_07_04) plot(str_07_05) plot(str_07_06) par(mfrow = c(1, 3)) plot(str_07_07) plot(str_07_08) plot(str_07_09) par(mfrow = c(1, 3)) plot(str_07_10) plot(str_07_11) plot(str_07_12)
アルゴリズムによって少しずつ推定された構造が違うのですが、low -> smokeおよびrace -> smokeのエッジは概ね共通して見られるようです。しかし残念ながらこのデータは「出生児の低体重の有無」と「母体環境」であるため、エッジの向きは逆ですね。。。
WAICを計算してみる
ある統計モデルの予測精度を表す指標としてAIC(Akaike Information Criterion)というものがよく使われますが、stanを用いてモデリングした場合にはAICを計算することができません。しかし2009年に渡辺澄夫先生が発表したWAIC(Widely Applicable Information Criterion、広く使える情報量規準)はAICをより一般化したものとして定義され、MCMCによる事後分布から計算することができます。
現在のところ{rstan}では直接WAICを計算することができませんが、{loo}というパッケージを用いることでWAICを得ることができます。この記事では、stanの結果からWAICを計算により求め、それがlooの結果と一致することを確認します。
WAICの定義
まず始めにWAICの定義を確認しましょう*1。WAICは以下のように定義されます:
ここで
であり、また
です。およびは、それぞれ観測値iの事後分布から得られる対数尤度の平均および対数尤度の分散を表してます。なおlppdはlog-pointwise-predictive-densityです。
シミュレーション
簡単な数値例を発生させ、そのWAICを計算してみましょう。パラメータの推定には当然ながら{rstan}を利用します。
# シミュレーションデータの発生 set.seed(123) N <- 100 # サンプルサイズ b <- 1.2 # 回帰係数 X <- rnorm(N, 0, 1) # 説明変数 E <- rnorm(N, 0, 2) # 誤差項 Y <- b * X + E D <- data.frame(Y, X) # データフレーム
回帰係数1.2、誤差標準偏差2として発生させたデータに対し、{rstan}を用いて線形回帰モデルを当てはめてみましょう。いつものように{rstan}をロードします。
library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())
stanに渡すデータを以下のように定義します。今回の例では単純な線形モデルを用いているので、渡すパラメータもほとんどありません。
dat_stan <- list( N = N, X = D$X, Y = D$Y )
定義したdat_stanをstanに渡し、フィッティングします。
fit_01 <- stan(file = './WAIC.stan', data = dat_stan, iter = 3000, chains = 4, seed = 1234)
ここでWAIC.stanの内容は以下の通りです。簡単です。
data { int<lower=1> N; vector[N] X; vector[N] Y; } parameters { real b0; real b1; real<lower=0> sigma; } model { Y ~ normal(b0 + b1 * X, sigma); }
フィッティングが終わると、以下のような結果が得られます:
> summary(fit_01)$summary mean se_mean sd 2.5% 25% b0 -0.201859 0.002533264 0.1962258 -0.5791061 -0.3360031 b1 1.092886 0.002755139 0.2134122 0.6655672 0.9504853 sigma 1.964687 0.001920081 0.1421800 1.7080842 1.8637401 lp__ -116.169685 0.021008195 1.2014656 -119.2035323 -116.7161937 50% 75% 97.5% n_eff Rhat b0 -0.2018323 -0.06931008 0.1877817 6000.000 1.000041 b1 1.0919604 1.23311123 1.5082127 6000.000 1.000551 sigma 1.9570918 2.05945046 2.2612886 5483.250 1.000177 lp__ -115.8524219 -115.28245822 -114.7994731 3270.734 1.000839
色々と書いてあって情報が多いので、必要なところを抜き出しましょう。真の回帰係数および誤差標準偏差はそれぞれ1.2と2でした。
> summary(fit_01)$summary[c("b1", "sigma"), c("mean", "50%")] mean 50% b1 1.092886 1.091960 sigma 1.964687 1.957092
回帰係数が少し低いようですが、まぁ良しとしましょう。
WAICの計算
それではこの結果を用いて、WAICを計算してみましょう。まずはlppdからですが、lppdは前述の通り以下によって計算されます:
まずはstanの結果からposterior sampleを得ます。stanを実行するときにwarmupを指定していないためsampleサイズ(ite)の半分が切り捨てられた結果、残る6,000サンプル(1,500 * 4 chain)がposterior sampleとして保存されています。
post_samples <- extract(fit_01) > str(post_samples) List of 4 $ b0 : num [1:6000(1d)] -0.146 -0.13 -0.195 -0.535 -0.162 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL $ b1 : num [1:6000(1d)] 1.036 0.928 1.004 1.309 1.005 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL $ sigma: num [1:6000(1d)] 2.2 2.01 1.64 2.03 1.86 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL $ lp__ : num [1:6000(1d)] -116 -115 -118 -117 -115 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL
この6,000サンプルのそれぞれをパラメータとして与えた時の各の対数尤度を求めると:
Des <- cbind(1, X) # 計画行列(Design Matrix) B <- cbind(post_samples$b0, post_samples$b1) # パラメータ(Beta) tmp <- matrix(NA, length(post_samples$b0), N) # 6000行、 100列の行列 for (i in 1:N) { tmp[, i] <- dnorm(Y[i], mean = B %*% Des[i, ], sd = post_samples$sigma) }
となります。ここからlppdは:
lppd <- sum(log(colMeans(tmp)))
で得られます。
同じようにを求めてみましょう。の定義は以下のようでした:
ここから、
pwaic <- sum(apply(tmp, 2, var))
として求められます。
最後に、これらの値を用いてWAICを求めます:
WAIC <- -2 * (lppd - pwaic)
これでようやくWAICを得ることがことができましたので、これらの値を出力してみると:
> cat("WAIC:", WAIC, "\n", + "lppd:", lppd, "\n", + "pwaic:", pwaic, "\n") WAIC: 420.8022 lppd: -207.2305 pwaic: 3.170572
となります。ちなみに同じモデルについて、lmを用いてフィッティングした時のAICを求めると:
fit_02 <- lm(Y ~ X, D) > AIC(fit_02) [1] 420.4523
であり、非常に近い値となります。
パッケージ{loo}を使ってWAICを求める。
では、今度は上記の行程で求めた数値をパッケージで得た値と比較してみましょう。WAICの計算には{loo}を用います。
library(loo)
looによりWAICを計算させるためには、stanのスクリプトにgenerated quantitiesブロックを追加しておく必要があります。以下をmodelブロックの後に追加します。
generated quantities { vector[N] log_lik; for (n in 1:N) log_lik[n] = normal_lpdf(Y[n] | b0 + b1 * X[n], sigma); }
normal_lpdfというのはRで言うところのdnormですね。
このスクリプトをWAIC_02.stanとして保存し、stanを実行します。
fit_03 <- stan(file = './WAIC_02.stan', data = dat_stan, iter = 3000, chains = 4, seed = 1234)
結果は以下のようになります。前の結果と同じであることを確認してください。
> summary(fit_03)$summary[c("b1", "sigma"), c("mean", "50%")] mean 50% b1 1.092886 1.091960 sigma 1.964687 1.957092
では、今度は{loo}のextract_log_likを使ってposterior sampleを抽出し、loo::waicによりWAICを計算します:
tmp02 <- extract_log_lik(fit_03) WAIC_02 <- waic(tmp02) > cat("WAIC by loo:", WAIC_02$waic, "\n", + "WAIC by manual:", WAIC, "\n") WAIC by loo: 420.8022 WAIC by manual: 420.8022
ちゃんと同じ値が得られていますね!
終わりに
分析内容は以上となります。この記事ではWAICの定義を確認した上で、その計算のための方法を示し、既存のパッケージによる数値と一致することを確かめました。WAICは階層モデルのような複雑なものであっても適用可能で、stanを使っていく上では抑えておきたい指標なので今後も勉強を続けていきます。
*1:Statistical RethinkingのChapter 6から
Stanで疑似相関を見抜く
ベイズ統計をちゃんと勉強しようと思い、最近Statistical Rethinkingという本を読んでいます。
- 作者: Richard McElreath
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 2016/02/19
- メディア: ハードカバー
- この商品を含むブログを見る
その中でタイトルの通り面白い話があったので紹介がてら実行してみます。以下は5章1節からの引用ですが、x_spurとx_realはそれぞれyと相関の「ない」変数と相関の「ある」変数を表しています:
But when you include both x variables in a linear regression predicting y, the posterior mean for the association between y and x_spur will be close to zero, while the comparable mean for x_real will be closer to 1.
この記事の内容を端的に言えば「yと相関のあるx1と、x1と相関のあるx2(つまりyとは直接相関がない)が説明変数として与えられている時、stanを用いることで不要な変数(この場合x2)を除外できる」というものになります。
またこちらの本では全体を通して{rethinking}パッケージによってパラメータの推定が行われていますが、ここではせっかくなので{rstan}でやってみます。
前準備
まずは{rstan}を読み込みます。
library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())
続いてy、x1、x2をそれぞれ以下のように定義します。
## 乱数のシードの設定 set.seed(15) N <- 100 # サンプルサイズ x1 <- rnorm(N) # yと相関のある(真の)説明変数 x2 <- rnorm(N, x_real) # x1と相関することでyと相関する(偽の)説明変数 y <- rnorm(N, x_real) # x1と相関する目的変数 ## データフレームにする dat <- data.frame(y, x_real, x_spur)
定義からx2はyと相関しないため、本来であれば統計モデルから除外されるべきものです。しかし通常はデータを取得した段階で上記のような情報は得られないため、モデリングした結果を見て判断することになります。ちょっとやってみましょう。
## 回帰分析の実行(一部省略) > summary(res.lm <- lm(y ~ ., data = dat)) Call: lm(formula = y ~ ., data = dat) Residuals: Min 1Q Median 3Q Max -4.1262 -0.5336 0.0269 0.6975 1.8557 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.10912 0.10485 -1.041 0.301 x_real 1.27236 0.14407 8.832 4.45e-14 *** x_spur -0.19652 0.09488 -2.071 0.041 *
このデータではx1、x2ともにyに対して有意に(P < 0.05)影響していることとなり*1、本来であれば不要であるx2をモデルから落とすことができません。当然かもしれませんが、両者ともに有意であるためstep関数などを用いた変数選択も意味がありません。
## stepによる変数選択の実行 > step(res.lm) Start: AIC=10.25 y ~ x_real + x_spur Df Sum of Sq RSS AIC <none> 104.34 10.248 - x_spur 1 4.615 108.95 12.576 - x_real 1 83.901 188.24 67.255 Call: lm(formula = y ~ x_real + x_spur, data = dat) Coefficients: (Intercept) x_real x_spur -0.1091 1.2724 -0.1965
Stanによる推定
それでは上記のデータセットを用いてstanで推定を実行した場合はどのようになるか試してみましょう。まずはstanに渡すデータを定義します。
## stanのdataブロックに渡すデータの定義 dat_stan <- list( N = nrow(dat), D = 2, X = dat[, 2:3], Y = dat[, 1] )
stanのスクリプトは、以下のような内容になります:
data { int<lower=0> N; int<lower=0> D; matrix[N, D] X; vector[N] Y; } parameters { real a; vector[D] beta; real<lower=0> sig; } model { for (d in 1:D) { beta[d] ~ normal(0, 1); } Y ~ normal(a + X * beta, sig); }
特に工夫した点などもなく、単に重回帰を推定するためのスクリプトとなっていますが、このスクリプトを用いてstan()を実行してみましょう*2。
## フィッティング fit_01 <- stan(file = './StanModel/Detect_Confound_Factor.stan', data = dat_stan, iter = 3000, chains = 4, seed = 1234)
結果はこちらとなります*3
> summary(fit_01)$summary[, c("mean", "2.5%", "50%", "97.5%", "Rhat")] mean 2.5% 50% 97.5% Rhat a -0.1055297 -0.3129441 -0.1052853 0.102064821 0.9998115 beta[1] 1.2470319 0.9587653 1.2454155 1.530978156 1.0001602 beta[2] -0.1827657 -0.3682030 -0.1828429 0.003176939 1.0005802 sig 1.0497539 0.9158756 1.0444427 1.210491329 0.9995408 lp__ -54.9435610 -58.5631720 -54.6322411 -53.160356164 1.0010219
ここでbeta[2]に注目すると95%ベイズ信用区間に0が含まれており、モデルから変数を除外することを検討できるようになりました。
おわりに
分析内容は以上となります。冒頭で紹介したStatistical Rethinkingには他にも多重共線性が生じている時の挙動など、関心を引く話題が多く散りばめられており、非常に読み応えのある本なので、しばらくこれで勉強しようと思います。
Ad-Stock効果を推定しつつ回帰を回したい④
これまでにMarketing Mix Modelingにおける広告の累積効果の推定について、以下のような記事を書いてきた。
ushi-goroshi.hatenablog.com
ushi-goroshi.hatenablog.com
ushi-goroshi.hatenablog.com
今回は遺伝的アルゴリズムを用いてパラメータの推定を行ったので、その結果を示しておく。
シミュレーションデータの作成〜GAによるフィッティングの実行
前回同様、必要なパッケージをインストールしてから読み込む。インストールするパッケージは遺伝的アルゴリズムの実装として非常にメジャーな{GA}。
install.packages("GA") library(GA)
続いてシミュレーションデータ作成用の関数を定義する。これは先日の記事と全く同じである。
simulate_y <- function(pars) { ## Simulation parameters n <- pars[1] # num of observation mu <- pars[2] # intercept var_e <- pars[3] # residual variance beta_01 <- pars[4] # regression coefficient of X1 to be esitmated lambda_01 <- pars[5] # decay rate of Ad-Stock effect of X1 beta_02 <- pars[6] # regression coefficient of X2 to be esitmated lambda_02 <- pars[7] # decay rate of Ad-Stock effect of X2 ## Create true Ad-Stock variables X_01_raw <- rgamma(n, 3) * ifelse(runif(n) > 0.7, 0, 1) X_01_fil <- stats::filter(X_01_raw, lambda_01, "recursive") X_02_raw <- rgamma(n, 2) * ifelse(runif(n) > 0.8, 0, 1) X_02_fil <- stats::filter(X_02_raw, lambda_02, "recursive") ## Create residuals error <- rnorm(n, 0, sqrt(var_e)) ## Create observations y <- mu + beta_01 * X_01_fil + beta_02 * X_02_fil + error ## Return dataset dat <- data.frame( "Y" = y, "X_01" = X_01_raw, "X_02" = X_02_raw, "X_01_Fil" = X_01_fil, "X_02_Fil" = X_02_fil, "Y_lag" = dplyr::lag(y, 1), "True_Error" = error) return(dat) }
この関数を用いてサンプルデータを発生させる。乱数のシードを固定しているので、前回と同じデータが得られる。推定対象となるパラメータも同様に前回と同じ。
set.seed(123) init_par <- array(c(100, 5, 2, 0.5, 0.6, 0.8, 0.5)) dat_Ana <- na.omit(simulate_y(init_par))
さらに、GAで用いる評価関数(目的関数)を最小二乗法で定義する。
OLS <- function(dat, pars) { X_01_Fil_E <- stats::filter(dat$X_01, pars[3], "recursive") X_02_Fil_E <- stats::filter(dat$X_02, pars[5], "recursive") Y_hat <- pars[1] + pars[2] * X_01_Fil_E + pars[4] * X_02_Fil_E SSE <- t(dat$Y - Y_hat) %*% (dat$Y - Y_hat) return(SSE) }
結果
上記で準備が整ったので、以下のようにGAを回してみる。
res_GA <- ga(type = 'real-valued', min = c(0, 0, 0, 0, 0), max = c(20, 2, 1, 2, 1), popSize = 500, maxiter = 500, names = c('Intercept', 'b_01', 'l_01', 'b_02', 'l_02'), keepBest = T, fitness = function(b) -OLS(dat_Ana, b))
{stan}で実行した時と同様に、{GA}を用いた場合では各パラメータに対してminとmaxを指定することでパラメータスペースに矩形制約をかけることができる。
結果を見てみよう。
> summary(res_GA)$solution Intercept b_01 l_01 b_02 l_02 [1,] 4.729741 0.4938017 0.6421675 0.718373 0.5559099 > init_par[c(-1, -3)] # サンプルサイズと誤差分散を除外している [1] 5.0 0.5 0.6 0.8 0.5
推定精度はかなり良さそうだ。
なおこれらの数値は{stan}で推定したものとほとんど同様であり、前回の記事と同じくb_02の数値の誤差が大きいようである。しかしサンプルデータを発生させた時の誤差によるかもしれないので、lmでフィットした時の結果も確認しておこう。
> lm(Y ~ X_01_Fil + X_02_Fil, dat_Ana) Call: lm(formula = Y ~ X_01_Fil + X_02_Fil, data = dat_Ana) Coefficients: (Intercept) X_01_Fil X_02_Fil 5.0846 0.5170 0.7277
やはり、なかなかの精度で推定できていたようだ。誤差分散も見てみる。
> -summary(res_GA)$fitness / (nrow(dat_Ana) - 6) [1] 2.064139 > init_par[3] [1] 2
誤差分散の推定値は{GA}の返り値(SSE)をサンプルサイズ - 6で除して求めた。ここで6はパラメータ数である。誤差分散も極めて精度よく推定できているようだ。
まとめ
これまでの結果と同様に、{GA}を用いたAd-Stock効果の推定はなかなか実用的であると言えそうだ。特にMarketing Mix Modelingによって広告効果の検証を行う場合には回帰係数を0以上に制約をかけたいことが多いため、パラメータの範囲を指定できる点はパルダ・モデルやoptimと比較すると大きなメリットとなるだろう*1。{GA}のインストールも簡単で、{stan}と異なり躓くポイントも少ない。
一方、{stan}同様に分析に要する時間がパルダ・モデルやoptimと比較して長いため、分析の初期段階での使用は向いていないようだ。ある程度モデルの形が見えてきた段階でパラメータを精度よく推定したい場合に用いるのが良いだろうか。
これまで様々な手法を用いてAd-Stock効果の推定を試してきたが、いずれの方法も単純なモデルでは十分に実用的と言える精度でパラメータを推定することができていた。そのため今後はより複雑・より現実を反映したデータを用いた時に、各手法の推定精度がどの程度悪化するかを検証していく必要がある。
*1:念のために言っておくと、optimでも「L-BFGS-B法」を指定することで矩形制約を与えることができる
Bayesian t-testを書いてみたい
IBM SPSS Statisticsにベイズ推論が実装されるらしい。
これ以外にも最近は身の回りでやけにBayesian t-testに関する話題を耳にすることが多いので、t-検定をベイズ化することのどこにニーズがあるのかを考えていたのだけれど、どうやらA/B testにおいて結果を途中で確認(Peeking)したり、サンプルサイズを変更したいという状況が結構多いそうだ。
古典的な、そして計画立てられたt-検定では分析の事前に必要なサンプルサイズを求めておく。通常は実験の開始以降にデータを追加したりなどしないのだが、ビジネスではA/Bテストの結果の如何によってデータを追加したり、あるいは途中経過を見て早めに施策の実施判断を行いたいと状況がよくあるため、そういったニーズにマッチしたのだろう(実際のところ伝統的なt-検定を使う場合であってもデータの追加や結果の覗き見は行われていると思うが)。
さて、Rにおけるベイジアンモデリングのためのパッケージとしては{MCMCpack}があり、より複雑なモデルを構築したいのであればStanを呼び出す方法があげられる。しかしt-検定ぐらい自分で書いてみようと思ったので、以下にスクリプトを示す。
データの準備
# パッケージの読み込み library(tidyverse) library(ggplot2)
サンプルデータとして以下のような正規分布に従う2標本を用意する。数値に意味はないので関心があれば標本サイズなどをそれぞれ変更してみると良いが、目的は標本の平均の差の分布について確認することである。
N_01 <- 100 mu_01 <- 100 sig_01 <- 10 N_02 <- 100 mu_02 <- 103 sig_02 <- 10 set.seed(123) Pop_01 <- rnorm(N_01, mu_01, sig_01) Pop_02 <- rnorm(N_02, mu_02, sig_02)
以下のようなサンプルが得られるだろう:
> Pop_01 %>% as_data_frame() %>% mutate("Group" = "01") %>% bind_rows(as_data_frame(Pop_02) %>% mutate("Group" = "02")) %>% ggplot(aes(value, fill = Group)) + geom_histogram(alpha = 0.5, position = "identity") + theme_bw()
このデータに対して古典的なt-検定を実施すると、平均値の差が0であるという帰無仮説はp = 0.44で棄却されない。ちなみにt.testはデフォルトでは等分散性を仮定しないWelchのt-検定となる。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
関数の定義
このデータに対して何らかのパラメータ(ここでは平均と標準偏差)を与えた時の尤度を求める関数を以下のように定義する:
likelihood <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dnorm(Pop_01, mu_01, sig_01, log = TRUE), dnorm(Pop_02, mu_02, sig_02, log = TRUE) ) }
データの尤度は0 ~ 1の数字の乗算となるため標本サイズが大きい場合には極めて小さな数値となり、桁あふれが発生してしまう。そのためdnormにおいてlog = TRUEとし、対数尤度の総和を求めることとする。
この関数を用いて最尤法で平均と分散を求めてみる。optimを使用するが、デフォルトは最小化なのでfnscaleを負の数として最大化を行う:
> pars <- c(100, 10, 100, 20) > optim(par = optim(par = pars, fn = likelihood, control=list(fnscale = -1))$par, + fn = likelihood, control=list(fnscale = -1)) $par [1] 100.904769 9.082697 101.923277 9.620816 $value [1] -730.8205 $counts function gradient 101 NA $convergence [1] 0 $message NULL
続いて事前分布を定義する。事前分布の設定は色々と工夫すべきポイントではあるが、簡単のため-10000 ~ 10000の一様分布とした。尤度はデータと同じく対数としてある。
prior <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dunif(mu_01, -10000, 10000, log = TRUE), dunif(mu_02, -10000, 10000, log = TRUE), dunif(sig_01, -10000, 10000, log = TRUE), dunif(sig_02, -10000, 10000, log = TRUE) ) }
事後分布は事前分布と尤度の積であるが、対数を取ってあるためここでは和となる。
posterior <- function(pars) { likelihood(pars) + prior(pars) }
MCMCの実行
以上で準備ができたので以下のようにMCMCを実行する。なおサンプラーとしてGibbs Samplerを実装したつもりだが、Gibbs Samplingにおいてパラメータの採択確率があったかは記憶が不確かなので、この辺りは後日確認しておくとしてひとまず進めることとする。サンプルの数は30,000としたが、実際には割と早く収束しているようなのでこんなには必要ない。
N_iter <- 30000 pars <- c(50, 30, 20, 50) results <- matrix(0, nrow = N_iter, ncol = 8) results[1, ] <- c(pars, 0, 0, 0, 0) for (i in 2:N_iter) { for (j in 1:4) { candidate <- pars candidate[j] <- candidate[j] + rnorm(1, sd = 0.5) ratio <- exp(posterior(candidate) - posterior(pars)) t <- runif(1) if (t < ratio) { pars <- candidate } } results[i, ] <- c(pars, ratio, t, posterior(pars), posterior(candidate)) }
以下のような結果となる。
> par(mfrow = c(2, 2)) > plot(results[, 1]) > plot(results[, 2]) > plot(results[, 3]) > plot(results[, 4])
要約統計量も見ておこう。burn-inとして最初の15,000サンプルを切り捨て、自己回帰性を小さくするために10刻みでサンプルを取得する。
> Res_MCMC <- results[seq(15000, N_iter, 10), ] > summary(Res_MCMC[, c(1:4)]) V1 V2 V3 V4 Min. : 97.79 Min. : 7.198 Min. : 98.41 Min. : 7.997 1st Qu.:100.24 1st Qu.: 8.732 1st Qu.:101.34 1st Qu.: 9.266 Median :100.86 Median : 9.191 Median :101.98 Median : 9.753 Mean :100.86 Mean : 9.233 Mean :101.98 Mean : 9.774 3rd Qu.:101.45 3rd Qu.: 9.675 3rd Qu.:102.63 3rd Qu.:10.218 Max. :104.16 Max. :11.959 Max. :105.20 Max. :12.632
単純な推定なので精度よく推定できているようである。
平均値の差の分布
では最後に、この分析の目的であった2標本の平均の差の分布について確認してみよう。
par(mfrow = c(1, 1)) hist(Res_MCMC[, 1] - Res_MCMC[, 3])
これを見ると、平均値の差の分布が比較的中心に近い位置で0を跨いでおり、伝統的なt-検定で言う所の有意であるとは言えなさそうだ。実際に平均値の差の経験分布を確認すると、95%確信区間は-3.76 ~ 1.62である。
> quantile(Res_MCMC[, 1] - Res_MCMC[, 3], seq(0, 1, 0.025)) 0% 2.5% 5% 7.5% 10% 12.5% 15% 17.5% -5.62647895 -3.76661174 -3.35462735 -3.11314731 -2.89345678 -2.70960554 -2.56386382 -2.43985754 20% 22.5% 25% 27.5% 30% 32.5% 35% 37.5% -2.30487729 -2.18744684 -2.07270632 -1.98450427 -1.89198288 -1.77273198 -1.67356108 -1.57205147 40% 42.5% 45% 47.5% 50% 52.5% 55% 57.5% -1.47598439 -1.40161516 -1.30787797 -1.21558957 -1.12424331 -1.03598088 -0.95253546 -0.85994024 60% 62.5% 65% 67.5% 70% 72.5% 75% 77.5% -0.78004294 -0.68704363 -0.59508922 -0.49850858 -0.42349299 -0.30716666 -0.19718374 -0.10822330 80% 82.5% 85% 87.5% 90% 92.5% 95% 97.5% 0.01189818 0.14609443 0.30742495 0.45296435 0.69891710 0.92152549 1.18895282 1.62042273 100% 3.58972078
改めてt-検定の結果を見てみると、95%信頼区間として同様の結果が得られていることがわかるだろう。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
おまけ
以上でt-検定のベイズ化は終わりだが、ベイズ化した分析のメリットとして以下のような疑問に答えることが可能となる:1つ目の標本の平均が2つ目の標本の平均よりも大きい確率はどの程度だろうか?
> mean(Res_MCMC[, 3] - Res_MCMC[, 1] < 0) [1] 0.2031979
上記のような質問は、古典的な仮説検定(頻度論的アプローチ)では答えることができない。仮説検定は常に「帰無仮説が棄却されるか否か」(例えば平均値と平均値に差が0であるか否か)であり、「何%の確率で棄却される」とは表現できないのだ。しかしこれをベイズ化することでより直感的な表現が可能となる。これはビジネスへの応用において確かにメリットになると感じた。