『マンガでわかる統計学 回帰分析編』第4章 ロジスティック回帰分析 をR(R言語)で

Rコード

#
# 『マンガでわかる統計学 回帰分析編』 第4章 ロジスティック回帰分析
#




# 水曜or土曜or日曜 1が水曜or土曜or日曜、 0はそれ以外
# ノルンスペシャルの販売状況 1は売れたということ、0は売れなかったということ
df <- data.frame(水曜or土曜or日曜 = c(0, 0, 1, 0, 0,
                                    1, 1, 0, 0, 1,
                                    0, 0, 1, 1, 0,
                                    0, 1, 0, 0, 1, 1),
                最高気温 = c(28, 24, 26, 24, 23,
                            28, 24, 26, 25, 28,
                            21, 22, 27, 26, 26,
                            21, 21, 27, 23, 22, 24),
                ノルンスペシャルの販売状況 = c(1, 0, 0, 0, 0,
                                        1, 0, 1, 0, 1,
                                        0, 0, 1, 1, 0,
                                        0, 1, 0, 0, 0, 1))

df


水曜or土曜or日曜 <- df$水曜or土曜or日曜
水曜or土曜or日曜

最高気温 <- df$最高気温
最高気温

ノルンスペシャルの販売状況 <- df$ノルンスペシャルの販売状況
ノルンスペシャルの販売状況




summary(df)




# 散布図と相関係数
par(family = "Osaka") # 文字化けするなら。

plot(ノルンスペシャルの販売状況, 水曜or土曜or日曜)
# 重なってしまうので、sunflowerplot()
sunflowerplot(ノルンスペシャルの販売状況, 水曜or土曜or日曜)
cor.test(ノルンスペシャルの販売状況, 水曜or土曜or日曜)


plot(ノルンスペシャルの販売状況, 最高気温)
# 重なってしまうので、sunflowerplot()
sunflowerplot(ノルンスペシャルの販売状況, 最高気温)
cor.test(ノルンスペシャルの販売状況, 最高気温)




# ロジスティック回帰分析
glm1 <- glm(ノルンスペシャルの販売状況 ~ 水曜or土曜or日曜+最高気温, family = binomial(link = "logit"), data = df)
glm1
summary(glm1)




予測値 <- glm1$fitted.values
data.frame(水曜or土曜or日曜, 最高気温, ノルンスペシャルの販売状況, 予測値)




# 以下、予測する
#
# 数式 y = 1 / (1 + exp(-(a1x1 + a2x2 + b)))
#
# 数式に利用するため、変数に付値(代入)
b <- round(coef(glm1)[[1]], 2)
b

a1 <- round(coef(glm1)[[2]], 2)
a1

a2 <- round(coef(glm1)[[3]], 2)
a2


# 予測するための条件
# 今日が日曜日、最高気温23度
# 0.5以上なら、売れる見込みがあるといえるかもしれない
#
y1 <- 1 / (1 + exp(-(a1*1 + a2*23 + b)))
y1
round(y1, 2)


# exp()を少数第2位まで丸めて、同じ計算
e <- round(exp(1), 2)
y2 <- 1 / (1 + e^(-(a1*1 + a2*23 + b)))
y2
round(y2, 2)




#
# 結論 今日は売れる見込みがあるとはいえないかもしれない
#




#
# なお、少数について
# 書籍(マンガ)では、0.44となっているが、
# R(R言語)では、0.4158095で約0.42、手持ちの関数電卓でも0.4158094771で約0.42(round(0.4158094771, 2))となった。
#
# 書籍(マンガ)冒頭に、少数が一致しないことの記載がある。
#

R Console

> #
> # 『マンガでわかる統計学 回帰分析編』 第4章 ロジスティック回帰分析
> #
> 
> 
> 
> 
> # 水曜or土曜or日曜 1が水曜or土曜or日曜、 0はそれ以外
> # ノルンスペシャルの販売状況 1は売れたということ、0は売れなかったということ
> df <- data.frame(水曜or土曜or日曜 = c(0, 0, 1, 0, 0,
+                                     1, 1, 0, 0, 1,
+                                     0, 0, 1, 1, 0,
+                                     0, 1, 0, 0, 1, 1),
+                 最高気温 = c(28, 24, 26, 24, 23,
+                             28, 24, 26, 25, 28,
+                             21, 22, 27, 26, 26,
+                             21, 21, 27, 23, 22, 24),
+                 ノルンスペシャルの販売状況 = c(1, 0, 0, 0, 0,
+                                         1, 0, 1, 0, 1,
+                                         0, 0, 1, 1, 0,
+                                         0, 1, 0, 0, 0, 1))
> 
> df
   水曜or土曜or日曜 最高気温 ノルンスペシャルの販売状況
1                 0       28                          1
2                 0       24                          0
3                 1       26                          0
4                 0       24                          0
5                 0       23                          0
6                 1       28                          1
7                 1       24                          0
8                 0       26                          1
9                 0       25                          0
10                1       28                          1
11                0       21                          0
12                0       22                          0
13                1       27                          1
14                1       26                          1
15                0       26                          0
16                0       21                          0
17                1       21                          1
18                0       27                          0
19                0       23                          0
20                1       22                          0
21                1       24                          1
> 
> 
> 水曜or土曜or日曜 <- df$水曜or土曜or日曜
> 水曜or土曜or日曜
 [1] 0 0 1 0 0 1 1 0 0 1 0 0 1 1 0 0 1 0 0 1 1
> 
> 最高気温 <- df$最高気温
> 最高気温
 [1] 28 24 26 24 23 28 24 26 25 28 21 22 27 26 26 21 21 27 23 22 24
> 
> ノルンスペシャルの販売状況 <- df$ノルンスペシャルの販売状況
> ノルンスペシャルの販売状況
 [1] 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 0 1 0 0 0 1
> 
> 
> 
> 
> summary(df)
 水曜or土曜or日曜    最高気温     ノルンスペシャルの販売状況
 Min.   :0.0000   Min.   :21.00   Min.   :0.000             
 1st Qu.:0.0000   1st Qu.:23.00   1st Qu.:0.000             
 Median :0.0000   Median :24.00   Median :0.000             
 Mean   :0.4286   Mean   :24.57   Mean   :0.381             
 3rd Qu.:1.0000   3rd Qu.:26.00   3rd Qu.:1.000             
 Max.   :1.0000   Max.   :28.00   Max.   :1.000             
> 
> 
> 
> 
> # 散布図と相関係数
> par(family = "Osaka") # 文字化けするなら。
> 
> plot(ノルンスペシャルの販売状況, 水曜or土曜or日曜)
> # 重なってしまうので、sunflowerplot()
> sunflowerplot(ノルンスペシャルの販売状況, 水曜or土曜or日曜)
> cor.test(ノルンスペシャルの販売状況, 水曜or土曜or日曜)

    Pearson's product-moment correlation

data:  ノルンスペシャルの販売状況 and 水曜or土曜or日曜
t = 2.5812, df = 19, p-value = 0.01831
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.09978641 0.77151321
sample estimates:
      cor 
0.5095247 

> 
> 
> plot(ノルンスペシャルの販売状況, 最高気温)
> # 重なってしまうので、sunflowerplot()
> sunflowerplot(ノルンスペシャルの販売状況, 最高気温)
> cor.test(ノルンスペシャルの販売状況, 最高気温)

    Pearson's product-moment correlation

data:  ノルンスペシャルの販売状況 and 最高気温
t = 2.4031, df = 19, p-value = 0.02663
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.06457691 0.75676592
sample estimates:
      cor 
0.4828045 

> 
> 
> 
> 
> # ロジスティック回帰分析
> glm1 <- glm(ノルンスペシャルの販売状況 ~ 水曜or土曜or日曜+最高気温, family = binomial(link = "logit"), data = df)
> glm1

Call:  glm(formula = ノルンスペシャルの販売状況 ~ 水曜or土曜or日曜 + 
    最高気温, family = binomial(link = "logit"), data = df)

Coefficients:
     (Intercept)  水曜or土曜or日曜          最高気温  
        -15.2035            2.4426            0.5445  

Degrees of Freedom: 20 Total (i.e. Null);  18 Residual
Null Deviance:      27.91 
Residual Deviance: 17.8    AIC: 23.8
> summary(glm1)

Call:
glm(formula = ノルンスペシャルの販売状況 ~ 水曜or土曜or日曜 + 
    最高気温, family = binomial(link = "logit"), data = df)

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 *
水曜or土曜or日曜   2.4426     1.2405   1.969   0.0489 *
最高気温           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

> 
> 
> 
> 
> 予測値 <- glm1$fitted.values
> data.frame(水曜or土曜or日曜, 最高気温, ノルンスペシャルの販売状況, 予測値)
   水曜or土曜or日曜 最高気温 ノルンスペシャルの販売状況     予測値
1                 0       28                          1 0.51066015
2                 0       24                          0 0.10570583
3                 1       26                          0 0.80159346
4                 0       24                          0 0.10570583
5                 0       23                          0 0.06417107
6                 1       28                          1 0.92310444
7                 1       24                          0 0.57621937
8                 0       26                          1 0.25992396
9                 0       25                          0 0.16926163
10                1       28                          1 0.92310444
11                0       21                          0 0.02255708
12                0       22                          0 0.03825829
13                1       27                          1 0.87443878
14                1       26                          1 0.80159346
15                0       26                          0 0.25992396
16                0       21                          0 0.02255708
17                1       21                          1 0.20978135
18                0       27                          0 0.37710404
19                0       23                          0 0.06417107
20                1       22                          0 0.31394534
21                1       24                          1 0.57621937
> 
> 
> 
> 
> # 以下、予測する
> #
> # 数式 y = 1 / (1 + exp(-(a1x1 + a2x2 + b)))
> #
> # 数式に利用するため、変数に付値(代入)
> b <- round(coef(glm1)[[1]], 2)
> b
[1] -15.2
> 
> a1 <- round(coef(glm1)[[2]], 2)
> a1
[1] 2.44
> 
> a2 <- round(coef(glm1)[[3]], 2)
> a2
[1] 0.54
> 
> 
> # 予測するための条件
> # 今日が日曜日、最高気温23度
> # 0.5以上なら、売れる見込みがあるといえるかもしれない
> #
> y1 <- 1 / (1 + exp(-(a1*1 + a2*23 + b)))
> y1
[1] 0.4158095
> round(y1, 2)
[1] 0.42
> 
> 
> # exp()を少数第2位まで丸めて、同じ計算
> e <- round(exp(1), 2)
> y2 <- 1 / (1 + e^(-(a1*1 + a2*23 + b)))
> y2
[1] 0.4157573
> round(y2, 2)
[1] 0.42
> 
> 
> 
> 
> #
> # 結論 今日は売れる見込みがあるとはいえないかもしれない
> #
> 
> 
> 
> 
> #
> # なお、少数について
> # 書籍(マンガ)では、0.44となっているが、
> # R(R言語)では、0.4158095で約0.42、手持ちの関数電卓でも0.4158094771で約0.42(round(0.4158094771, 2))となった。
> #
> # 書籍(マンガ)冒頭に、少数が一致しないことの記載がある。
> #

散布図のスクリーンショット

f:id:my_notes:20170623162215p:plain

f:id:my_notes:20170623162233p:plain

f:id:my_notes:20170623162249p:plain

f:id:my_notes:20170623162303p:plain

参考文献

マンガでわかる統計学 回帰分析編

マンガでわかる統計学 回帰分析編

Rによる統計解析

Rによる統計解析