機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan：階層モデルに関する考察

Stan

はじめに


今回もStanを使って階層モデルを勉強していく。

今回は階層モデル（というよりベイズ学習そのもの）の基本的なことに関する疑問の解消を目指す。

最尤推定ではなく、ベイズ推論を用いる目的の１つとして「ベイズ推論による機械学習」では以下のように書かれている。

ベイズ推論の事後分布は、Nを大きくしていけば漸近的に最尤推定の結果に近づきます。
しかし、「データ数Nが十分に大きい」という仮定自体が現実問題にアプローチするうえで適切ではありません。
（中略）
データ数Nが十分であると思うのであれば、解析対象をもっと詳細にするべきです。

下記でも同様なことが書かれている。
andrewgelman.com

以下はベイズ統計の扱い方に関する基礎的な論文（メモ）
Philosophy and the practice of Bayesian statistics

「なるほど、そらそうだ」と思ったわけだが、勉強し始めた当時はここである疑問が生じた。
詳細に解析するってことはデータを分割することになるのだが、「その分割基準をカテゴリ変数として含めるのと何が違うのか」と思った。
少し分かりにくいので以下で例を示す。
①「年収とA党支持（0：不支持、1：支持）の関係」を回帰分析したい場合のデータは

x（年収） y（A党支持）
500 1
300 0
1200 1
850 1
250 0

をロジスティック回帰すると、Nが十分な場合はどの程度関係があるかが容易に分かる。
これはいわゆるビッグデータを用いた分析であり最尤法で問題ないとされる。
例えば、このデータから分かることは以下である。

  • 年収が高くなれば、A党を支持しやすくなる。

②次に以下を考える。
先ほどのデータを取得した際、「性別」（0：男性、1：女性）データも同時に取得しており、データは

x1（年収） x2（性別） y（A党支持）
500 0 1
400 1 0
300 0 0
500 1 0
1200 0 1
850 0 1
250 1 0

であったとする。
これもNが十分な場合は「性別」をダミー変数として、ロジスティック回帰すると容易に分析できる。
このデータであればから、例えば以下の情報を得れる。

  • 年収が高くなれば、A党を支持しやすくなる。
  • 男性の方が、A党を支持しやすい。

②では、データを分割せずに①よりも詳細な情報が得られている。
これは単に説明変数を増やしているだけであり、ベイズ推論とは関係がない。

試してみる


まず、以下のようなデータを用意。

data <- read.csv("bayes_logit.csv")
head(data)
          x1 x2 y
1  0.5850919  1 0
2 -3.0039250  0 1
3  1.1599543  1 0
4  2.6738048  0 0
5  0.8788801  1 0
6 -0.4721039  0 1

※テストデータはGithubより取得
repo/data at master · kento1109/repo · GitHub

GLM

一般線形モデルで比べてみる。

①年収のみで解析
モデル式は以下の通り。
P(y_i=1)=logit^{-1}(\beta_0 + \beta_1x_1)

summary(fit)

Call:
glm(formula = y ~ x1, family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7917  -0.7198   0.0035   0.7690   2.7342  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.09074    0.05770  -1.573    0.116    
x1          -1.23773    0.05565 -22.242   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2772.6  on 1999  degrees of freedom
Residual deviance: 1853.1  on 1998  degrees of freedom
AIC: 1857.1

Number of Fisher Scoring iterations: 5p

係数が有意であることが分かる。

分類器の指標としてAUCを評価する。

pred_y <- predict(fit,data,type="response")

library("ROCR")
pred <- prediction(pred_y, data$y)
performance(pred,"auc")@y.values
[[1]]
[1] 0.860869

※今回は、訓練データ自体で精度評価

②年収＋性別で解析
モデル式は以下の通り。
P(y_i=1)=logit^{-1}(\beta_0 + \beta_1x_1 + \beta_2x_2)

fit <- glm(y~.,data,family="binomial")
summary(fit)

Call:
glm(formula = y ~ ., family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9677  -0.3726   0.0100   0.4470   3.6465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.39849    0.10171   13.75   <2e-16 ***
x1          -1.58768    0.07566  -20.98   <2e-16 ***
x2          -3.39827    0.17427  -19.50   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2772.6  on 1999  degrees of freedom
Residual deviance: 1261.4  on 1997  degrees of freedom
AIC: 1267.4

Number of Fisher Scoring iterations: 6

両係数とも有意であることが分かる。

pred_y <- predict(fit,data,type="response")

pred <- prediction(pred_y, data$y)
performance(pred,"auc")@y.values
[[1]]
[1] 0.938541

AUCは②の方が高く、性別は政党の決定に重要な変数であることが分かった。

Stanで階層モデル

①のモデル式
P(y_i=1)=inv\_logit(\beta_0 + \beta_1x_1)
の問題点は「性別差」が考慮されていないことであった。
すなわち、

  • 一郎さん（年収500万の男性）
  • 花子さん（年収500万の女性）

のどっちが、A党を支持しやすいのか分からないということである。
これではより詳細な分析ができないのでもったいない。

さて、今回は階層モデルにおけるアプローチとして、「性別における係数」をモデル式に追加する。
モデル式は以下の通り。
P(y_i=1)=logit^{-1}(\beta_0 + \beta_1x_1 + \beta_s[\rm{KID}])
SIDは各データの性別インデックスを表す。
この係数が「政党支持にかんする男女差」を表している。

Stanコードは以下の通り。

data {
  int N;
  int K;
  real X[N];
  int<lower=0,upper=1> y[N];
  int<lower=1,upper=K> KID[N];
}
parameters {
  real beta[2];
  real b_s[K];
  real<lower=0> s_beta;
}

transformed parameters{
  real z[N];
  for (n in 1:N)
    z[n] = inv_logit(beta[1] + X[n]*beta[2] + b_s[KID[n]]);
}

model {
  s_beta ~ cauchy(0,5);
  for (k in 1:K){
    b_s[k] ~ normal(0, s_beta);
  }
  for (n in 1:N) 
    y[n] ~ bernoulli(z[n]);
}

少し嵌ったがのが、性別差の事前分布の分散s_beta。
元々、無情報事前分布（一様分布）にしていたが、うまく収束しなかった。
階層モデルの分散パラメータの事前分布に関しては、以下を目安とすると良いらしい。

  • グループ数>5 の場合は、一様分布で問題ない。
  • それ以下の場合、半コーシー分布が良い。

詳しくはこちらを参照

www.slideshare.net

実行

以下の通り、Rでデータを作って実行した。

data <- read.csv("bayes_logit.csv")
dat <- list(N=nrow(data),K=2,X=data$x1,y=data$y,KID=data$x2+1)
stan.fit <- stan(file = "logistic.stan", data=dat)
推定結果

サンプリングは5分ほどで終わったが、以下のような警告が表示された。

Warning messages:
1: There were 6 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
2: Examine the pairs() plot to diagnose sampling problems

Brief Guide to Stan’s Warnings
によると、maximum treedepthに達したパスがあった場合に表示されるエラーらしい。
とりあえず、言われるままにpairs()で結果を可視化してみた。

pairs(stan.fit_, pars = c("b_s[1]", "b_s[2]", "lp__"), las = 1)

f:id:kento1109:20180620121917p:plain

maximum treedepthに達したパスは黄色くプロットされるらしい。
見えにくいが、数か所黄色い点が存在する。
これ自体は警告であり、開発チームも深刻なものではないと言っているので、いったん無視する。

推定結果自体は問題なさそうである。

summary(stan.fit)
$summary
                 mean      se_mean           sd          2.5%           25%
beta[1] -6.038971e-01 3.957173e-01 6.117218e+00 -1.205801e+01 -1.986270e+00
beta[2] -1.594435e+00 1.873761e-03 7.506944e-02 -1.747818e+00 -1.643088e+00
b_s[1]   2.004916e+00 3.955603e-01 6.117507e+00 -7.941019e+00 -5.643430e-02
b_s[2]  -1.401408e+00 3.959316e-01 6.119758e+00 -1.137238e+01 -3.450435e+00

                  50%           75%         97.5%     n_eff     Rhat
beta[1] -2.166488e-01  1.453795e+00  9.395985e+00  238.9669 1.015113
beta[2] -1.594364e+00 -1.544348e+00 -1.447651e+00 1605.0841 1.001643
b_s[1]   1.622581e+00  3.381019e+00  1.346151e+01  239.1793 1.015064
b_s[2]  -1.771588e+00 -5.230465e-02  1.004941e+01  238.9066 1.015121

plot(stan.fit, pars = c("beta[1]", "beta[2]", "b_s[1]", "b_s[2]"))

f:id:kento1109:20180620122915p:plain:w500

ここで、男女の係数に関して注目する。

ms <- rstan::extract(stan.fit)
mean(ms$b_s[,1])
[1] 2.004916
mean(ms$beta[,2])
[1] -1.594435

男女に平均に違いがあることが分かる。
一応、t検定でも確認する。

t.test(ms$b_s[,1],ms$b_s[,2],var.equal = T)

	Two Sample t-test

data:  ms$b_s[, 1] and ms$b_s[, 2]
t = 24.897, df = 7998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 3.138127 3.674521
sample estimates:
mean of x mean of y 
 2.004916 -1.401408

→有意に平均に差があるという結果

男性の場合の係数が2.0、女性の場合の係数が-1.6なので、
「男性の場合はA政党を支持しやすく、女性の場合はA政党を支持しにくい」ことが分かる。
その差は3.6なので、オッズでいうと36.6（exp(3.6)）ということが分かる。
※②（年収＋性別）でGLMで解析した時のオッズとも近い推定値であることも分かる。

さいごに


上記の例では男女差を階層モデルに組み込んでベイズ推論を行った。
結果は、男女をダミー変数としてGLMするのと似たような解析となった。
男女差のような取り得る値が２値と分かっている場合、それを階層モデルで推定する意味はあんまりない。
（データ数が十分にあれば、GLMで全く問題ないと思う。）
そうではなく、もっとグループが多い変数（会社差、個人差など）を扱おうとすると、当然各グループのデータ数は少なくなってくる。そうなると、ダミー変数を用意したGLMの推定が難しくなってくる。
そんなときに役に立つのがその差を組み込んだ階層モデルなのだと思う。