ベイジアンネットワークで変数間の関連性を見る
実務でデータ分析を行うときにしばしば変数間の因果関係の有無およびその強さについて問われることがあるのですが、その回答としてベイジアンネットワークを使いたいと思っていたので調べてみました。
この記事ではベイジアンネットワークを扱うためのパッケージ （{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
分析実行その１
構造学習
それではベイジアンネットワークによる関連性の推定に移ります。ベイジアンネットワークはデータから変数間の関連性を推定する構造学習と，推定された構造を元にパラメータを推定するフィッティングの二段階で実施されます。
今回、構造学習のアルゴリズムには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
同じ値となっていますね。
分析実行その２
構造学習
次にサンプルデータの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
今度は上手くいったようですね！
分析実行その３
続いて直接効果と間接効果の分離を試みます。再度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への直接の効果はなさそうです。
分析実行その４
構造学習
今度は構造の推定に{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)
分析実行その５
さらに続いて、データを標準化した場合の結果も確認してみます。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)
何と、標準化すると構造が変わってしまいました。。。フィッティングは省略します。
分析実行その６
最後に、ライブラリ{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
- メディア: ハードカバー
- この商品を含むブログを見る
その中でタイトルの通り面白い話があったので紹介がてら実行してみます。以下は５章１節からの引用ですが、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)
サンプルデータとして以下のような正規分布に従う２標本を用意する。数値に意味はないので関心があれば標本サイズなどをそれぞれ変更してみると良いが、目的は標本の平均の差の分布について確認することである。
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
単純な推定なので精度よく推定できているようである。
平均値の差の分布
では最後に、この分析の目的であった２標本の平均の差の分布について確認してみよう。
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-検定のベイズ化は終わりだが、ベイズ化した分析のメリットとして以下のような疑問に答えることが可能となる：１つ目の標本の平均が２つ目の標本の平均よりも大きい確率はどの程度だろうか？
> mean(Res_MCMC[, 3] - Res_MCMC[, 1] < 0) [1] 0.2031979
上記のような質問は、古典的な仮説検定（頻度論的アプローチ）では答えることができない。仮説検定は常に「帰無仮説が棄却されるか否か」（例えば平均値と平均値に差が0であるか否か）であり、「何%の確率で棄却される」とは表現できないのだ。しかしこれをベイズ化することでより直感的な表現が可能となる。これはビジネスへの応用において確かにメリットになると感じた。
Random Forestの結果をSQLに落としたい
Random Forestは数ある機械学習アルゴリズムの中でも、高い精度を得やすく分散処理が容易であることから頻繁に用いられるものの一つであると思う。
Rでは{randomForest}パッケージが有名だが、最近では{ranger}や{Rborist}といった新しいパッケージが出ていたり、{edarf}のようなRandom Forestによるデータの解析や解釈を主眼にしたものが出ていたりして動きとして面白い。
そんなRandom Forestであるが、せっかくRで良いモデルを作ることができたとしても実際には使えないことがある。これは、例えばSQLに落とすことが容易ではないためシステムに組み込めないからだ。
ビジネスにおいてデータをRで分析する場面は数多くあるが、Rの結果を（文字通り）そのままシステムに埋め込むことはほとんどないと言って良いと思う。Rを使って得られたモデルは、SQLなどの他の環境に合わせて再構築した上で利用することが多いだろう。そこで今回はRの{randomForest}パッケージで作成したモデルをSQLに落とすための方法を紹介したい。
なおこの記事は基本的に以下の記事を参考に作成したが、ひとまずRの{randomForest}で作成したモデルからSQLを得るための流れを示すことが目的であったため大幅に簡略化してある。またSQLのインデントも私の好みに合わせたので、元記事も参照することをお薦めする。
https://gist.github.com/shanebutler/96f0e78a02c84cdcf558
処理の流れ
Random Forestというアルゴリズムは、大雑把に言えば簡易な決定木を大量に作成し、その結果の多数決によって最終的な判断を行うというものである。よってSQLに落とす場合も手順としては以下のようになる：
ただし実際には1と2は同時に行われるだろう。
1. 決定木の結果をSQLとして書き下す
それでは実装に移ろう。まずはSQLに書き下すところから。サンプルとしてirisを用い、簡単のためrandomForestの木の数は2としたモデルを作成する。
library(randomForest) set.seed(123) res <- randomForest(Species ~ ., data = iris, ntree = 2)
この処理で得られた結果はresオブジェクトに保存されるが、特にSQLを書くために必要となる決定木の分岐はgetTreeで取り出すことができる。取り出したい木が何番目であるかを指定し、ラベルの表示をTRUEにすることで：
> getTree(res, k = 1, labelVar = TRUE) left daughter right daughter split var split point status prediction 1 2 3 Petal.Width 0.80 1 <NA> 2 0 0 <NA> 0.00 -1 setosa 3 4 5 Petal.Width 1.65 1 <NA> 4 6 7 Petal.Width 1.35 1 <NA> 5 8 9 Sepal.Length 6.05 1 <NA> 6 0 0 <NA> 0.00 -1 versicolor 7 10 11 Petal.Length 4.95 1 <NA> 8 12 13 Petal.Length 4.85 1 <NA> 9 0 0 <NA> 0.00 -1 virginica 10 0 0 <NA> 0.00 -1 versicolor 11 0 0 <NA> 0.00 -1 virginica 12 14 15 Sepal.Length 5.40 1 <NA> 13 0 0 <NA> 0.00 -1 virginica 14 0 0 <NA> 0.00 -1 virginica 15 0 0 <NA> 0.00 -1 versicolor
このようなテーブルが得られる。決定木は親となるノードから次々に分岐していく形で木が形成されるが、このテーブルはそのためのルールを示したものとなっており、「left daughter」が該当する場合に次に参照する行を、「right daughter」が該当しなかった場合に参照する行を示す。
例えばこの場合、1行目を読むことで「Petal.Widthが0.8未満であるか」を基準に分岐が発生することがわかり、また該当する場合（left daughter）には2行目を参照することでsetosaになることを示している。
このようなテーブルを各決定木について得ることができるため、SQLに落とすための流れとしては：
- テーブルをSQLにするための処理を書く
- 決定木の数だけループを回す
ことが考えられるだろう。そこでこのテーブルをSQLに落とすための関数を以下のように定義する：
Make_SQL <- function(Rule_Table, cnt, ind = 1) { Rule <- Rule_Table[cnt, ] Indent <- paste(rep(" ", ind), collapse = "") var <- as.character(Rule[, "split var"]) val <- Rule[, "split point"] if(Rule[, "status"] != -1) { cat(paste0("\n", Indent, "CASE\n", Indent, " WHEN ", var, " <= ", val, " THEN")) Make_SQL(Rule_Table, Rule[, "left daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, " ELSE")) Make_SQL(Rule_Table, Rule[, "right daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, "END")) } else { cat(paste0(" '", Rule[, "prediction"], "'")) } }
ここでRule_Tableは先ほどのテーブルを示している。また関数の序盤で定義しているIndentはSQLを書くときのインデントを私の好みに合わせるために修正してある。
全体的に、SQLはCASE文で形成されることがわかる。これは決定木がIf THEN Elseで表現可能なルールによって作成されることからも自明だろう。なおRule_Tableのstatusには、当該ノードが末端（ターミナルノード）であるかが示されているため、ここを参照することで予測値が決まるか（それ以上の分岐がないか）がわかる。もしターミナルノードでないなら分岐した先で同様の処理を行うため、再帰的な呼び出しが必要であることがわかる。分岐した場合、該当するものについてはleft daughterの行を参照し、そうでないものについてはright daughterを参照するようになっている。
これを実行すると、例えば以下のようなSQL文が得られる：
> Make_SQL(Rule_Table, 1, 1) CASE WHEN Petal.Width <= 0.75 THEN 'setosa' ELSE CASE WHEN Petal.Width <= 1.75 THEN CASE WHEN Sepal.Width <= 2.25 THEN CASE WHEN Petal.Length <= 4.75 THEN 'versicolor' ELSE 'virginica' END ELSE CASE WHEN Petal.Width <= 1.35 THEN 'versicolor' ELSE CASE WHEN Sepal.Width <= 2.65 THEN CASE WHEN Petal.Width <= 1.45 THEN 'virginica' ELSE 'versicolor' END ELSE 'versicolor' END END END ELSE CASE WHEN Sepal.Length <= 5.95 THEN CASE WHEN Sepal.Width <= 3 THEN 'virginica' ELSE 'versicolor' END ELSE 'virginica' END END END
これをRandom Forestのモデル作成に用いた決定木の数だけ作成すれば良いということになる。なお決定木の数を大きなものとした場合にSQL文が大きくなってしまうため、実行の際には注意が必要である。
このようにして得られたSQLはこのままだと表示して終わりであるため、何らかのテーブルに格納しておく必要がある。そこでINSERT文を前に置いておくことにする。すなわち：
cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1)
このような書き方になる。ここでidは予測する場合の粒度を示すことになる。irisの場合は1行が1データ（個体）になるため行ごとに連番を振れば良いが、ここでは省略する。
これで一つの決定木についてSQLによる結果をテーブルに格納することができるようになったため、上記のスクリプトをループで回すと以下のようになるだろう：
for (i in 1:(res$ntree)) { Rule_Table <- getTree(res, k = i, labelVar = TRUE) cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1) cat(paste0(" as tree", i, "\nFROM\n", " input_data", ";\n\n")) }
res$ntreeにはモデルの作成に用いた決定木の数が格納されているため、その分だけループを回すことになる。またFROM句では決定木を当てはめるためのデータ（この例ではiris）を指定することになる。
1つのテーブルに各SQLの結果を集約する
上記の処理によりtbl_rfには各決定木の結果が格納されているだろう。これを用いて最終的な判断を行うには、各決定木の結果を集約する必要がある。以下のようになるだろう：
> cat(paste0("INSERT INTO rf_predictions\n", + "SELECT\n a.id,\n a.pred\n", + "FROM\n (\n SELECT id as id, pred, COUNT(*) as cnt \n", + " FROM tbl_rf\n GROUP BY id, pred\n ) a\n", + " INNER JOIN\n (\n ", + "SELECT id, MAX(cnt) as cnt\n", + " FROM\n", + " (\n", + " SELECT id as id, pred, COUNT(*) as cnt\n", + " FROM tbl_rf\n", + " GROUP BY id, pred\n", + " )\n", + " GROUP BY id\n ) b\n", + " ON a.id = b.id AND a.cnt = b.cnt;\n\n")) INSERT INTO rf_predictions SELECT a.id, a.pred FROM ( SELECT id as id, pred, COUNT(*) as cnt FROM tbl_rf GROUP BY id, pred ) a INNER JOIN ( SELECT id, MAX(cnt) as cnt FROM ( SELECT id as id, pred, COUNT(*) as cnt FROM tbl_rf GROUP BY id, pred ) GROUP BY id ) b ON a.id = b.id AND a.cnt = b.cnt;
irisのSpeciesを用いる場合はclassificationとなるため、各決定木の判断を集約し、最も多かったクラスを最終的な予測とする。上記はそのためのクエリとなっていることがわかる。
最後に、これまでのスクリプトを一つにまとめたものを以下に示す：
library(randomForest) set.seed(123) res <- randomForest(Species ~ ., data = iris, ntree = 2) Make_SQL <- function(Rule_Table, cnt, ind = 1) { Rule <- Rule_Table[cnt, ] Indent <- paste(rep(" ", ind), collapse = "") var <- as.character(Rule[, "split var"]) val <- Rule[, "split point"] if(Rule[, "status"] != -1) { cat(paste0("\n", Indent, "CASE\n", Indent, " WHEN ", var, " <= ", val, " THEN")) Make_SQL(Rule_Table, Rule[, "left daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, " ELSE")) Make_SQL(Rule_Table, Rule[, "right daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, "END")) } else { cat(paste0(" '", Rule[, "prediction"], "'")) } } for (i in 1:(res$ntree)) { Rule_Table <- getTree(res, k = i, labelVar = TRUE) cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1) cat(paste0(" as tree", i, "\nFROM\n", " input_data", ";\n\n")) } cat(paste0("INSERT INTO rf_predictions\n", "SELECT\n a.id,\n a.pred\n", "FROM\n (\n SELECT id as id, pred, COUNT(*) as cnt \n", " FROM tbl_rf\n GROUP BY id, pred\n ) a\n", " INNER JOIN\n (\n ", "SELECT id, MAX(cnt) as cnt\n", " FROM\n", " (\n", " SELECT id as id, pred, COUNT(*) as cnt\n", " FROM tbl_rf\n", " GROUP BY id, pred\n", " )\n", " GROUP BY id\n ) b\n", " ON a.id = b.id AND a.cnt = b.cnt;\n\n"))
最後に
Random ForestはDeep Learningが注目されるようになった現在においても強力なアルゴリズムの一つである。またDeep Forestのように多段構造とすることでDeep Learningを上回る性能を得る場合もある。分散処理に強く、人による解釈も容易であるためビジネスにおける実用性は相当に高いアルゴリズムであるため、実装への障害によって敬遠されるのは非常に勿体無いと思う。この記事がそれらの障害を乗り越えるための一助となれば幸いである。
現状ではDeep LearningよりもRandom Forestを使いこなせる方が実用においては有益だろう。曲線的な分類に弱いなどの欠点もあるが、まだまだ発展のポテンシャルが高いアルゴリズムだと思うので今後も注目していきたい。
ggvisによるRで可視化
面白いパッケージを発見したので紹介。可視化のための機能を提供してくれるもので、{ggvis}というもの。もしかすると{ggplot2}を超えるかもしれない。
使い方が非常にシンプルで、思想として{ggplot2}に似ているので慣れた人にはスイッチしやすいと思う。dplyrユーザーにはお馴染みの%>%演算子を使用可能というのもよい。
ひとまず試してみる。
## Library installation install.packages("ggvis") ## Load library library(ggvis)
ggvisの実行。
mtcars %>% ggvis(~mpg, ~wt, fill = ~factor(cyl)) %>% layer_points()
このように可視化に用いるデータを%>%で後続の処理に渡し、後から可視化に用いるグラフを指定する形式は{ggplot2}と同じ。
では{ggvis}の何が面白いかというと、Shinyのようなインタラクティブな操作が可能なグラフを以下のように簡単に作成できること。
以下をお試しあれ。
data.frame(x = rnorm(1000)) %>% ggvis(~x, fill := input_select(c("red", "blue", "green"))) %>% layer_histograms(width = input_slider(0.01, 3, 0.1, 0.01, label="bin_width"))
連続量のヒストグラムを作成するとき、ビンの幅を決めるのは以外に難しい。R標準のパッケージである{graphics}のhistや、{MASS}のtruehistは所定のアルゴリズム（ScottとかFDとか）に基づいて自動的に決定してくれるが、実際にデータを分析する立場からするとビンを様々に変えながら眺めたいところ。
そのような時にShinyのようにインタラクティブにビン幅を変えられると嬉しいのだが、Shinyは準備に手間がかかるのでなかなか手が出ない。それを{ggvis}では先ほどのように簡単なコマンドでインタラクティブなグラフを作成できる。
他にも
mtcars %>% ggvis(x = ~wt, y = ~mpg) %>% layer_points() %>% layer_smooths(span = input_slider(0, 1, 0.5, 0.1, label = "span"))
このような使い方ができる。
このパッケージはまだ開発の途上のようで、demoで見ようとしたいくつかの機能がエラーで見れなかったり、作成したグラフが保存できなかったりする（だからこのブログでもグラフの画像を貼っていない）。
今後の修正に大きく期待できる、面白いパッケージだと思う。
