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

# 記事の目的
 
アヒル本５章の3節 ロジスティック回帰をRとTableauの自習として挑戦したものです
アヒル本とは[『StanとRでベイズ統計モデリング』](http://www.kyoritsu-pub.co.jp/bookdetail/9784320112421)のことです

アヒル本ありきの内容で目新しい情報はありません
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|図によりモデルのチェック|モデルがうまく当てはまっていそうか図でチェックする|

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

## ステップ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の値を数値に変換している
1. Chapter5の練習問題(4)のモデルはWeatherの値をダミー変数に変換している


### Weatherの値を数値化したモデル
Weatherは、離散的な文字列データなので、数値に変換する



```math: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の晴れを０として、曇りまたは雨の場合に効果を加える

```math: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にする

```splus: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値がある程度小さいので問題なさそう

```splus: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曲線をみて分類性能を確認する

```splus: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](https://qiita-image-store.s3.amazonaws.com/0/47064/e61428e4-3df0-e988-3ce3-1e79272d0a29.png)

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

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

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

```splus: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曲線をみて分類性能を確認する

```splus: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](https://qiita-image-store.s3.amazonaws.com/0/47064/a52a79d9-a4a0-536a-7079-93d8993a50cf.png)

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でベイズ統計モデリング』のサポートページ](https://github.com/MatsuuraKentaro/RStanBook/)を参考にした 

### Weatherの値を数値化したモデル
+ Stanコード
[model5-5.stan](https://github.com/MatsuuraKentaro/RStanBook/blob/master/chap05/model/model5-5.stan)
+ Rコード
[run-model5-5.R](https://github.com/MatsuuraKentaro/RStanBook/blob/master/chap05/run-model5-5.R)


#### 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')
```

図で見ても収束していることがわかる
<img width="413" alt="fit-traceplot.pdf" src="https://qiita-image-store.s3.amazonaws.com/0/47064/ef5fb07e-85ee-28d8-2846-8e32c9b4087e.png">

#### 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)に変換している
 

+ Stanコード
[ex4.stan](https://github.com/MatsuuraKentaro/RStanBook/blob/master/chap05/exercise/ex4.stan)
+ Rコード
[ex4.R](https://github.com/MatsuuraKentaro/RStanBook/blob/master/chap05/exercise/ex4.R)

### 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')
```

図で見ても収束していることがわかる
<img width="472" alt="fit-traceplot2.pdf" src="https://qiita-image-store.s3.amazonaws.com/0/47064/d6daaeb7-c200-b1f2-7b55-fffe676130a1.png">

### 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に推定結果を当てはめる

```math: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に推定結果を当てはめる

```math: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