ベイズ統計モデリング備忘録
はじめに
今回は順序ロジスティック回帰をやってみる。
順序ロジスティック回帰分析
使用するパッケージの読み込み
library(tidyverse)## ─ Attaching packages ────────────────── tidyverse 1.2.1 ─## ✔ ggplot2 3.2.0 ✔ purrr 0.3.3
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0## ─ Conflicts ──────────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(rstan)## Loading required package: StanHeaders## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extractlibrary(ggmcmc)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(ordinal)##
## Attaching package: 'ordinal'## The following object is masked from 'package:dplyr':
##
## sliceデータ作成
set.seed(1234)
n <- 1000
X1 <- rnorm(n = n, mean = 50, sd = 10)
X2 <- runif(n = n, min = 0, max = 1) %>% round(0)
X3 <- runif(n = n, min = 1, max = 4) %>% round(0)
Y <- 0.1 * X1 + 0.4 * X2 + 0.2 * X3 + rnorm(n = n, mean = 0, sd = 2)
threshold1 <- quantile(Y, 0.3) %>% as.numeric()
threshold2 <- quantile(Y, 0.6) %>% as.numeric()
Y <- ifelse(Y < threshold1, 1, ifelse(Y < threshold2, 2, 3)) %>% as.factor()
df <- tibble(X1 = X1, X2 = X2, X3 = X3, Y = Y)
print(df)## # A tibble: 1,000 x 4
## X1 X2 X3 Y
## <dbl> <dbl> <dbl> <fct>
## 1 37.9 0 4 1
## 2 52.8 0 2 2
## 3 60.8 1 1 3
## 4 26.5 1 2 2
## 5 54.3 0 2 1
## 6 55.1 1 2 1
## 7 44.3 1 2 1
## 8 44.5 0 3 2
## 9 44.4 1 2 1
## 10 41.1 0 4 1
## # … with 990 more rows通常の順序ロジスティク回帰分析
res <- clm(Y ~ X1 + X2 + X3, data = df)
summary(res)## formula: Y ~ X1 + X2 + X3
## data: df
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 1000 -1005.30 2020.60 4(0) 1.44e-09 2.0e+05
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## X1 0.081319 0.006829 11.907 <2e-16 ***
## X2 0.219841 0.121319 1.812 0.070 .
## X3 0.069365 0.063128 1.099 0.272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 3.3592 0.3750 8.958
## 2|3 4.7907 0.3908 12.260ベイズによるアプローチ
順序ロジットは従属変数\(y\)（1〜3とする）の背後にある何らかの潜在変数\(y^*\)を仮定してやり、それが何らかの閾値1より低ければ1、閾値1〜閾値2の間にあれば2、閾値2より大きければ3となっていると仮定する。そして、この潜在変数\(y^*\)を独立変数で予測してやる。なお、その時の誤差は標準化したロジスティック分布に従うとする。すなわち、以下のように表現できる。なお、ここでは簡単のため独立変数を1つにしている。 \[ y^* = beta_1 X + \epsilon \\ \begin{equation} \left\{ \begin{aligned} if \ y^*< threshold1 \ then \ y =1 \\ if \ threshold1 < y^* < threshold2 \ then \ y =2 \\ if \ threshold2 < y^* \ then \ y =3 \\ \end{aligned} \right. \end{equation} \] \[ \epsilon 〜 Logistic Distribution(\mu = 0, s = 1) \] よって、以下のように表現できる。 \[ y^* - (\beta_1X) 〜 Logistic Distribution(\mu = 0, s = 1) \] なお、\(y = 1\)、\(y = 2\)、\(y = 3\)となる確率はそれぞれ以下のように表現できる。 $$ \[\begin{eqnarray} Pr(y = 1) = \\ Pr(y^* < threshold1) = \\ Pr(\beta_1X + \epsilon < threshold1) = \\ Pr(\epsilon < threshold1 - (\beta_1X)) \\ Pr(y = 2) = \\ Pr(threshold1 < y^* < threshold2) = \\ Pr(threshold1 < \beta_1X + \epsilon < threshold2) = \\ Pr(threshold1 - \beta_1X < \epsilon < threshold2- \beta_1X) \\ Pr(y = 3) = \\ Pr(threshold2 < y^* ) = \\ Pr(threshold2 <\beta_1X + \epsilon) = \\ Pr(threshold2 - \beta_1X < \epsilon ) \\ \end{eqnarray}\] \[ ここで、$\epsilon$は標準化されたロジスティック分布に従うと仮定しているので、上記のそれぞれの確率は標準化されたロジスティック分布の累積分布関数を用いて以下のように表現できる。 \] Pr(y = 1) = - 0 \ Pr(y = 2) = - \ Pr(y = 3) = 1 - \[ これらについて、各個人の確率を用いて最尤法でパラメータを求めるのが順序ロジスティック回帰分析となる。上記の確率はロジスティック分布を用いながらその順序性を反映させた分布とも捉えることができる。これを仮に順序ロジスティック分布と呼ぶ。stanにはこのような順序ロジスティック分布を簡単に表現してくれる関数`ordered_logistic()`が用意されている。これは$c$（閾値）と$\eta$（上記の$\beta X$の部分）という2つのパラメータを有する。よって、これらのパラメータを推定してやることを考える。そこで、まず以下のように仮定する。 \] = X \[ 以上から従属変数$y$は以下のように順序ロジスティク分布($OrderedLogistic$)に従うと仮定できる。 \] Y〜OrderedLogistic(c, x) \[ なお、$k$（閾値の数）は従属変数$y$のカテゴリー数（$N_c$とする）から1を引いたものなので、以下のように考えられる。 \] k = N_c - 1 $$ 以上を踏まえてstan側のファイルに以下のように記述する。なお、これはordered_logit.stanとして保存しておく。
// ordered logistic regression
data {
int N;
int k;
real X1[N];
int<lower =0, upper = 1> X2[N];
int<lower =0, upper = 4> X3[N];
int Y[N];
}
parameters {
real beta1;
real beta2;
real beta3;
ordered[k - 1] c;
}
model {
beta1 ~ normal(0, 1000);
beta2 ~ normal(0, 1000);
beta3 ~ normal(0, 1000);
c ~ normal(0, 1000);
for(i in 1:N)
Y[i] ~ ordered_logistic(beta1 * X1[i] + beta2 * X2[i] + beta3 * X3[i], c);
}なお、求められるパラメータは\(\beta_1\)、\(\beta_2\)、\(\beta_3\)、\(c\)（\(k\)個）である。これらの事前分布は以下のように設定してある。 \[ \beta_1 〜 Nromal(0, 1000) \\ \beta_2 〜 Nromal(0, 1000) \\ \beta_3 〜 Nromal(0, 1000) \\ c 〜 Normal(0, 1000) \] それではMCMCサンプルを生成して結果を確認する。コードは以下の通りである。
data <- list(N = nrow(df), X1 = df$X1, X2 = df$X2, X3 = df$X3, Y = as.numeric(df$Y), k = 3)
fit <- stan(file = "ordered_logit.stan", data = data, seed = 1234)## Warning in readLines(file, warn = TRUE): '/Users/YONE/Desktop/
## ordered_logit.stan' で不完全な最終行が見つかりました##
## SAMPLING FOR MODEL 'ordered_logit' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001951 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 19.51 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 93.6716 seconds (Warm-up)
## Chain 1: 52.4389 seconds (Sampling)
## Chain 1: 146.111 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'ordered_logit' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001513 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 15.13 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 90.7078 seconds (Warm-up)
## Chain 2: 55.065 seconds (Sampling)
## Chain 2: 145.773 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'ordered_logit' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.00201 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 20.1 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 84.0923 seconds (Warm-up)
## Chain 3: 49.0269 seconds (Sampling)
## Chain 3: 133.119 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'ordered_logit' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.00161 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 16.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 83.3436 seconds (Warm-up)
## Chain 4: 48.1706 seconds (Sampling)
## Chain 4: 131.514 seconds (Total)
## Chain 4:print(fit)## Inference for Stan model: ordered_logit.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5%
## beta1 0.08 0.00 0.01 0.07 0.08 0.08 0.09 0.10
## beta2 0.22 0.00 0.12 -0.01 0.14 0.22 0.31 0.47
## beta3 0.07 0.00 0.06 -0.05 0.03 0.07 0.11 0.19
## c[1] 3.38 0.01 0.37 2.68 3.12 3.37 3.62 4.12
## c[2] 4.82 0.01 0.39 4.08 4.54 4.81 5.07 5.61
## lp__ -1007.47 0.04 1.59 -1011.51 -1008.29 -1007.17 -1006.31 -1005.37
## n_eff Rhat
## beta1 2182 1
## beta2 3230 1
## beta3 3002 1
## c[1] 1825 1
## c[2] 1770 1
## lp__ 1487 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 19 16:34:33 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).mcmc_sample <- rstan::extract(fit)traceplotでも収束を確認する。
g_mcmc <- fit %>%
ggs() %>%
ggs_traceplot()
print(g_mcmc)
最後にMCMCサンプルを用いて\(\beta_1\)〜\(\beta_3\)の事後分布を確認する。なお、緑の縦ラインが5パーセンタイル、95パーセンタイルである。赤い縦ラインが2.5パーセンタイル、97.5パーセンタイルである。また青い縦ラインが0.5パーセンタイル、99.5パーセンタイルである。黒の点線が0を意味する。
mcmc_df <- tibble(beta1 = mcmc_sample$beta1, beta2 = mcmc_sample$beta2,
beta3 = mcmc_sample$beta3)
get_quantile <- function(x){
x %>% quantile(c(0.05, 0.95, 0.025, 0.975, 0.005, 0.995)) %>% return()
}
print("beta1のパーセンタイル")## [1] "beta1のパーセンタイル"mcmc_df$beta1 %>% get_quantile() %>% print()## 5% 95% 2.5% 97.5% 0.5% 99.5%
## 0.07091238 0.09292705 0.06892764 0.09516128 0.06528206 0.09870584print("beta2のパーセンタイル")## [1] "beta2のパーセンタイル"mcmc_df$beta2 %>% get_quantile() %>% print()## 5% 95% 2.5% 97.5% 0.5% 99.5%
## 0.02154156 0.42590232 -0.01408207 0.46986018 -0.09379231 0.53982164print("beta3のパーセンタイル")## [1] "beta3のパーセンタイル"mcmc_df$beta3 %>% get_quantile() %>% print()## 5% 95% 2.5% 97.5% 0.5% 99.5%
## -0.03369016 0.17307138 -0.05324751 0.19313210 -0.09245136 0.22729446get_posterior_distribution <- function(data, x, xlab = "", colour = "blue"){
xlab = xlab
colour = colour
data %>%
ggplot() +
geom_density(aes(x = x), fill = colour, alpha = 0.3) +
geom_vline(xintercept = 0, colour = 'black', linetype = 'dashed') +
geom_vline(xintercept = quantile(x, 0.975), colour = 'red') +
geom_vline(xintercept = quantile(x, 0.025), colour = 'red') +
geom_vline(xintercept = quantile(x, 0.995), colour = 'blue') +
geom_vline(xintercept = quantile(x, 0.005), colour = 'blue') +
geom_vline(xintercept = quantile(x, 0.05), colour = 'green') +
geom_vline(xintercept = quantile(x, 0.95), colour = 'green') +
xlab(xlab) %>%
return()
}
g_beta1 <- mcmc_df %>%
get_posterior_distribution(mcmc_df$beta1, xlab = "beta1", colour = "blue")
g_beta2 <- mcmc_df %>%
get_posterior_distribution(mcmc_df$beta2, xlab = "beta2", colour = "red")
g_beta3 <- mcmc_df %>%
get_posterior_distribution(mcmc_df$beta3, xlab = "beta3", colour = "gray")
grid.arrange(g_beta1, g_beta2, g_beta3)