R
Tableau
Stan

アヒル本(StanとRでベイズ統計モデリング)5.3節 ロジスティック回帰をやってみた(RStan)

記事の目的

アヒル本5章の3節 ロジスティック回帰をRとTableauの自習として挑戦したものです
アヒル本とは『StanとRでベイズ統計モデリング』のことです

アヒル本ありきの内容で目新しい情報はありません
5.3節 ロジスティック回帰は5.1節や5.2節ですでに説明済みの箇所は端折っていると感じたので挑戦しただけです

記事のまとめ

  • アヒル本の5.3節を著者の推奨する分析のステップに従って実施した
  • 5.3節のモデルより、Chapter 5の練習問題4のモデルのほうがあてはまりがよいことがわかった

著者の推奨する分析のステップ

著者の推奨する分析のステップは、 「Chapter3 統計モデリングを始める前に」の章で説明

# 手順 概要
1 解析の目的 データで何を知りたいんか、何を主張したいのか
2 データの分布の確認 ヒストグラムの作成、散布図の作成、クロス集計
3 メカニズムの想像 データを生成するメカニズム、データをつなぐメカニズムを考える
4 モデル式の記述 メカニズムをモデル式に落とし込む
5 Rでシミュレーション シミュレーションでモデル式がどのような形式か確認する
6 Stanで実装 Stanのコードを書く、パラメータの推定を実行する
7 推定結果の解釈 推定結果やベイズ信頼区間をもとに解釈したり
8 図によりモデルのチェック モデルがうまく当てはまっていそうか図でチェックする

アヒル本5章 ロジスティック分布の例題を分析する

ステップ1.解析の目的

「解析の目的」は5.3.1に記載されている

  • T大学のある三年生が授業に出席するかを予測するモデルをつくる
  • 解析した結果のモデルの評価を発表する必要がある

このモデルによって利用することで欠席しそうや生徒を予測し、出席を促すことができる可能性がある
モデルの評価は、欠席しそうや生徒を予測できるかとする(TBD)

ステップ2.データの分布の確認

「データの分布の確認」は5.3.2に記載されているに足りない部分があったので追加した
※足りない部分はChapter3の練習問題(3)の課題となっている

データの内容

# PersonID A Score Weather Y
1 1 0 69 B 1
2 1 0 69 A 1
3 1 0 69 C 1
4 1 0 69 A 1
5 1 0 69 B 1
6 1 0 69 B 1
... ... ... ... ... ...
2391 50 1 99 C 0
2392 50 1 99 B 1
2393 50 1 99 A 1
2394 50 1 99 A 1
2395 50 1 99 C 1
2396 50 1 99 C 0
d3 <- read.csv(file='https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap05/input/data-attendance-3.txt')
head(d3)
  PersonID A Score Weather Y
1        1 0    69       B 1
2        1 0    69       A 1
3        1 0    69       C 1
4        1 0    69       A 1
5        1 0    69       B 1
6        1 0    69       B 1

変数

変数名 説明 とる値
PersonID 学生のID
A アルバイトがすきか 0:好きではない、1:好き
Score 学問への興味の強さを数値化したもの 200点満点
Weather 天気 a:晴れ、b:曇り、c:雨
Y 授業に出席したかどうか 0:欠席、1:出席

分布の確認

カテゴリ値や二値の場合は集計でデータの分布を確認する

天気と出席・欠席の集計

カテゴリ 欠席 出席
晴れ 306 953
曇り 230 500
138 269

アルバイトが好きかと出席・欠席の集計

カテゴリ 欠席 出席
アルバイトが好き 288 994
アルバイトが好きではない 386 728
aggregate(Y~ Weather,data=d3,FUN=table)
  Weather Y.0 Y.1
1       A 306 953
2       B 230 500
3       C 138 269
aggregate(Y~ A,data=d3,FUN=table)
  A Y.0 Y.1
1 0 288 994
2 1 386 728

ステップ3.メカニズムの想像

「メカニズムの想像」は5.3.4に記載されている。解釈を追記した

  • 天気(Weather)が悪くなると出席率が悪くなる
  • アルバイト(A)が好きだと出席率が悪くなる
  • 学問の興味が強い(Score)と出席率が良くなる
  • 予測したい変数は、0,1の二値である
  1. Weather,A,Scoreを予測変数としロジスティック関数で出席確率に変換する
  2. 出席確率をベルヌーイ分布からY(出席か欠席)を決める

ステップ4.モデル式の記述

「モデル式の記述」は5.3.5と、Chapter3の練習問題(4)の課題がある

  1. 5.3.5のモデルはWeatherの値を数値に変換している
  2. Chapter5の練習問題(4)のモデルはWeatherの値をダミー変数に変換している

Weatherの値を数値化したモデル

Weatherは、離散的な文字列データなので、数値に変換する

5-5
q[i] = inv\_logit(b_1 + b_2A[i] + b_3Score[i] + b_4Score Weather[i] )   i = 1,...,I

\\

Y[i] = Bernoulli(q[i])   i = 1,...,I

Weatherの値と数値の対応

カテゴリ名 文字列 数値
晴れ A 0
曇り B 0.2
C 1

Weatherの値をカテゴリ化したモデル

Weatherの晴れを0として、曇りまたは雨の場合に効果を加える

5-5
q[i] = inv\_logit(b_1 + b_2A[i] + b_3Score[i] + bw_2Score Weather_Cloudy[i] + bw_3Score Weather_Rain[i]  )   i = 1,...,I

\\
 bw_2 : 曇りの効果 / Weather_Cloudy:曇りを表すダミー変数\\
 bw_3 :雨の効果 / Weather_Rain :雨を表すダミー変数 \\
\\

Y[i] = Bernoulli(q[i])   i = 1,...,I

ステップ5.Rでシミュレーション

4.4.4を参考にした

ベルヌーイ分布を使ったglmを見てみる

Weatherの値を数値化したモデル

Weatherを数値に変換したパラメータWにする

glm_numeric.R
conv <- c(0, 0.2, 1)
names(conv) <- c('A', 'B', 'C')
data <- list(I=nrow(d3), A=d3$A, Score=d3$Score/200, W=conv[d3$Weather], Y=d3$Y)

f <- glm(Y~A+Score+W, family=binomial,data=data)

予測変数のp値がある程度小さいので問題なさそう

summary_numeric.R
summary(f)

Call:
glm(formula = Y ~ A + Score + W, family = binomial, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0533  -1.3475   0.7020   0.8415   1.2578  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.17108    0.22838   0.749 0.453792    
A           -0.61999    0.09281  -6.681 2.38e-11 ***
Score        1.95408    0.36527   5.350 8.81e-08 ***
W           -0.46106    0.12362  -3.730 0.000192 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2847.3  on 2395  degrees of freedom
Residual deviance: 2762.1  on 2392  degrees of freedom
AIC: 2770.1

Number of Fisher Scoring iterations: 4

ROC曲線とAUCの確認

ROC曲線をみて分類性能を確認する

roc_numeric.R
prob=predict(f,type=c("response"))
data$prob=prob
install.packages("pROC")
library(pROC)
roc(Y ~ prob, data = data)

Call:
roc.formula(formula = Y ~ prob, data = data)

Data: prob in 674 controls (Y 0) < 1722 cases (Y 1).
Area under the curve: 0.61

rocの結果をplot

AICをみると0.61なので、分類性能はよくない(0.8位が目標とのこと)

Weatherの値をカテゴリ化したモデル

glm関数はカテゴリカルデータを自動的にダミー変数に変換してくれるらしいので、Weatherをそのまま渡す

glm_summary.R
summary(glm(Y~A+Score+Weather, family=binomial,data=d4))

Call:
glm(formula = Y ~ A + Score + Weather, family = binomial, data = d4)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0966  -1.3038   0.6916   0.8408   1.2322  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.265917   0.231294   1.150 0.250270    
A           -0.624011   0.092980  -6.711 1.93e-11 ***
Score        0.009807   0.001829   5.362 8.22e-08 ***
WeatherB    -0.379258   0.104972  -3.613 0.000303 ***
WeatherC    -0.495455   0.125793  -3.939 8.19e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2847.3  on 2395  degrees of freedom
Residual deviance: 2754.2  on 2391  degrees of freedom
AIC: 2764.2

Number of Fisher Scoring iterations: 4

ROC曲線とAUCの確認

ROC曲線をみて分類性能を確認する

category_roc.R
f4<-glm(Y~A+Score+Weather, family=binomial,data=d4)
d4$prob=predict(glm(Y~A+Score+Weather, family=binomial,data=d4),type=c("response"))
roc(Y ~ prob, data = d4)
Call:
roc.formula(formula = Y ~ prob, data = d4)

Data: prob in 674 controls (Y 0) < 1722 cases (Y 1).
Area under the curve: 0.6241

roc

AUCは0.6241
Weatherの値を数値化したモデルの0.61よりよいことがわかる

AICの比較

説明変数をかえてAICもみてみる
AICが小さいモデルがうまく当てはまっている
AICが一番良いモデルは”全ての説明変数を使ったWeatherの値をカテゴリ化したモデル”となる

カテゴリ名 AIC
Y ~ A + Score + W 2770.1 Weatherの値を数値化したモデル
Y ~ A + Score + Weather 2764.2 Weatherの値をカテゴリ化したモデル
Y ~ A + Score 2781.7 Weatherを組み込まないモデル
Y ~ Score + W 2797.2 Aを組み込まない、Weatherの値を数値化したモデル
Y ~ A + W 2813.3 Scoreを組み込まない、Weatherの値を数値化したモデル

ステップ6.Stanで実装

『StanとRでベイズ統計モデリング』のサポートページを参考にした

Weatherの値を数値化したモデル

Rの出力結果 

q[I]は、個体値なので省略

fit
Inference for Stan model: model5-5.
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% n_eff Rhat
b[1]        0.17    0.01 0.22    -0.26     0.02     0.17     0.33     0.63  1995    1
b[2]       -0.62    0.00 0.09    -0.81    -0.69    -0.62    -0.56    -0.44  2354    1
b[3]        1.95    0.01 0.36     1.26     1.71     1.95     2.19     2.66  1982    1
b[4]       -0.46    0.00 0.12    -0.71    -0.55    -0.47    -0.38    -0.22  2728    1
q[1]        0.68    0.00 0.02     0.63     0.66     0.68     0.70     0.73  2105    1
...略

収束の確認

結果からパラメータの要約を見る
以下を満たしていることから収束しているとみなす

  • n_eff が大きい
  • Rhatは1.1より小さい
パラメータ 説明 mean n_eff Rhat
b[1] 切片 0.17 1995 1
b[2] A -0.62 2354 1
b[3] Score 1.95 1982 1
b[4] W(天気を数値化したもの) -0.46 2728 1

ggmcmcを使い収束診断ファイルの出力する
パラメータbのtraceplotのみを出力するRスクリプト
lp__も出力したほうがいいかも

ggmcmc(
  ggs(fit,inc_warmup=TRUE,family="b",stan_include_auxiliar=FALSE),
  file='fit-traceplot.pdf',
  plot='traceplot')

図で見ても収束していることがわかる
fit-traceplot.pdf

glmと比較

glmとMCMCサンプルの推定結果に違いはなさそう

パラメータ 説明 MCMCサンプルのmean glmの結果  
b[1] 切片 0.17 0.17108
b[2] A -0.62 -0.61999
b[3] Score 1.95 1.95408
b[4] W -0.46 -0.46106

Weatherの値をカテゴリ化したモデル

カテゴリカルデータ(A,B,C)を自動的にダミー変数(1,2,3)に変換している

Rの出力結果 

省略

収束の確認

結果からパラメータの要約を見る
以下の値から収束しているとみなす

  • n_eff が大きい
  • Rhatは1.1より小さい
パラメータ 説明 mean n_eff Rhat
b[1] 切片 0.26 1961 1
b[2] A -0.63 3183 1
b[3] Score 1.97 1934 1
bw2 曇りの効果を表すダミー変数 -0.38 2905 1
bw2 雨の効果を表すダミー変数 -0.49 2906 1

ggmcmcを使い収束診断ファイルの出力する
lp__も出力したほうがいいかも

ggmcmc(
  ggs(fit,inc_warmup=TRUE,stan_include_auxiliar=FALSE),
  file='fit-traceplot2.pdf',
  plot='traceplot')

図で見ても収束していることがわかる
fit-traceplot2.pdf

glmと比較

Step.7 推定結果の解釈

「推定結果の解釈」は5.2.6にある二項ロジスティック回帰をベースに5.1.7も参考にした

Weatherの値を数値化したモデル

回帰係数の解釈

パラメータ 説明 mean 2.5% 50% 97.5%
b[1] 切片 0.17 -0.26 0.17 0.63
b[2] A -0.62 -0.8 -0.62 -0.44
b[3] Score 1.95 1.26 1.95 2.66
b[4] W(天気を数値化したもの) -0.46 -0.70 -0.47 -0.21

モデル式5-5に推定結果を当てはめる

5-5
q[i] = inv\_logit( 0.17 -0.62A[i] + 1.95Score[i]  -0.46Score Weather[i] )   i = 1,...,I
\\

Y[i] = Bernoulli(q[i])   i = 1,...,I

b[2]とb[4]は負の効果、 b[3]は正の効果があることがわかるが、線形回帰じゃないので解釈できるのは以下くらい 

  • アルバイトが好きな生徒はアルバイトが好きな生徒にくらべて出席率が悪い傾向
  • 天気(W)が良いとが出席率が悪くなる傾向
  • 学問の興味が強い(Score)と出席率が良くなる傾向

パラメータ幅の確認

線形回帰じゃないので解釈できない可能性はあるが、どのパラメータの95%信頼区間に0を含んでいないので、良いまたは悪い傾向に有意性はありそう

オッズ比

MCMCサンプルの平均とglmの結果に違いがないのでオッズ比はglmの結果から求めた
晴れの時のオッズと雨の時のオッズを比較すると、オッズ比は0.53になる

exp(f$coefficients)
(Intercept)           A       Score           W 
  1.1865853   0.5379483   7.0574164   0.6306160 

Weatherの値をカテゴリ化したモデル

回帰係数の解釈

パラメータ 説明 mean 2.5% 50% 97.5%
b[1] 切片 0.26 -0.20 0.26 0.71
b[2] A -0.63 -0.80 -0.62 -0.45
b[3] Score 1.97 1.28 1.97 2.68
bw2 曇りの効果を表すダミー変数 -0.37 -0.588 -0.37 -0.17
bw3 雨の効果を表すダミー変数 -0.49 -0.74 -0.49 -0.24

モデル式5-5に推定結果を当てはめる

5-5
q[i] = inv\_logit( 0.17 -0.62A[i] + 1.95Score[i] +  0.37 Weather_Cloudy[i] -0.49  Weather_Rain[i] )   i = 1,...,I
\\

Y[i] = Bernoulli(q[i])   i = 1,...,I

b[2]とb[4]は負の効果、 b[3]は正の効果があることがわかるが、線形回帰じゃないので解釈できるのは以下くらい 

  • アルバイトが好きな生徒はアルバイトが好きな生徒にくらべて出席率が悪い傾向
  • 天気(W)が良いとが出席率が悪くなる傾向
  • 学問の興味が強い(Score)と出席率が良くなる傾向

パラメータ幅の確認

省略

オッズ比

省略

Step8. 図によるモデルのチェック

「図によるモデルのチェック」は 5.3.7にある、可視化は Tableauで作成する

応答変数とqの信頼区間を重ねる

TBD

確率と実測値のプロット

TBD

ROC曲線

TBD