階層ベイズと状態空間モデルで広告効果を推定したい
背景
これまでMarketing Mix Modeling(MMM)におけるAdStock効果の推定について色々と記事を書いてきましたが、その他にも試したいと思っているモデルがいくつかあります。その一つが階層ベイズモデルと状態空間モデルを同時に取り扱うものです。
例えば「地域別の売上推移のデータ」が手元にあると考えてみましょう。地域ではなく人や商品でも構いませんが、ある要因の各水準がそれぞれ時系列データを持っている状況(いわゆるパネルデータ)で、ひとまずここでは地域とします。このようなデータはあらゆる会社で保有していることでしょう。 今、各地域についてMMMにより広告効果を推定することを考えたとき、どのようなモデリングが可能でしょうか?
シンプルに考えれば、地域ごとに一つずつモデルを作るという方法が挙げられます。例えば地域の数が2つ3つしかなかったり、モデルの作成に時間をかけることが可能であればこの方法は有効かもしれません。しかし問題もあります。地域ごとに独立してモデルを作成するということは、各地域のパラメータは互いに独立であるとの仮定を置くことになるのですが、それは事実でしょうか?
確かにTVや新聞には地方欄やローカル局があり、地域限定のコンテンツや提供される場合もあるでしょう。しかしどちらかと言えば全国で同じコンテンツが提供される割合の方が大きいでしょうし、購買行動に地域による差がそれほど認められるとは、経験的にも思えません。もし仮に地域による異質性というものがそれほど強くないのであれば、むしろそういった情報を積極的に活かした方が良いのではないでしょうか?
目的
そのような時に、パラメータ間に緩やかな縛りをかける方法として階層ベイズモデルという手法があります。詳しくは書籍1を参考にして頂くとして、ここでは「各地域の広告効果を表すパラメータを生成する分布を仮定する」方法と説明しておきます。この分布の幅(バラつき、分散)が十分に大きければ広告効果は地域によって様々な値を取り得ますし、逆に分布の幅が極端に狭ければ地域ごとの広告効果に差はないことを意味します。
この階層構造を持つモデルを、時系列データの分析でお馴染みの状態空間モデルと合わせて取り扱いたいというのが今回の試みです。具体的には、以下のようなモデルを想定しています。
添字のA
およびt
はそれぞれ地域(Area)と時点(time)を意味しており、は「ある地域A
のt
時点における観測値」となります。
ポイントは4行目で、地域ごとの広告効果に対して、その生成元となる分布(ここでは正規分布)を仮定しています。もしが大きければは様々な値を取り得ますが、値が小さければ非常に似通った値になることが理解できると思います。
ライブラリの読み込み
今回の分析では{rstan}
を使用します。階層ベイズ + 状態空間モデルといった非常に複雑なモデルを表現できるのが強みです。
library(tidyverse) library(ggplot2) library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE)
シミュレーションデータの生成
分析用のデータを生成しますが、検証のポイントとしては以下のようになります:
- 5地域(
num_region
)から4年分(num_year
)の月次(num_month
)データを収集し、合計240点の観測値を得たと想定 - 回帰係数(広告効果)は平均
0.5
、分散0.01
の正規分布にしたがって発生 - 状態変数のパラメータは各地域で共通
## シードの固定 set.seed(123) ## 観測値の数 num_region <- 5 num_month <- 12 num_year <- 4 data_length <- num_month * num_year ## 回帰係数 mu_beta <- 0.5 var_beta <- 0.01 beta_ad <- rnorm(num_region, mu_beta, sqrt(var_beta)) beta_ad_all <- rep(beta_ad, each = data_length) ## 分布の設定 # 状態変数の初期値 state_t0 <- 3 var_state_t0 <- 1 # 状態変数の分散 var_state <- 0.01 # 誤差項 var_error <- 0.01 ## 説明変数X scale_x <- 5 zero_per <- 0.2 X <- ifelse(runif(num_region * data_length) > zero_per, 1, 0) * rpois(num_region * data_length, scale_x)
上記の条件に従いデータを発生させます。状態変数を除けば先に示したモデルの通り簡単な線形モデルです。Area_ID
は、後ほどStanでフィッティングを行うために必要となる変数で、YM
は描画用です。
## 状態変数 State <- matrix(0, nrow = num_region * data_length) for (i in 1:num_region) { for (j in 1:data_length) { if (j == 1) { ## 1, 49, 97, 145, 193行目が各地域の先頭となる State[(i-1) * data_length + j] <- rnorm(1, state_t0, sqrt(var_state_t0)) } else { State[(i-1) * data_length + j] <- State[(i-1) * data_length + j - 1] + rnorm(1, 0, sqrt(var_state)) } } } ## 目的変数 error <- rnorm(num_region * data_length, 0, sqrt(var_error)) Y <- State + X * beta_ad_all + error ## 地域ごとのID Area_ID <- 1:num_region DF <- data.frame( "Area_ID" = rep(Area_ID, each = data_length), "YM" = rep(1:data_length, time = num_region), "Y" = Y, "X" = X, "True_s" = State, "True_e" = error)
ここで観測値と状態変数が各地域でどのように推移しているか確認しておきます。
DF %>% gather("Var", "Val", -c(Area_ID, YM)) %>% filter(Var %in% c("Y", "True_s")) %>% ggplot(., aes(x = YM, y = Val, color = Var)) + geom_line() + facet_wrap(~Area_ID)
また、この時点で地域をまとめてlmで推定した$\beta$がどのようになるかも確認しておきましょう。
> coef(lm(Y ~ X, data = DF)) (Intercept) X 3.8088184 0.5003114
また地域ごとにデータを分割した場合も試してみます。
results <- list() for (i in Area_ID) { results[[as.character(i)]] <- coef(lm(Y ~ X, data = DF, subset = Area_ID == i)) }
> print(cbind(do.call("rbind", results), beta_ad)) (Intercept) X beta_ad 1 5.066906 0.4292258 0.4439524 2 2.700195 0.4837768 0.4769823 3 3.764447 0.6628493 0.6558708 4 2.953617 0.4988667 0.5070508 5 4.051693 0.5314325 0.5129288
そこそこ近い値が推定されていますね(汗
Stanによるフィッティング
気を取り直して、Stanを用いたフィッティングを実行してみます。Stanに渡すデータは以下の通りです。
dat_Stan <- list(N = data_length, K = num_region, Y = DF$Y, X = DF$X, Area_ID = DF$Area_ID)
またStanのスクリプトは以下のようになります。
data { int N; // 地域ごとの観測値の数(data_length) int K; // 地域の数(num_region) vector[N*K] Y; // 観測値のベクトル vector[N*K] X; // 説明変数のベクトル int<lower=1, upper=K> Area_ID[N*K]; } parameters { real state_t0; // 状態変数の0期目の平均 vector[N*K] state; // 状態変数のベクトル real beta_0; // 回帰係数の事前分布の平均 vector[K] beta; // 地域ごとの回帰係数のベクトル real<lower=0> var_state_t0; // 状態変数の0期目の分散 real<lower=0> var_state; // 状態変数の分散 real<lower=0> var_beta_0; // 回帰係数の事前分布の分散 real<lower=0> var_error; // 誤差分散 } model { // 状態変数をサンプリング for (k in 1:K) { // 1期目の値は0期目の分布からサンプルする state[1 + (k-1)*N] ~ normal(state_t0, var_state_t0); // 2期目以降は前期の値を平均とした分布からサンプルする for(i in 2:N) { state[i + (k-1)*N] ~ normal(state[i-1 + (k-1)*N], var_state); } } // 回帰係数をサンプリング for (k in 1:K) { beta[k] ~ normal(beta_0, var_beta_0); } // Yをサンプリング for(i in 1:(N*K)) { Y[i] ~ normal(state[i] + beta[Area_ID[i]] * X[i], var_error); } }
このスクリプトをHB_SSM_Sim.stan
として保存し、フィッティングを行います。
fit_01 <- stan(file = '/YourDirectory/HB_SSM_Sim.stan', data = dat_Stan, iter = 1000, chains = 4, seed = 123)
結果の確認
フィッティングが終わったので、結果を見てみましょう。fit_01
からサンプルを抽出して加工します。
## サンプルを抽出する res_01 <- rstan::extract(fit_01) ## 該当するパラメータを取り出す ests <- summary(fit_01)$summary t0_rows <- rownames(ests)[grep("state_t0", rownames(ests))] state_rows <- rownames(ests)[grep("state\\[", rownames(ests))] b0_rows <- rownames(ests)[grep("beta_0", rownames(ests))] beta_rows <- rownames(ests)[grep("beta\\[", rownames(ests))] ## 状態変数 state_par <- ests %>% data.frame %>% select(mean) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% state_rows) %>% mutate("Area" = DF$Area_ID) state_cmp <- data.frame( True = State, Est = state_par$mean, Area = as.factor(state_par$Area), YM = DF$YM )
まずは状態変数の推定結果を見てみましょう。実際の値と推定値を並べてみます。 (なお以降の作業の前にサンプルのtraceプロットとヒストグラムを確認し、収束していると判断しています)
state_cmp %>% gather("Var", "Val", -c(Area, YM)) %>% ggplot(., aes(x = YM, y = Val, colour = Var)) + geom_line() + facet_wrap(~Area)
おぉ、かなり精度良く推定できているようですね!状態変数の初期値や分散の推定値はどうでしょうか?
> ests %>% data.frame %>% select(mean) %>% rename("Estimated" = mean) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% t0_rows) %>% mutate("True" = c(state_t0, var_state_t0)) %>% mutate("Simulated" = c( mean(DF[c(1, 49, 97, 145, 193), "True_s"]), var(DF[c(1, 49, 97, 145, 193), "True_s"]) )) %>% select(Par, True, Simulated, Estimated) Par True Simulated Estimated 1 state_t0 3 3.7429556 3.731479 2 var_state_t0 1 0.8521877 1.328192
状態変数初期値の平均および分散の真の値がそれぞれ3
および1
であったのに対し、生成されたデータでは3.7
および0.85
でした。推定された値は3.7
および1.3
で、分散がやや過大に推定されているようです。ただし以下に示すように、推定値の95%信用区間も結構広いため、逸脱しているとまでは言えないようです。
> ests %>% data.frame %>% select(X2.5., X97.5.) %>% rename("Lower95" = X2.5., "Upeer95" = X97.5.) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% t0_rows) %>% select(Par, everything()) Par Lower95 Upeer95 1 state_t0 2.3975622 5.027106 2 var_state_t0 0.5549254 3.525860
続いて回帰係数はどうでしょうか。
beta_par <- ests %>% data.frame %>% select(mean) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% beta_rows) beta_cmp <- data.frame( True = beta_ad, Est = beta_par$mean ) ggplot(beta_cmp, aes(x = True, y = Est)) + geom_point() + coord_fixed()
こちらも、事前に設定した回帰係数と推定値が似通っているようです。数値を確認しても、良い精度で推定できています。
ests %>% data.frame %>% select(mean) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% beta_rows) %>% bind_cols("True" = beta_ad) %>% select(Par, True, mean) Par True mean 1 beta[1] 0.4439524 0.4520343 2 beta[2] 0.4769823 0.4825661 3 beta[3] 0.6558708 0.6604340 4 beta[4] 0.5070508 0.5001494 5 beta[5] 0.5129288 0.5355829
最後に、回帰係数の事前分布についても見ておきましょう。
> ests %>% data.frame %>% select(mean) %>% mutate("Par" = rownames(ests)) %>% filter(Par %in% b0_rows) %>% bind_cols("True" = c(mu_beta, var_beta)) %>% select(Par, True, mean) Par True mean 1 beta_0 0.50 0.5184981 2 var_beta_0 0.01 0.1385978
回帰係数の事前分布の分散はちょっと大きく推定されているようですが、平均は近しい値となっているようです。
終わりに
以上、「階層ベイズと状態空間モデルを合わせて取り扱いたい」という試みでしたが、結果としては概ね満足の行くものになったと思います。序盤に書いた通り、階層ベイズ + 状態空間モデルのような非常に複雑なモデルであっても、stanを使えば非常に簡単に推定が可能です。
これまで広告効果の推定に関していくつか記事を書いてきましたが、実はゴールとなるモデルとしてはこれを想定していました。あとはこのモデルに、これまで紹介してきたようなAdStock効果の推定を組み込めば試してみたかったモデルは一通り完了することになります。
これらについても追って紹介したいと思います。
スタイン推定量を確かめてみる
いきなりですが、最近購入したベイズモデリングの世界を読んでいて非常に面白い話題を発見しました。
- 作者: 伊庭幸人
- 出版社/メーカー: 岩波書店
- 発売日: 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法」を指定することで矩形制約を与えることができる