My Notes

統計学とかR(R言語)とかPython3の覚え書きとか走り書きとか。 座右の銘にしたい: All work and no play makes Jack a dull boy.

『多変量解析がわかる(ファーストブック)』 第2章 回帰分析 単回帰分析をR(R言語)で

Rコード

# 目的変量 婚姻率 y
# 説明変量 老年人口の割合 x
# 回帰方程式 y = a + bx
# aを切片、bを単回帰係数(回帰係数)




# 47都道府県ごとの婚姻率(人口千人当たりの結婚数)と老年人口の割合(65歳以上の割合)
df <- data.frame(都道府県 = c("北海道", "青森県", "岩手県", "宮城県", "秋田県",
                            "山形県", "福島県", "茨城県", "栃木県", "群馬県",
                            "埼玉県", "千葉県", "東京都", "神奈川県", "新潟県",
                            "富山県", "石川県", "福井県", "山梨県", "長野県",
                            "岐阜県", "静岡県", "愛知県", "三重県", "滋賀県",
                            "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県",
                            "鳥取県", "島根県", "岡山県", "広島県", "山口県",
                            "徳島県", "香川県", "愛媛県", "高知県", "福岡県",
                            "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県",
                            "鹿児島県", "沖縄県"),
                老年人口 = c(23.6, 24.4, 26.3, 21.5, 28.4,
                            26.6, 24.2, 21.3, 21.1, 22.5,
                            19.1, 20.1, 20.2, 19.2, 25.5,
                            25.2, 22.9, 24.3, 23.7, 25.5,
                            22.9, 22.6, 19.2, 23.1, 19.7,
                            22.4, 21.2, 22.1, 22.6, 26.1,
                            25.5, 28.6, 24.3, 23.0, 26.9,
                            26.1, 24.9, 25.6, 27.8, 21.4,
                            23.9, 25.2, 25.1, 25.9, 25.2,
                            26.0, 17.2),
                婚姻率 = c(5.17, 4.55, 4.66, 5.46, 4.00,
                        4.56, 4.92, 5.25, 5.52, 5.14,
                        5.68, 5.86, 6.99, 6.36, 4.65,
                        4.69, 5.12, 5.05, 5.08, 5.11,
                        5.08, 5.56, 6.38, 5.29, 5.65,
                        5.30, 5.90, 5.45, 4.90, 4.87,
                        4.80, 4.38, 5.19, 5.62, 4.93,
                        4.69, 5.22, 5.03, 4.54, 5.83,
                        4.90, 4.80, 5.17, 5.25, 5.47,
                        5.05, 6.28))

df


老年人口の割合x <- df$老年人口
老年人口の割合x


婚姻率y <- df$婚姻率
婚姻率y




# 以下、分析していく
# 文字化けするなら
par(family = "Osaka")


# 散布図
plot(老年人口の割合x, 婚姻率y)


# 相関係数
cor(老年人口の割合x, 婚姻率y)
# 相関係数の検定(ピアソン)
cor.test(老年人口の割合x, 婚姻率y)


# 回帰式
lm1 <- lm(婚姻率y ~ 老年人口の割合x, data = df)
lm1


# 回帰直線
abline(lm1)


summary(lm1)


# 使用したデータ、予測値(推測値)、残差
# 予測値(推測値)
予測値 <- predict(lm1)
# 残差
残差 <- residuals(lm1)
data.frame(df, 予測値, 残差)


# 回帰診断図
par(mfrow = c(2, 2), oma = c(2, 2, 2, 1), mar = c(4, 4, 3, 2))
plot(lm1)




#
# 負の相関なので、以下、手計算上の値の取り扱い方にちょっと気をつけることがある
#




# 回帰方程式に具体的な値を付値(代入)
# y = a - bx
# y = 9.53 - 0.18x
# 切片a
a <- round(coef(lm1)[[1]], 2)
a
# 単回帰係数b
b <- round(coef(lm1)[[2]], 2)
b
# y = a - b*x となる


# 例えば、老年人口の割合が25%の県の婚姻率はどれくらいか
y1 <- a - (-b)*25 # 負の値は正の値にしておく。この場合、bは負だから、(-b)とした
y1
# abs()を使用した場合
y2 <- a - abs(b)*25
y2


# 26%の場合
y3 <- a - abs(b)*26
y3


# xが1%増加するとyは0.18下がる
y2 - y3
y3 - y2


# 30%の場合
y4 <- a - abs(b)*30
y4




#
# 決定係数等の手計算は割愛
#




#
# 結論 この単回帰分析は、税収の見込みや人口予測に利用できることもあるかもしれない
#

R Console

> # 目的変量 婚姻率 y
> # 説明変量 老年人口の割合 x
> # 回帰方程式 y = a + bx
> # aを切片、bを単回帰係数(回帰係数)
> 
> 
> 
> 
> # 47都道府県ごとの婚姻率(人口千人当たりの結婚数)と老年人口の割合(65歳以上の割合)
> df <- data.frame(都道府県 = c("北海道", "青森県", "岩手県", "宮城県", "秋田県",
+                             "山形県", "福島県", "茨城県", "栃木県", "群馬県",
+                             "埼玉県", "千葉県", "東京都", "神奈川県", "新潟県",
+                             "富山県", "石川県", "福井県", "山梨県", "長野県",
+                             "岐阜県", "静岡県", "愛知県", "三重県", "滋賀県",
+                             "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県",
+                             "鳥取県", "島根県", "岡山県", "広島県", "山口県",
+                             "徳島県", "香川県", "愛媛県", "高知県", "福岡県",
+                             "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県",
+                             "鹿児島県", "沖縄県"),
+                 老年人口 = c(23.6, 24.4, 26.3, 21.5, 28.4,
+                             26.6, 24.2, 21.3, 21.1, 22.5,
+                             19.1, 20.1, 20.2, 19.2, 25.5,
+                             25.2, 22.9, 24.3, 23.7, 25.5,
+                             22.9, 22.6, 19.2, 23.1, 19.7,
+                             22.4, 21.2, 22.1, 22.6, 26.1,
+                             25.5, 28.6, 24.3, 23.0, 26.9,
+                             26.1, 24.9, 25.6, 27.8, 21.4,
+                             23.9, 25.2, 25.1, 25.9, 25.2,
+                             26.0, 17.2),
+                 婚姻率 = c(5.17, 4.55, 4.66, 5.46, 4.00,
+                         4.56, 4.92, 5.25, 5.52, 5.14,
+                         5.68, 5.86, 6.99, 6.36, 4.65,
+                         4.69, 5.12, 5.05, 5.08, 5.11,
+                         5.08, 5.56, 6.38, 5.29, 5.65,
+                         5.30, 5.90, 5.45, 4.90, 4.87,
+                         4.80, 4.38, 5.19, 5.62, 4.93,
+                         4.69, 5.22, 5.03, 4.54, 5.83,
+                         4.90, 4.80, 5.17, 5.25, 5.47,
+                         5.05, 6.28))
> 
> df
   都道府県 老年人口 婚姻率
1    北海道     23.6   5.17
2    青森県     24.4   4.55
3    岩手県     26.3   4.66
4    宮城県     21.5   5.46
5    秋田県     28.4   4.00
6    山形県     26.6   4.56
7    福島県     24.2   4.92
8    茨城県     21.3   5.25
9    栃木県     21.1   5.52
10   群馬県     22.5   5.14
11   埼玉県     19.1   5.68
12   千葉県     20.1   5.86
13   東京都     20.2   6.99
14 神奈川県     19.2   6.36
15   新潟県     25.5   4.65
16   富山県     25.2   4.69
17   石川県     22.9   5.12
18   福井県     24.3   5.05
19   山梨県     23.7   5.08
20   長野県     25.5   5.11
21   岐阜県     22.9   5.08
22   静岡県     22.6   5.56
23   愛知県     19.2   6.38
24   三重県     23.1   5.29
25   滋賀県     19.7   5.65
26   京都府     22.4   5.30
27   大阪府     21.2   5.90
28   兵庫県     22.1   5.45
29   奈良県     22.6   4.90
30 和歌山県     26.1   4.87
31   鳥取県     25.5   4.80
32   島根県     28.6   4.38
33   岡山県     24.3   5.19
34   広島県     23.0   5.62
35   山口県     26.9   4.93
36   徳島県     26.1   4.69
37   香川県     24.9   5.22
38   愛媛県     25.6   5.03
39   高知県     27.8   4.54
40   福岡県     21.4   5.83
41   佐賀県     23.9   4.90
42   長崎県     25.2   4.80
43   熊本県     25.1   5.17
44   大分県     25.9   5.25
45   宮崎県     25.2   5.47
46 鹿児島県     26.0   5.05
47   沖縄県     17.2   6.28
> 
> 
> 老年人口の割合x <- df$老年人口
> 老年人口の割合x
 [1] 23.6 24.4 26.3 21.5 28.4 26.6 24.2 21.3 21.1 22.5 19.1 20.1 20.2 19.2 25.5 25.2 22.9 24.3 23.7
[20] 25.5 22.9 22.6 19.2 23.1 19.7 22.4 21.2 22.1 22.6 26.1 25.5 28.6 24.3 23.0 26.9 26.1 24.9 25.6
[39] 27.8 21.4 23.9 25.2 25.1 25.9 25.2 26.0 17.2
> 
> 
> 婚姻率y <- df$婚姻率
> 婚姻率y
 [1] 5.17 4.55 4.66 5.46 4.00 4.56 4.92 5.25 5.52 5.14 5.68 5.86 6.99 6.36 4.65 4.69 5.12 5.05 5.08
[20] 5.11 5.08 5.56 6.38 5.29 5.65 5.30 5.90 5.45 4.90 4.87 4.80 4.38 5.19 5.62 4.93 4.69 5.22 5.03
[39] 4.54 5.83 4.90 4.80 5.17 5.25 5.47 5.05 6.28
> 
> 
> 
> 
> # 以下、分析していく
> # 文字化けするなら
> par(family = "Osaka")
> 
> 
> # 散布図
> plot(老年人口の割合x, 婚姻率y)
> 
> 
> # 相関係数
> cor(老年人口の割合x, 婚姻率y)
[1] -0.845768
> # 相関係数の検定(ピアソン)
> cor.test(老年人口の割合x, 婚姻率y)

    Pearson's product-moment correlation

data:  老年人口の割合x and 婚姻率y
t = -10.634, df = 45, p-value = 7.328e-14
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9115425 -0.7377940
sample estimates:
      cor 
-0.845768 

> 
> 
> # 回帰式
> lm1 <- lm(婚姻率y ~ 老年人口の割合x, data = df)
> lm1

Call:
lm(formula = 婚姻率y ~ 老年人口の割合x, data = df)

Coefficients:
    (Intercept)  老年人口の割合x  
         9.5319          -0.1825  

> 
> 
> # 回帰直線
> abline(lm1)
> 
> 
> summary(lm1)

Call:
lm(formula = 婚姻率y ~ 老年人口の割合x, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.52767 -0.21052 -0.05371  0.21240  1.14563 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      9.53186    0.40793   23.37  < 2e-16 ***
老年人口の割合x -0.18255    0.01717  -10.63 7.33e-14 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3066 on 45 degrees of freedom
Multiple R-squared:  0.7153, Adjusted R-squared:  0.709 
F-statistic: 113.1 on 1 and 45 DF,  p-value: 7.328e-14

> 
> 
> # 使用したデータ、予測値(推測値)、残差
> # 予測値(推測値)
> 予測値 <- predict(lm1)
> # 残差
> 残差 <- residuals(lm1)
> data.frame(df, 予測値, 残差)
   都道府県 老年人口 婚姻率   予測値         残差
1    北海道     23.6   5.17 5.223708 -0.053708384
2    青森県     24.4   4.55 5.077669 -0.527669222
3    岩手県     26.3   4.66 4.730826 -0.070826211
4    宮城県     21.5   5.46 5.607061 -0.147061185
5    秋田県     28.4   4.00 4.347473 -0.347473410
6    山形県     26.6   4.56 4.676062 -0.116061526
7    福島県     24.2   4.92 5.114179 -0.194179012
8    茨城県     21.3   5.25 5.643571 -0.393570976
9    栃木県     21.1   5.52 5.680081 -0.160080766
10   群馬県     22.5   5.14 5.424512 -0.284512232
11   埼玉県     19.1   5.68 6.045179 -0.365178672
12   千葉県     20.1   5.86 5.862630 -0.002629719
13   東京都     20.2   6.99 5.844375  1.145625176
14 神奈川県     19.2   6.36 6.026924  0.333076223
15   新潟県     25.5   4.65 4.876865 -0.226865374
16   富山県     25.2   4.69 4.931630 -0.241630060
17   石川県     22.9   5.12 5.351493 -0.231492651
18   福井県     24.3   5.05 5.095924 -0.045924117
19   山梨県     23.7   5.08 5.205453 -0.125453489
20   長野県     25.5   5.11 4.876865  0.233134626
21   岐阜県     22.9   5.08 5.351493 -0.271492651
22   静岡県     22.6   5.56 5.406257  0.153742663
23   愛知県     19.2   6.38 6.026924  0.353076223
24   三重県     23.1   5.29 5.314983 -0.024982861
25   滋賀県     19.7   5.65 5.935649 -0.285649300
26   京都府     22.4   5.30 5.442767 -0.142767128
27   大阪府     21.2   5.90 5.661826  0.238174129
28   兵庫県     22.1   5.45 5.497532 -0.047531814
29   奈良県     22.6   4.90 5.406257 -0.506257337
30 和歌山県     26.1   4.87 4.767336  0.102663998
31   鳥取県     25.5   4.80 4.876865 -0.076865374
32   島根県     28.6   4.38 4.310964  0.069036380
33   岡山県     24.3   5.19 5.095924  0.094075883
34   広島県     23.0   5.62 5.333238  0.286762244
35   山口県     26.9   4.93 4.621297  0.308703160
36   徳島県     26.1   4.69 4.767336 -0.077336002
37   香川県     24.9   5.22 4.986395  0.233605255
38   愛媛県     25.6   5.03 4.858610  0.171389522
39   高知県     27.8   4.54 4.457003  0.082997218
40   福岡県     21.4   5.83 5.625316  0.204683919
41   佐賀県     23.9   4.90 5.168944 -0.268943698
42   長崎県     25.2   4.80 4.931630 -0.131630060
43   熊本県     25.1   5.17 4.949885  0.220115045
44   大分県     25.9   5.25 4.803846  0.446154207
45   宮崎県     25.2   5.47 4.931630  0.538369940
46 鹿児島県     26.0   5.05 4.785591  0.264409103
47   沖縄県     17.2   6.28 6.392022 -0.112021683
> 
> 
> # 回帰診断図
> par(mfrow = c(2, 2), oma = c(2, 2, 2, 1), mar = c(4, 4, 3, 2))
> plot(lm1)
>
>
>
> 
> #
> # 負の相関なので、以下、手計算上の値の取り扱い方にちょっと気をつけることがある
> #
> 
> 
> 
> 
> # 回帰方程式に具体的な値を付値(代入)
> # y = a - bx
> # y = 9.53 - 0.18x
> # 切片a
> a <- round(coef(lm1)[[1]], 2)
> a
[1] 9.53
> # 単回帰係数b
> b <- round(coef(lm1)[[2]], 2)
> b
[1] -0.18
> # y = a - b*x となる
> 
> 
> # 例えば、老年人口の割合が25%の県の婚姻率はどれくらいか
> y1 <- a - (-b)*25 # 負の値は正の値にしておく。この場合、bは負だから、(-b)とした
> y1
[1] 5.03
> # abs()を使用した場合
> y2 <- a - abs(b)*25
> y2
[1] 5.03
> 
> 
> # 26%の場合
> y3 <- a - abs(b)*26
> y3
[1] 4.85
> 
> 
> # xが1%増加するとyは0.18下がる
> y2 - y3
[1] 0.18
> y3 - y2
[1] -0.18
> 
> 
> # 30%の場合
> y4 <- a - abs(b)*30
> y4
[1] 4.13
> 
> 
> 
> 
> #
> # 決定係数等の手計算は割愛
> #
> 
> 
> 
> 
> #
> # 結論 この単回帰分析は、税収の見込みや人口予測に利用できることもあるかもしれない
> #

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

f:id:my_notes:20170615153122p:plain

回帰直線のスクリーンショット

f:id:my_notes:20170615153149p:plain

回帰診断図のスクリーンショット

f:id:my_notes:20170615181038p:plain

参考文献

多変量解析がわかる (ファーストブック)

多変量解析がわかる (ファーストブック)

Rコードの参考文献

Rによるデータサイエンス データ解析の基礎から最新手法まで

Rによるデータサイエンス データ解析の基礎から最新手法まで

新版が出ている

Rによるデータサイエンス データ解析の基礎から最新手法まで

Rによるデータサイエンス データ解析の基礎から最新手法まで

あまりうまくいかなかった単回帰分析の記事

my-notes.hatenablog.com