R 単回帰分析と重回帰分析をやってみる

Rを使って、単回帰分析と重回帰分析をやってみます。
サンプルデータは、WEB+DB PRESS vol.84で使われていたメジャーリーグの統計情報を使ってみます。
サンプルソースはこちら
https://github.com/hamuhamu/web_db_84_mlb

git clone git@github.com:hamuhamu/web_db_84_mlb.git

##### 環境構成 * R 3.1.1 <br> ##### 回帰分析とは？ 単回帰分析の説明に入る前に、相関係数の話をします。 相関係数は、２変数の関係性を表したもので、仮に正の相関が出ていても因果の向きはわかりません。 例えば、身長と体重に正の相関が出ていた場合、身長が高いことが原因で体重も重くなるのか体重が重いことが原因で身長も高くなるのか わかりません。 <br> 回帰分析は、因果の向きを仮定し１つの被説明変数を１つ又は複数の説明変数でどのくらい説明できるか評価するものです。 #### 被説明変数と説明変数とはなにか? 被説明変数: 　家賃 <br> 説明変数: 　駅からの距離 　部屋の広さ 　築年数 <br> 非説明変数である家賃がいくらかを説明するために説明変数を使って回帰式を作りどの説明変数がどのくらい家賃に影響を与えているかがわかります。 <br> 説明変数が１つなら単回帰、複数なら重回帰となります。 回帰式は、こんな１次方程式 単回帰　y = a + bx 重回帰　y = a + b1x1 + b2x2 + ・・・ <br> * y：非説明変数 * a：切片 * b：回帰係数 * x：説明変数 この回帰式により、yを予測することができます。 <br> 最後に決定係数と呼ばれるyの予測値の精度を評価します。 決定係数は、説明変数によって説明される分散の割合で、決定係数が1に近いほど精度は高いと言えます。 <br> ##### 単回帰分析 MLBのデータを使って単回帰分析してみます。 年度ごとの各チームの得点や打数などの通算の値が記録されています。 このデータを使い、何をするかというと、得点に対して打率と出塁率のどちらがより得点に寄与しているか調べます。 非説明変数： 　得点 説明変数： 　打率 　出塁率 <br> #### 打率と得点 まずは、打率と得点の散布図を出してみます。 <script type="syntaxhighlighter" class="brush: text" title="sample.R"> # 打率と得点 df <- read.table("4-1.table", header=T) head(df) ## 年度 リーグ 球団 得点 打数 安打 二塁打 三塁打 本塁打 打点 犠打 犠飛 死球 四球 故意四球 三振 盗塁 盗塁刺 併殺打 打撃妨害 打率 得点.1 ## 1 2013 AL ANA 733 5588 1476 270 39 164 696 37 64 48 523 40 1221 82 34 150 0 0.264 733 ## 2 2013 AL BAL 745 5620 1460 298 14 212 719 27 45 36 416 36 1125 79 29 105 0 0.260 745 ## 3 2013 AL BOS 853 5651 1566 363 29 178 819 24 50 72 581 51 1308 123 19 137 4 0.277 853 ## 4 2013 AL CHA 598 5563 1385 237 19 148 574 19 48 34 411 23 1207 105 42 124 2 0.249 598 ## 5 2013 AL CLE 745 5465 1391 290 23 171 712 31 56 51 562 19 1283 117 36 106 0 0.255 745 ## 6 2013 AL DET 796 5735 1625 292 23 176 767 32 47 43 531 40 1073 35 20 146 0 0.283 796 ## 打率と得点の相関関係 p <- ggplot(df, aes(x=打率, y=得点)) p <- p + geom_point() + theme_bw(base_family="Osaka") p <- p + geom_smooth(method="lm") p </script> <br> <a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi12U-05j0ZkkGGoLRbBkTMYMUXOArvsGXj9-d7fHVAb0YvT3Dm5rzSf-IPfO21cg_KRZUK73IPwQMUU0fosuvbJ_NVpZfTvCDJpCeqoQpKGDAVotNWmCUtRNJkZMB161djPGsi0aWHnSvx/s1600/%E6%89%93%E7%8E%87_%E5%BE%97%E7%82%B9.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi12U-05j0ZkkGGoLRbBkTMYMUXOArvsGXj9-d7fHVAb0YvT3Dm5rzSf-IPfO21cg_KRZUK73IPwQMUU0fosuvbJ_NVpZfTvCDJpCeqoQpKGDAVotNWmCUtRNJkZMB161djPGsi0aWHnSvx/s1600/%E6%89%93%E7%8E%87_%E5%BE%97%E7%82%B9.png" /></a> 次に決定係数を出します。 <script type="syntaxhighlighter" class="brush: text" title="sample.R"> ## 決定係数 df.reg <- lm(得点 ~ 出塁率, data=df) summary(df.reg) Call: lm(formula = 得点 ~ 打率, data = df) Residuals: Min 1Q Median 3Q Max -140.068 -33.592 -1.922 31.005 154.001 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -764.16 53.39 -14.31 <2e-16 *** 打率 5767.24 203.28 28.37 <2e-16 *** Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 49.05 on 418 degrees of freedom Multiple R-squared: 0.6582, Adjusted R-squared: 0.6574 F-statistic: 804.9 on 1 and 418 DF, p-value: < 2.2e-16 </script> <br> lm関数を使い回帰分析します。 summaryした結果を順番に見ていきます。 <br> #### Residuals 残差の四分位数 <pre> Residuals: Min 1Q Median 3Q Max -140.068 -33.592 -1.922 31.005 154.001 </pre> <br> #### Coefficients * Estimate：係数 * Std.Error：標準誤差 * t value：t値 * pr(>|t|)：p値 <pre> Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -764.16 53.39 -14.31 <2e-16 *** 打率 5767.24 203.28 28.37 <2e-16 *** Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 </pre> 5767.24が、回帰係数で-764.16が切片。 t値の絶対値が２以上であれば説明変数としてみなせる。 <br> #### Multiple R-squared 決定係数。 Multiple R-squaredの値が、１に近いほど精度が高い。 <pre> Multiple R-squared: 0.6582, Adjusted R-squared: 0.6574 </pre> <br> #### 回帰式 以上の結果から、回帰式はこうなる。 y(得点) = -764.16 + 打率 * 5767.24 この式がどのくらいの当てはまるか確認してみる。 2013年のANAとBALとBOSで確認した。 <br> | 球団 | 予測値 | 実測値 | | ANA | 758.3914 | 733 | | BAL | 735.3224 | 745 | | BOS | 833.3655 | 853 | <br> #### 出塁率と得点 次に出塁率と得点の散布図を出してみます。 <script type="syntaxhighlighter" class="brush: text" title="sample.R"> # 出塁率と得点 df <- read.table("4-2.table", header=T) head(df) ## 年度 リーグ 球団 得点 打数 安打 二塁打 三塁打 本塁打 打点 犠打 犠飛 死球 四球 故意四球 三振 盗塁 盗塁刺 併殺打 打撃妨害 出塁率 得点.1 ## 1 2013 AL ANA 733 5588 1476 270 39 164 696 37 64 48 523 40 1221 82 34 150 0 0.329 733 ## 2 2013 AL BAL 745 5620 1460 298 14 212 719 27 45 36 416 36 1125 79 29 105 0 0.313 745 ## 3 2013 AL BOS 853 5651 1566 363 29 178 819 24 50 72 581 51 1308 123 19 137 4 0.349 853 ## 4 2013 AL CHA 598 5563 1385 237 19 148 574 19 48 34 411 23 1207 105 42 124 2 0.302 598 ## 5 2013 AL CLE 745 5465 1391 290 23 171 712 31 56 51 562 19 1283 117 36 106 0 0.327 745 ## 6 2013 AL DET 796 5735 1625 292 23 176 767 32 47 43 531 40 1073 35 20 146 0 0.346 796 ## 出塁率と得点の相関関係 p <- ggplot(df, aes(x=出塁率, y=得点)) p <- p + geom_point() + theme_bw(base_family="Osaka") p <- p + geom_smooth(method="lm") p </script> <br> <a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhj0tS4R0F8q51Ae3yzvSF5cqcFGg7_P1YwARRJ-dKnO4pNSYXz3XudZim61huDfg1d6xml_bDOtzYhH6DjMRLLnWsA7J4yfcUxvihcKHWvveNXfbUmRdzZ782q0lw1o1fkmZoQe3fgNHcx/s1600/%E5%87%BA%E5%A1%81%E7%8E%87_%E6%89%93%E7%8E%87.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhj0tS4R0F8q51Ae3yzvSF5cqcFGg7_P1YwARRJ-dKnO4pNSYXz3XudZim61huDfg1d6xml_bDOtzYhH6DjMRLLnWsA7J4yfcUxvihcKHWvveNXfbUmRdzZ782q0lw1o1fkmZoQe3fgNHcx/s1600/%E5%87%BA%E5%A1%81%E7%8E%87_%E6%89%93%E7%8E%87.png" /></a> 次に決定係数を出します。 <script type="syntaxhighlighter" class="brush: text" title="sample.R"> ## 決定係数 df.reg <- lm(得点 ~ 打率, data=df) summary(df.reg) Call: lm(formula = 得点 ~ 出塁率, data = df) Residuals: Min 1Q Median 3Q Max -102.771 -23.821 -0.789 22.687 123.387 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1021.71 42.42 -24.09 <2e-16 *** 出塁率 5359.65 128.28 41.78 <2e-16 *** Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 36.88 on 418 degrees of freedom Multiple R-squared: 0.8068, Adjusted R-squared: 0.8064 F-statistic: 1746 on 1 and 418 DF, p-value: < 2.2e-16 </script> <br> 決定係数は、打率と得点より高い0.8068です。 #### 回帰式 以上の結果から、回帰式はこうなる。 y(得点) = -1021.71 + 出塁率 * 5359.65 この式がどのくらいの当てはまるか確認してみる。 2013年のANAとBALとBOSで確認した。 <br> | 球団 | 予測値 | 実測値 | | ANA | 741.6148 | 733 | | BAL | 655.8604 | 745 | | BOS | 848.8078 | 853 | <br> BALの精度は、下がってしまいましたがその他チームの精度は上がったようです。 <br> ##### 重回帰分析 単回帰分析は、説明変数一つで行いましたが重回帰分析は複数の説明変数で行います。 <br> 非説明変数： 　得点 説明変数： 　単打 　二塁打 　三塁打 　本塁打 　非故意四球 　故意四球 　盗塁 　盗塁刺 　凡打 　三振 　併殺打 　犠飛 　犠打 <br> <script type="syntaxhighlighter" class="brush: text" title="sample.R"> # 重回帰分析 df <- read.table("4-7.table", header=T) head(df) ## 年度 リーグ 球団 得点 単打 二塁打 三塁打 本塁打 非故意四球 故意四球 盗塁 盗塁刺 凡打 三振 併殺打 犠飛 犠打 ## 1 2013 AL ANA 733 1003 270 39 164 531 40 82 34 2891 1221 150 64 37 ## 2 2013 AL BAL 745 936 298 14 212 416 36 79 29 3035 1125 105 45 27 ## 3 2013 AL BOS 853 996 363 29 178 602 51 123 19 2777 1308 137 50 24 ## 4 2013 AL CHA 598 981 237 19 148 422 23 105 42 2971 1207 124 48 19 ## 5 2013 AL CLE 745 907 290 23 171 594 19 117 36 2791 1283 106 56 31 ## 6 2013 AL DET 796 1134 292 23 176 534 40 35 20 3037 1073 146 47 32 ## 重回帰に使用しない列を除去 df <- df[,c(-1, -2, -3)] head(df) ## 得点 単打 二塁打 三塁打 本塁打 非故意四球 故意四球 盗塁 盗塁刺 凡打 三振 併殺打 犠飛 犠打 ## 1 733 1003 270 39 164 531 40 82 34 2891 1221 150 64 37 ## 2 745 936 298 14 212 416 36 79 29 3035 1125 105 45 27 ## 3 853 996 363 29 178 602 51 123 19 2777 1308 137 50 24 ## 4 598 981 237 19 148 422 23 105 42 2971 1207 124 48 19 ## 5 745 907 290 23 171 594 19 117 36 2791 1283 106 56 31 ## 6 796 1134 292 23 176 534 40 35 20 3037 1073 146 47 32 ## 得点 ~ . 得点を得点以外の変数で説明する ## 決定係数 df.reg <- lm(得点 ~ ., data=df) summary(df.reg) Call: lm(formula = 得点 ~ ., data = df) Residuals: Min 1Q Median 3Q Max -62.982 -12.240 -0.844 13.684 75.920 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 347.40793 127.71200 2.720 0.00680 ** 単打 0.51247 0.02239 22.884 < 2e-16 *** 二塁打 0.70584 0.04432 15.927 < 2e-16 *** 三塁打 1.12746 0.12847 8.776 < 2e-16 *** 本塁打 1.40526 0.03717 37.807 < 2e-16 *** 非故意四球 0.34896 0.01862 18.742 < 2e-16 *** 故意四球 0.05836 0.08710 0.670 0.50326 盗塁 0.07171 0.04392 1.633 0.10327 盗塁刺 -0.15144 0.11806 -1.283 0.20033 凡打 -0.16903 0.02832 -5.968 5.24e-09 *** 三振 -0.20036 0.02870 -6.980 1.21e-11 *** 併殺打 -0.38957 0.08861 -4.396 1.41e-05 *** 犠飛 0.45688 0.15073 3.031 0.00259 ** 犠打 -0.17369 0.06657 -2.609 0.00941 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 20.88 on 406 degrees of freedom Multiple R-squared: 0.9398, Adjusted R-squared: 0.9379 F-statistic: 488 on 13 and 406 DF, p-value: < 2.2e-16 </script> <br> 決定係数は、出塁率と得点より高い0.9398です。 回帰式を書くのは大変なので、書きませんが単回帰分析の回帰式に説明変数の数が増えただけです。 回帰係数をみることで、何が一番得点に影響を与えているかわかります。 本塁打が1.40526と高く一番、得点に影響を与えていることがわかります。 <br> ##### 最近読んだ本 <iframe src="http://rcm-fe.amazon-adsystem.com/e/cm?t=hamuhamu06-22&o=9&p=8&l=as1&asins=4774169552&ref=qf_sp_asin_til&fc1=000000&IS2=1<1=_blank&m=amazon&lc1=0000FF&bc1=000000&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe> <iframe src="http://rcm-fe.amazon-adsystem.com/e/cm?t=hamuhamu06-22&o=9&p=8&l=as1&asins=4774170577&ref=qf_sp_asin_til&fc1=000000&IS2=1<1=_blank&m=amazon&lc1=0000FF&bc1=000000&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe> <iframe src="http://rcm-fe.amazon-adsystem.com/e/cm?t=hamuhamu06-22&o=9&p=8&l=as1&asins=4897979536&ref=qf_sp_asin_til&fc1=000000&IS2=1<1=_blank&m=amazon&lc1=0000FF&bc1=000000&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe>