読者です 読者をやめる 読者になる 読者になる

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が低く「より良い予測」をするモデルとして選ばれました。

 

以上でおしまい。

広告を非表示にする