Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第49回 ロジスティック回帰分析(R実践)

第49回はロジスティック回帰分析をRで実践していきます。データセットのサンプルは、以前の記事でも使わせていただいた「マンガでわかる統計学〜回帰分析〜」をもとにしております。書籍では、ノルンという喫茶店の「スペシャルケーキ」が売れるかどうかをロジスティック回帰分析で予測しています。では、始めていきましょう。 

・glm関数でのロジスティック回帰分析

 glm関数でロジスティック回帰分析を行う場合、目的変数の与え方はいくつかあります。例えば、割合と分母で示す方法や、一般的な0−1の変数にして回帰するなどです。ここでは、一般的な0−1の目的変数を利用します。以下のデータでは、ケーキが売れた場合=1、売れない場合=0で表現してあります。glm関数では、以下のように引数を指定します。

glm(cake ~ day + temp, data=log, family=binomial(link="logit"))

family=binomialは、誤差構造は「二項分布」を利用することを指定、

link="logit"は、リンク関数は「ロジット」を利用することを指定という意味です。

・ロジスティック回帰分析の説明変数の解釈

log <- read.csv("log.csv",header=TRUE)
head(log)
  day temp cake
1   0   28    1
2   0   24    0
3   1   26    0
4   0   24    0
5   0   23    0
6   1   28    1
attach(log)
par(mfrow = c(1, 2))
plot(jitter(log$cake,0.1),jitter(log$day,0.1))
plot(jitter(log$cake,0.1),jitter(log$temp,0.1))

 

f:id:teruaki-sugiura:20150829010033p:plain

 cor(log)
           day           temp         cake
day   1.0000000 0.2013562 0.5095247
temp 0.2013562 1.0000000 0.4828045
cake 0.5095247 0.4828045 1.0000000
result.all <-  glm(cake ~ day + temp, data=log, family=binomial(link="logit"))
summary(result.all)
Call:
glm(formula = cake ~ day + temp, family = binomial(link = "logit"), data = log)
Deviance Residuals:
    Min       1Q       Median       3Q      Max 
-1.7986  -0.6090  -0.2793   0.5180   1.7673 
Coefficients:
                   Estimate    Std. Error   z value    Pr(>|z|) 
(Intercept)   -15.2035    7.6567       -1.986      0.0471 *
day              2.4426       1.2405        1.969      0.0489 *
temp            0.5445       0.2969        1.834      0.0666 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 27.910  on 20  degrees of freedom
Residual deviance: 17.802  on 18  degrees of freedom
AIC: 23.802
Number of Fisher Scoring iterations: 5
exp(2.44) #dayのodds
[1] 11.47304
exp(0.54) #tempのodds
[1] 1.716007
round(exp(result.all$coefficients),2)
(Intercept)         day        temp
       0.00       11.50        1.72
以上の結果より、他の変数が一定という条件で各変数が 1 増加した時、つまり該当する曜日の場合、そうでない曜日に比べodds比が11倍高くなり、tempの場合は1.72倍高くなる。tempが5度上昇すると、odds比はexp{0.54*5}=14.87となります。
par(mfrow = c(2, 2))
plot(result.all)

f:id:teruaki-sugiura:20150821011612p:plain

・誤判別率の計算

fitted関数もしくはpredict関数を用いて、ケーキが売れる確率を求めます。その後に、ifelse関数で0.5以上(ケーキが売れる)は1、そうではない場合は0を割り当て、table関数から計算します。

fit<- fitted(result.all)
1:   0.51
2:   0.10
3:   0.80
4:   0.10
5:   0.06
6:   0.92
7:   0.57
8:   0.25
9:   0.16
10: 0.92
11: 0.02
12: 0.03
13: 0.87
14: 0.80
15: 0.25
16: 0.02
17: 0.20
18: 0.37
19: 0.06
20: 0.31
21: 0.57
log2 <- cbind(log2,fit)
predict <- ifelse(log2$fit>0.5,1,0) #0.5以上=1、それ以外=0
log2 <- cbind(log2, predict)
table(log2$cake, log2$predict)
    0   1
 0  11   2
 1   2   6

 実際に観測されたケーキの販売数と予測値(0.5以上は売れた)を比較すると、4/21(事実とは異なる/n)で計算すると、誤判別率がわかります。つまり誤判別率は0.19となります(正答率は0.81=81%)。この値が小さいほど、回帰モデルの精度が高いことになります。

・モデル比較

①day,tempを入れたモデル

②dayのみのモデル

③tempのみのモデル

④interceptのみ

result.all<- glm(cake ~ day + temp,data=log,family=binomial(link="logit"))
result.day <-  glm(cake ~ day , data=log, family=binomial(link="logit"))
result.temp <-  glm(cake ~ temp , data=log, family=binomial(link="logit"))
result.inter<- glm(cake ~1,data=log,family=binomial(link="logit"))
AIC(result.all,result.day,result.temp,result.inter)
                df      AIC
result.all      3     23.80209
result.day      2     26.27072
result.temp    2     26.53560
result.inter     1     29.91019 

上記の3つのモデルを比較してみた結果、①day+tempを入れたモデルのAICが低く「より良い予測」をするモデルとして選ばれました。

 

以上でおしまい。