Lesson3.Regression
Kosugitti
2015年11月6日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
野球のデータ，baseball2015.csvをb.dataというオブジェクトに入れているとします。 できてないひとは，次の通りやってみよう。
b.data <- read.csv("baseball2015.csv",na.strings="*")回帰分析
回帰分析とは，ある目的となる変数を説明する関数を当てはめることです。一般に，\[y=ax+b\] という線形関数を当てはめます。
身長と体重を例にして説明しましょう。 この二つの変数の散布図を書いてみます。
library(ggplot2)
g <- ggplot(b.data,aes(x=height,y=weight)) + geom_point()
plot(g)
右肩あがりの関係がありそうです。線を入れてみます。
g <- g + geom_smooth(method="lm")
plot(g)
この線の引き方が，線形回帰をするということです。回帰モデルをデータに当てはめる，とも言います。 lmは，liner model(線形モデル)という意味です。
この線形モデルを数式で表現するのが次の方法です。
reg <- lm(weight~height,data=b.data)
reg##
## Call:
## lm(formula = weight ~ height, data = b.data)
##
## Coefficients:
## (Intercept) height
## -141.764 1.265すなわち，体重＝1.265 * 身長 -141.764 という関係にあるということがわかります。言い換えると，身長が1cm高いと体重が1.2kg大きいということです(このデータの中では）。また。Intercept(切片)は身長が0だった時の(!)値です。
もちろん，すべてのデータがこの直線上に乗っているわけではありません。むしろほとんど外れていますが，この線が「当てはまっている」といえるのはどうしてでしょうか？ここでは，最小二乗基準 によってこの数字を決めています。最小二乗法による回帰分析，ということもあります。 推定値の決め方は他にもある，ということを頭の片隅に置いておいてください。
さて，最小二乗基準で推定値が得られはしましたが，誤差も多いので，これをどう評価すればよいでしょうか。 ちなみに，推定値や誤差は結果のオブジェクトに入っています。
reg$residuals #誤差。残差ともいう。## 1 2 3 4 5
## 2.46407731 -8.85879020 -0.44623122 3.34748929 0.81834229
## 6 7 8 9 10
## 9.87663630 -2.91708421 8.67035681 -7.38793721 -2.12336370
## 11 12 13 14 15
## -9.44623122 12.37888675 -11.85879020 -7.18165771 10.23090127
## 16 17 18 19 20
## 1.14120980 -9.97537822 3.58066534 -7.85879020 -8.74220218
## 21 22 23 24 25
## 11.67035681 11.23090127 -2.47312777 -5.29824574 4.23090127
## 26 27 28 29 30
## 4.23090127 -3.62111325 -10.94398076 -11.85879020 -1.65251071
## 31 32 33 34 35
## -7.12336370 3.08291579 -1.18165771 -1.32964319 -2.59421670
## 36 37 38 39 40
## -4.85879020 10.14120980 -13.38793721 -17.12336370 2.02462178
## 41 42 43 44 45
## -2.18165771 -7.80049619 2.61206279 -0.29824574 1.08291579
## 46 47 48 49 50
## -7.12336370 1.14120980 -6.32964319 -5.85879020 6.70175426
## 51 52 53 54 55
## 10.87663630 2.28919528 7.08291579 -20.62111325 -6.09196624
## 56 57 58 59 60
## -0.21305517 -5.41483376 -2.06506969 -9.29824574 0.93493031
## 61 62 63 64 65
## -5.85879020 5.19950381 -11.97537822 -0.65251071 5.52237133
## 66 67 68 69 70
## 0.37888675 -5.06506969 -0.44623122 -6.38793721 8.40578330
## 71 72 73 74 75
## -2.06506969 -3.32964319 -7.12336370 -5.97537822 2.02462178
## 76 77 78 79 80
## -20.82739274 -0.06506969 -3.85879020 -9.97537822 9.81834229
## 81 82 83 84 85
## 2.58066534 1.19950381 7.08291579 -2.06506969 8.34748929
## 86 87 88 89 90
## -0.44623122 8.46407731 9.14120980 -7.59421670 -0.85879020
## 91 92 93 94 95
## -2.12336370 -11.65251071 -1.59421670 0.34748929 -8.06506969
## 96 97 98 99 100
## -15.23995172 -1.80049619 19.34748929 -4.59421670 -5.00677568
## 101 102 103 104 105
## 11.14120980 6.81834229 26.43718076 -8.85879020 22.46407731
## 106 107 108 109 110
## -0.65251071 -6.71080472 4.14120980 0.81834229 9.93493031
## 111 112 113 114 115
## 7.81834229 14.81834229 15.46407731 6.02462178 2.93493031
## 116 117 118 119 120
## 3.34748929 -10.44623122 3.14120980 8.02462178 -0.85879020
## 121 122 123 124 125
## 2.14120980 4.28919528 10.34748929 3.14120980 1.61206279
## 126 127 128 129 130
## 9.37888675 -6.27134918 -2.65251071 7.67035681 -0.38793721
## 131 132 133 134 135
## -0.85879020 1.14120980 2.93493031 -2.80049619 -2.18165771
## 136 137 138 139 140
## -4.53592269 13.34748929 5.61206279 15.02462178 1.02462178
## 141 142 143 144 145
## -4.94848167 2.87663630 -2.59421670 -4.53592269 -11.85879020
## 146 147 148 149 150
## 15.90803376 -1.18165771 11.93493031 8.46407731 7.93493031
## 151 152 153 154 155
## 13.34748929 0.67035681 0.81834229 -1.00677568 -12.29824574
## 156 157 158 159 160
## -0.80049619 10.14120980 -0.27134918 0.19950381 7.14120980
## 161 162 163 164 165
## 7.93493031 -3.97537822 -5.47762867 -6.18165771 5.34748929
## 166 167 168 169 170
## -9.85879020 2.55376878 1.61206279 -12.59421670 4.46407731
## 171 172 173 174 175
## 2.46407731 -5.06506969 -9.85879020 -6.80049619 2.67035681
## 176 177 178 179 180
## -2.85879020 -0.85879020 18.23090127 -1.47762867 -0.53592269
## 181 182 183 184 185
## -14.06506969 -4.65251071 -1.32964319 -6.32964319 5.81834229
## 186 187 188 189 190
## -1.03367223 2.55376878 10.49547477 -2.38793721 4.72865082
## 191 192 193 194 195
## -11.65251071 -4.00677568 -4.00677568 1.19950381 -5.06506969
## 196 197 198 199 200
## -7.85879020 16.67035681 -2.12336370 -9.59421670 -2.32964319reg$fitted.values #推定値## 1 2 3 4 5 6 7
## 79.53592 85.85879 93.44623 89.65251 92.18166 87.12336 90.91708
## 8 9 10 11 12 13 14
## 83.32964 88.38794 87.12336 93.44623 108.62111 85.85879 92.18166
## 15 16 17 18 19 20 21
## 99.76910 85.85879 95.97538 69.41933 85.85879 75.74220 83.32964
## 22 23 24 25 26 27 28
## 99.76910 117.47313 102.29825 99.76910 99.76910 108.62111 114.94398
## 29 30 31 32 33 34 35
## 85.85879 89.65251 87.12336 90.91708 92.18166 83.32964 84.59422
## 36 37 38 39 40 41 42
## 85.85879 85.85879 88.38794 87.12336 95.97538 92.18166 80.80050
## 43 44 45 46 47 48 49
## 88.38794 102.29825 90.91708 87.12336 85.85879 83.32964 85.85879
## 50 51 52 53 54 55 56
## 102.29825 87.12336 94.71080 90.91708 108.62111 106.09197 73.21306
## 57 58 59 60 61 62 63
## 112.41483 82.06507 102.29825 82.06507 85.85879 80.80050 95.97538
## 64 65 66 67 68 69 70
## 89.65251 74.47763 108.62111 82.06507 93.44623 88.38794 84.59422
## 71 72 73 74 75 76 77
## 82.06507 83.32964 87.12336 95.97538 95.97538 104.82739 82.06507
## 78 79 80 81 82 83 84
## 85.85879 95.97538 92.18166 69.41933 80.80050 90.91708 82.06507
## 85 86 87 88 89 90 91
## 89.65251 93.44623 79.53592 85.85879 84.59422 85.85879 87.12336
## 92 93 94 95 96 97 98
## 89.65251 84.59422 89.65251 82.06507 97.23995 80.80050 89.65251
## 99 100 101 102 103 104 105
## 84.59422 77.00678 85.85879 92.18166 103.56282 85.85879 79.53592
## 106 107 108 109 110 111 112
## 89.65251 94.71080 85.85879 92.18166 82.06507 92.18166 92.18166
## 113 114 115 116 117 118 119
## 79.53592 95.97538 82.06507 89.65251 93.44623 85.85879 95.97538
## 120 121 122 123 124 125 126
## 85.85879 85.85879 94.71080 89.65251 85.85879 88.38794 108.62111
## 127 128 129 130 131 132 133
## 78.27135 89.65251 83.32964 88.38794 85.85879 85.85879 82.06507
## 134 135 136 137 138 139 140
## 80.80050 92.18166 79.53592 89.65251 88.38794 95.97538 95.97538
## 141 142 143 144 145 146 147
## 71.94848 87.12336 84.59422 79.53592 85.85879 106.09197 92.18166
## 148 149 150 151 152 153 154
## 82.06507 79.53592 82.06507 89.65251 83.32964 92.18166 77.00678
## 155 156 157 158 159 160 161
## 102.29825 80.80050 85.85879 78.27135 80.80050 85.85879 82.06507
## 162 163 164 165 166 167 168
## 95.97538 74.47763 92.18166 89.65251 85.85879 93.44623 88.38794
## 169 170 171 172 173 174 175
## 84.59422 79.53592 79.53592 82.06507 85.85879 80.80050 83.32964
## 176 177 178 179 180 181 182
## 85.85879 85.85879 99.76910 74.47763 79.53592 82.06507 89.65251
## 183 184 185 186 187 188 189
## 83.32964 83.32964 92.18166 101.03367 93.44623 98.50453 88.38794
## 190 191 192 193 194 195 196
## 78.27135 89.65251 77.00678 77.00678 80.80050 82.06507 85.85879
## 197 198 199 200
## 83.32964 87.12336 84.59422 83.32964この推定値と実際のデータとの相関係数が当てはまりの良さ(適合度)の一つの指標です。相関係数が高ければ高いほど，良い予測ができていると言えます。 相関係数は-1から1までの値を取りますが，適合度指標としては相関係数の二乗を使います(R-squared)。 この指標も実は結果のオブジェクトに既に入っています。
summary(reg)##
## Call:
## lm(formula = weight ~ height, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.8274 -5.1234 -0.5942 4.3329 26.4372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -141.76444 15.39438 -9.209 <2e-16 ***
## height 1.26457 0.08454 14.958 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.856 on 198 degrees of freedom
## Multiple R-squared: 0.5305, Adjusted R-squared: 0.5281
## F-statistic: 223.7 on 1 and 198 DF, p-value: < 2.2e-160.5ほどですね。これは実際のデータとしては，とても良い当てはまりだと言えます。
うんちく
- summary関数の出力において，推定値がどこにあるか確認してください。
- アスタリスクマークが出ていますが，これはこの回帰係数が「統計的に有意である」ことを示しています。
- R-squaredは\(R^2\)とも書き，説明される分散の大きさを表しています。
- F統計量が出ています。これは回帰分析が分散分析と数学的には同じであることの表れです！
重回帰分析
説明する変数が増えたら重回帰分析(Multiple Regression Analysis)と名前が変わります。 モデルの書き方はそれほど変わりません。
reg2 <- lm(pay~age+weight,data=b.data)
reg2##
## Call:
## lm(formula = pay ~ age + weight, data = b.data)
##
## Coefficients:
## (Intercept) age weight
## -13832.8 386.3 170.5summary(reg2)##
## Call:
## lm(formula = pay ~ age + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12895 -6013 -2714 1943 38730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13832.79 6997.60 -1.977 0.04946 *
## age 386.32 157.08 2.459 0.01478 *
## weight 170.48 61.11 2.790 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9809 on 197 degrees of freedom
## Multiple R-squared: 0.07236, Adjusted R-squared: 0.06294
## F-statistic: 7.683 on 2 and 197 DF, p-value: 0.0006125年齢が1年上だと年俸が386万円，体重が1kg高いと170万円高いようです。一応どちらも統計的にも有意であるようですが，どちらが重要なファクターなのでしょう？単位が違うから比べられないと思いませんか？
そこで，データの標準化を考えます。
b.data2 <- subset(b.data,select=c("pay","age","weight")) #データの一部抜き出し
b.data2z <- scale(b.data2) #標準化
b.data2z <- data.frame(b.data2z) #データフレーム型にしておく
summary(b.data2z)## pay age weight
## Min. :-0.9139 Min. :-2.3961 Min. :-1.8669
## 1st Qu.:-0.7363 1st Qu.:-0.5979 1st Qu.:-0.7302
## Median :-0.3662 Median :-0.1483 Median :-0.1180
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.2876 3rd Qu.: 0.7507 3rd Qu.: 0.5815
## Max. : 3.7048 Max. : 2.9984 Max. : 3.6420こうした上で，あらためて。
reg3 <- lm(pay~age+weight,b.data2z)
reg3##
## Call:
## lm(formula = pay ~ age + weight, data = b.data2z)
##
## Coefficients:
## (Intercept) age weight
## -1.093e-15 1.696e-01 1.924e-01summary(reg3)##
## Call:
## lm(formula = pay ~ age + weight, data = b.data2z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2726 -0.5934 -0.2679 0.1917 3.8222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.093e-15 6.845e-02 0.000 1.00000
## age 1.696e-01 6.897e-02 2.459 0.01478 *
## weight 1.924e-01 6.897e-02 2.790 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.968 on 197 degrees of freedom
## Multiple R-squared: 0.07236, Adjusted R-squared: 0.06294
## F-statistic: 7.683 on 2 and 197 DF, p-value: 0.0006125こうすると単位が整っていますから，大きさを比較することができます。年齢の係数は0.169,体重の係数は0.192なので，体重の方が年収に与える影響は大きいことがわかります。
うんちく
- 有意性検定，\(R^2\)は標準化したもの(reg3)としていないもの(reg2)とで変わりがないことを確認してください。
- 標準化したデータの係数(reg3)のことを標準化係数，していない生データの係数(reg2)のことを非標準化係数と言います。
- 説明変数が一つの回帰分析(reg1)の係数は回帰係数regression coefficientと言いますが，説明変数が二つ以上ある重回帰分析の係数は偏回帰係数partial regression coefficientと言います。
- partialとは，変数の影響を統計的に統制する作業のことを言います。
多重共線性の問題と変数選択
同じ説明変数を含んでいる回帰分析の場合，含める変数が多い方が常に良い適合度を示します。
reg4 <- lm(pay ~ age,data=b.data)
reg5 <- lm(pay ~ age+weight,b.data)
summary(reg4)##
## Call:
## lm(formula = pay ~ age, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13202 -6129 -3350 1941 36824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -165.2 5081.3 -0.033 0.97410
## age 430.4 158.9 2.708 0.00737 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9975 on 198 degrees of freedom
## Multiple R-squared: 0.03571, Adjusted R-squared: 0.03084
## F-statistic: 7.332 on 1 and 198 DF, p-value: 0.007367summary(reg5)##
## Call:
## lm(formula = pay ~ age + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12895 -6013 -2714 1943 38730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13832.79 6997.60 -1.977 0.04946 *
## age 386.32 157.08 2.459 0.01478 *
## weight 170.48 61.11 2.790 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9809 on 197 degrees of freedom
## Multiple R-squared: 0.07236, Adjusted R-squared: 0.06294
## F-statistic: 7.683 on 2 and 197 DF, p-value: 0.0006125こうなると何でもかんでも入れたくなりますが，これについて二つのトピックをお話しします。多重共線性の問題と変数選択です。
多重共線性multicollinearity
多重共線性 multicollinearity，通称マルチコは，説明変数同士の相関が高すぎて，推定値が適切な値にならない問題，を意味します。
例えば身長と体重で年俸を予測したとします。
reg6 <- lm(pay ~ height + weight,data=b.data)
summary(reg6)##
## Call:
## lm(formula = pay ~ height + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14611 -6389 -3336 3392 39219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26828.00 23219.79 1.155 0.24933
## height -205.41 155.72 -1.319 0.18868
## weight 271.76 89.69 3.030 0.00277 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9914 on 197 degrees of freedom
## Multiple R-squared: 0.05225, Adjusted R-squared: 0.04262
## F-statistic: 5.43 on 2 and 197 DF, p-value: 0.005065さて，身長と体重は高い相関を示すことがわかっています。
cor(b.data$height,b.data$weight)## [1] 0.7283599こういう時は回帰係数が適切に推定できていない可能性があります。 このエラーをチェックする指標がVIF(Variance Inflation Factor；分散拡大係数)です。この指標が2.0よりも大きければ，解釈に注意が必要です。
VIFはcarパッケージのvif関数で求められます。
# install.packages("car") #パッケージを持っていなければインストールする
library(car)## Warning: package 'car' was built under R version 3.2.2vif(reg6)## height weight
## 2.129962 2.129962変数選択
説明変数をいくらでも追加できるとしても，入れすぎたら「当てはまるけれども複雑すぎる」ことになり兼ねません。 統計的に有意な変数だけ入れて，適当に当てはまりの良いモデルを見繕ってくれたりしないかしら。
そういう願いに答えてくれるのが，変数選択法あるいはモデル選択と言います。 変数を選ぶ基準は，適合度です。適合度は\(R^2\)の他にもいろいろあります。ここで紹介する変数選択法の基準になっているのはAICという適合度基準で，「小さければ小さいほど当てはまりの良いモデル」ということを意味します。
reg.all <- lm(pay~year+age+height+weight,data=b.data)
reg.step <- step(reg.all) #変数選択## Start: AIC=3680.7
## pay ~ year + age + height + weight
##
## Df Sum of Sq RSS AIC
## - age 1 47462715 1.8748e+10 3679.2
## - height 1 55222462 1.8755e+10 3679.3
## - year 1 175833506 1.8876e+10 3680.6
## <none> 1.8700e+10 3680.7
## - weight 1 773686592 1.9474e+10 3686.8
##
## Step: AIC=3679.2
## pay ~ year + height + weight
##
## Df Sum of Sq RSS AIC
## - height 1 62199657 1.8810e+10 3677.9
## <none> 1.8748e+10 3679.2
## - year 1 616656059 1.9364e+10 3683.7
## - weight 1 1042136365 1.9790e+10 3688.0
##
## Step: AIC=3677.87
## pay ~ year + weight
##
## Df Sum of Sq RSS AIC
## <none> 1.8810e+10 3677.9
## - year 1 725482327 1.9535e+10 3683.4
## - weight 1 1379441303 2.0189e+10 3690.0summary(reg.step)##
## Call:
## lm(formula = pay ~ year + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13273 -6034 -2744 2390 38423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11408.89 6209.80 -1.837 0.067681 .
## year 382.88 138.90 2.756 0.006392 **
## weight 243.83 64.15 3.801 0.000192 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9771 on 197 degrees of freedom
## Multiple R-squared: 0.07938, Adjusted R-squared: 0.07004
## F-statistic: 8.493 on 2 and 197 DF, p-value: 0.0002896最終的に，年齢と体重で予測するモデルが最も当てはまりが良いことがわかりました(もっとも\(R^2\)はとても小さいものですが。)
本日の課題
- バッターのデータを使って，年俸をホームランの本数(HR)で予測する回帰分析を行ってください。
- 同じくバッターのデータで，年俸をホームランの本数(HR)，安打(hit)，二塁打(Double)，三塁打(Three)で予測する重回帰分析を行ってください。
- ホームランの本数をうまく説明する回帰分析をステップワイズ法で実施し，最適モデルをレポートしてください。