統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

階層ベイズと状態空間モデルで広告効果を推定したい

背景

これまでMarketing Mix ModelingMMM)におけるAdStock効果の推定について色々と記事を書いてきましたが、その他にも試したいと思っているモデルがいくつかあります。その一つが階層ベイズモデルと状態空間モデルを同時に取り扱うものです。

例えば「地域別の売上推移のデータ」が手元にあると考えてみましょう。地域ではなく人や商品でも構いませんが、ある要因の各水準がそれぞれ時系列データを持っている状況(いわゆるパネルデータ)で、ひとまずここでは地域とします。このようなデータはあらゆる会社で保有していることでしょう。 今、各地域についてMMMにより広告効果を推定することを考えたとき、どのようなモデリングが可能でしょうか?

シンプルに考えれば、地域ごとに一つずつモデルを作るという方法が挙げられます。例えば地域の数が2つ3つしかなかったり、モデルの作成に時間をかけることが可能であればこの方法は有効かもしれません。しかし問題もあります。地域ごとに独立してモデルを作成するということは、各地域のパラメータは互いに独立であるとの仮定を置くことになるのですが、それは事実でしょうか?

確かにTVや新聞には地方欄やローカル局があり、地域限定のコンテンツや提供される場合もあるでしょう。しかしどちらかと言えば全国で同じコンテンツが提供される割合の方が大きいでしょうし、購買行動に地域による差がそれほど認められるとは、経験的にも思えません。もし仮に地域による異質性というものがそれほど強くないのであれば、むしろそういった情報を積極的に活かした方が良いのではないでしょうか?

目的

そのような時に、パラメータ間に緩やかな縛りをかける方法として階層ベイズモデルという手法があります。詳しくは書籍1を参考にして頂くとして、ここでは「各地域の広告効果を表すパラメータ \betaを生成する分布を仮定する」方法と説明しておきます。この分布の幅(バラつき、分散)が十分に大きければ広告効果は地域によって様々な値を取り得ますし、逆に分布の幅が極端に狭ければ地域ごとの広告効果に差はないことを意味します。

この階層構造を持つモデルを、時系列データの分析でお馴染みの状態空間モデルと合わせて取り扱いたいというのが今回の試みです。具体的には、以下のようなモデルを想定しています。

 {
\begin{equation}
y_{A,t} = \mu_{A,t} + \beta_{A}X_{A,t} + \epsilon_{A,t} \
\end{equation}
}

 {
\begin{equation}
\mu_{A,t} \sim N(\mu_{A,t-1}, \sigma^{2}_{\mu}) \
\end{equation}
}

 {
\begin{equation}
\epsilon_{A,t} \sim N(0, \sigma^{2}_{\epsilon}) \
\end{equation}
}

 {
\begin{equation}
\beta_{A} \sim N(\beta_0, \sigma^{2}_{\beta})
\end{equation}
}

添字のAおよびtはそれぞれ地域(Area)と時点(time)を意味しており、 y_{A,t}は「ある地域At時点における観測値」となります。 ポイントは4行目で、地域ごとの広告効果 \beta_{A}に対して、その生成元となる分布(ここでは正規分布)を仮定しています。もし \sigma_{\beta}^{2}が大きければ \betaは様々な値を取り得ますが、値が小さければ非常に似通った値になることが理解できると思います。

ライブラリの読み込み

今回の分析では{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)

f:id:ushi-goroshi:20180220113808p:plain

また、この時点で地域をまとめて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)

f:id:ushi-goroshi:20180220114352p:plain

おぉ、かなり精度良く推定できているようですね!状態変数の初期値や分散の推定値はどうでしょうか?

> 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()

f:id:ushi-goroshi:20180220114731p:plain

こちらも、事前に設定した回帰係数と推定値が似通っているようです。数値を確認しても、良い精度で推定できています。

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効果の推定を組み込めば試してみたかったモデルは一通り完了することになります。

これらについても追って紹介したいと思います。


  1. 伊庭先生のベイズモデリングの世界、久保先生の緑本、松浦さんのアヒル本など良書はたくさんあります。

スタイン推定量を確かめてみる

いきなりですが、最近購入したベイズモデリングの世界を読んでいて非常に面白い話題を発見しました。

ベイズモデリングの世界

ベイズモデリングの世界

P118から始まるスタイン推定量についてなのですが、その直感的な意味について本文から引用すると:

それぞれの y_iを、全部の{ y_i}の平均値 \frac{1}{n}\Sigma_{i=1}^{n}y_iの方向に aだけ引っ張ってやる

ことで、 y_iそのものを推定量とするよりも、パラメータの真値 \hat{\theta_{i}}との二乗誤差の期待値が小さくなるそうです。これはいわゆる縮小推定量の一種で、引っ張りの程度 aをデータから適応的に求めるのがこの推定量のキモなのだそうです。

スタイン推定量 \hat{\theta_{i}^{S}}の定義自体はそれほど難しいものでもなく、以下のような式によって表されます:

 \displaystyle \hat{\theta_{i}^{S}} = (1-a)y_{i} + \frac{n}{a}\Sigma_{i=1}^{n}y_i

ここで、

 \displaystyle a = \frac{\sigma^{2}}{s^{2}}, \
s^{2} = \frac{1}{n-3}\Sigma_{i=1}^{n}(y_i - \frac{1}{n}\Sigma_{j=1}^{n}y_i)^{2}

です。この記事では、スタイン推定量が本当に y_iよりも二乗誤差が小さくなるのか、以下のように実験してみます。

実験内容

内容としては非常に簡単なもので、平均と分散が共に既知である正規分布からいくつかのサンプルを生成します。そのサンプルをそのまま用いた場合と、スタイン推定量を用いた場合とで、真のパラメータとの二乗誤差の期待値:

 \mathbb{E}[\Sigma_{i=1}^{n}(\theta_{i}^{*}({y_i}) - \theta_i)^{2}]

にどれほど差が生じるかを確認しています。

パラメータ設定

実験では、一回の試行につき平均が{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

書籍にはスタイン推定量が誤差を減少させる証明も付いていますが、まだ読み込めていませんので上記の理解が間違っているかもしれません。しかしこの推定量自体も非常に面白く、紹介したいと思ったので実験してみた次第です。


  1. なお「ベイズモデリングの世界」ではこれをバイアス-バリアンスのトレードオフとして解説しています。

Bayesian GLMで過適合を避ける

長らく知人に貸していた「みんなのR」が手元に戻ってきたのでパラパラと見ていたら、見逃していた面白いトピックを見つけたので紹介します。内容としてはこちらの記事とほぼ同じで、基本的には書籍の該当部分を再現しているだけですが、{purrr}のmapを使った書き方にしてみたり、少しだけ変更を加えています。

書籍はこちら。P320の「Bayesian 縮小」からです。

みんなのR -データ分析と統計解析の新しい教科書-

みんなのR -データ分析と統計解析の新しい教科書-

分析環境の準備

ライブラリの読み込み

まずは必要なライブラリの準備から。以下のライブラリを使用します。

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の各要素にlmglmの結果を持つオブジェクトを渡し、関心のある変数の回帰係数を比較することができます。

いずれも同じ結果となっています。

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))

f:id:ushi-goroshi:20180115164702p:plain
(同じグラフなので二枚は省略)

さてここで結果を見てみると、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()

f:id:ushi-goroshi:20180115165807p:plain

係数とその誤差がかなり改善されていますね!

事前分布の仮定が他の年の係数に影響を与えてはいないでしょうか?それも確認してみましょう。

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))

f:id:ushi-goroshi:20180115170220p:plain

> 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]))

f:id:ushi-goroshi:20180115170528p:plain

おまけ

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)))

f:id:ushi-goroshi:20180115170829p:plain

ほとんど一致していますね!

終わりに

この記事では、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)

f:id:ushi-goroshi:20180111224629p:plain

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)

f:id:ushi-goroshi:20180111225340p:plain

フィッティング

この構造を用いてパラメータの推定を行おうとするとエラーが出ます。どうやらグラフの構造が悪く、フィッティングできないようです。

> 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間のエッジが無向になっていることがわかります。

f:id:ushi-goroshi:20180111225340p:plain


そこでこの構造に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)

f:id:ushi-goroshi:20180111230018p:plain

では、この構造を用いて再度フィッティングしてみましょう。

> 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を追加

f:id:ushi-goroshi:20180111230351p:plain

フィッティング

続いてフィッティングです。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)

f:id:ushi-goroshi:20180111231012p:plain

どうもgsを用いた時と比べて構造がかなり異なるようですね。比較してみましょう。

par(mfrow = c(1, 2))
plot(str_02)
plot(str_05)

f:id:ushi-goroshi:20180111231120p:plain

"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)

f:id:ushi-goroshi:20180111231456p:plain

何と、標準化すると構造が変わってしまいました。。。フィッティングは省略します。

分析実行その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)

f:id:ushi-goroshi:20180111231826p:plain
f:id:ushi-goroshi:20180111231837p:plain
f:id:ushi-goroshi:20180111231930p:plain
f:id:ushi-goroshi:20180111231942p:plain

アルゴリズムによって少しずつ推定された構造が違うのですが、low -> smokeおよびrace -> smokeのエッジは概ね共通して見られるようです。しかし残念ながらこのデータは「出生児の低体重の有無」と「母体環境」であるため、エッジの向きは逆ですね。。。

終わりに

この記事ではベイジアンネットワークを扱うためのパッケージ およびその使い方を示しました。結果的にデータのみから因果関係を推定するのは困難であり、良い構造推定のためにはドメイン知識を持ったメンバーとの協業が必要となりそうですが、こういったデータの構造を確認しながらのディスカッションは非常に楽しそうですね。

WAICを計算してみる

ある統計モデルの予測精度を表す指標としてAICAkaike Information Criterion)というものがよく使われますが、stanを用いてモデリングした場合にはAICを計算することができません。しかし2009年に渡辺澄夫先生が発表したWAICWidely Applicable Information Criterion、広く使える情報量規準)はAICをより一般化したものとして定義され、MCMCによる事後分布から計算することができます。

現在のところ{rstan}では直接WAICを計算することができませんが、{loo}というパッケージを用いることでWAICを得ることができます。この記事では、stanの結果からWAICを計算により求め、それがlooの結果と一致することを確認します。

WAICの定義

まず始めにWAICの定義を確認しましょう*1。WAICは以下のように定義されます:

WAIC = -2(lppd - p_{WAIC})

ここで

\displaystyle lppd = \sum_{i = 1}^{N}\log{Pr(y_{i})}

であり、また

\displaystyle p_{WAIC} = \sum_{i = 1}^{N}V(y_{i})

です。Pr(y_{i})およびV(y_{i})は、それぞれ観測値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は前述の通り以下によって計算されます:

\displaystyle lppd = \sum_{i = 1}^{N}\log{Pr(y_{i})}

まずは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サンプルのそれぞれをパラメータとして与えた時の各 y_{i}の対数尤度を求めると:

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)))

で得られます。
同じように p_{waic}を求めてみましょう。 p_{waic}の定義は以下のようでした:

\displaystyle p_{WAIC} = \sum_{i = 1}^{N}V(y_{i})

ここから、

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という本を読んでいます。

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

その中でタイトルの通り面白い話があったので紹介がてら実行してみます。以下は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には他にも多重共線性が生じている時の挙動など、関心を引く話題が多く散りばめられており、非常に読み応えのある本なので、しばらくこれで勉強しようと思います。

*1:というかそうなるようにseedを色々と試しました

*2:ファイル名やフォルダは適当です

*3:ちなみに、最近私のMacではXcodeが更新された際に使用許諾を同意していなかったためコンパイルが通らず、しばらく悩みました

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法」を指定することで矩形制約を与えることができる