Rを使って重回帰分析を実践　野球選手の年俸には何が影響しているのか？

前回は線形回帰について勉強しましたね！
randpy.hatenablog.com

今回は習った線形回帰を使って、実際にデータを使って分析をしてみます。

線形回帰というシンプルなモデルですが、色々な分野に応用できます。
今回は、野球選手の年俸が何によって影響を受けて決定されているのかを、Rを使って線形回帰で分析したいと思います！

もしRの実行環境がまだない場合は、以下のリンクを参考にしながら整えてください。
randpy.hatenablog.com
randpy.hatenablog.com

Webからデータを収集するためのライブラリ

データ分析するためには、当然ですがデータを集める必要があります。

今回は
http://baseball-data.comさんに掲載されている野球選手の2016年の成績と、2017年の推定年俸を用いることにしました。

一つ一つ情報をコピペするのは大変なので、プログラムで自動化することで一気にデータを作りましょう！
このような技術を「スクレイピング」と言います。

それぞれの言語にて、スクレイピングするのに特化したライブラリが存在します。Rでは「rvest」というライブラリが使いやすく、多くのユーザーに用いられています。

まずはライブラリをインストールしましょう。コンソール上で

>install.packeges("rvest")

と入力すればOKです。

もしライブラリをインストールするのが初めての場合、

パッケージを○○中にインストールします
(‘lib’ が指定されていないため)
--- このセッションで使うために、CRAN のミラーサイトを選んでください」

というメッセージのあと、インストールをどこからするか選択する画面が現れます。

基本的にRでは、CRAN (Comprehensive R Archive Network)というR自身やライブラリを管理しているサイトからインストールを行います。
CRANは全世界にミラーサイトがあるため、そのうちどこ経由（日本なのかイギリスなのか……など）でインストールを行うか選ぶ必要があります。

以前に「日本のミラーサイトからだと○○というライブラリでインストールエラーが出た……」という話を聞いたことがあるのですが、どのRのバージョンだったか、どのライブラリだったか全く覚えていません！すいません！（反省はあまりしていません！）

ですが、後からでもどのミラーサイトにするかは変更できるので、一旦は日本のミラーサイトを選択するので良いでしょう。

その後

>library(rvest)

このコマンドで先ほどインストールしたライブラリが使えるようになります。

スクレイピングするためのコードを紹介します

コードは以下の通りです。

github.com

基本このまま実行してもらえれば全て完了するのではと思うのですが、もしかするとお使いの環境によってエラーが出る可能性があります。

その場合はお手数ですが、コメントいただけると助かります。

ざっくりと分割しながらコードの解説をしてきます。
もし「コードうんぬんよりも分析結果を早く見たい！」という方はこの章は飛ばしていただき、結果を確認に進んでもらえればと思います。

スクレイピング部分

まずライブラリを読み込みます。

先ほど紹介した「rvest」以外にも必要なライブラリがあるので、インストールがまだの方は前述した方法にてそれぞれインストールをお願いします。

次にrvestを使って各球団の情報を読み込み、選手データのあるページ（年俸や成績が載っているページ）を一覧で取得します。

例えばRのコンソールで「urlListForBsH」と実行すると

>urlListForBsH 
 [1] "http://baseball-data.com/player/bs/%E4%B8%AD%E5%B3%B6%E3%80%80%E5%AE%8F%E4%B9%8B"         
 [2] "http://baseball-data.com/player/bs/%EF%BC%B4%EF%BC%8D%E5%B2%A1%E7%94%B0"                  
 [3] "http://baseball-data.com/player/bs/%E3%83%AD%E3%83%A1%E3%83%AD"                           
 [4] "http://baseball-data.com/player/bs/%E5%B0%8F%E8%B0%B7%E9%87%8E%E3%80%80%E6%A0%84%E4%B8%80"
 [5] "http://baseball-data.com/player/bs/%E5%AE%89%E9%81%94%E3%80%80%E4%BA%86%E4%B8%80"         
 [6] "http://baseball-data.com/player/bs/%E3%83%A2%E3%83%AC%E3%83%AB"                           
 [7] "http://baseball-data.com/player/bs/%E4%BC%8A%E8%97%A4%E3%80%80%E5%85%89"                  
 [8] "http://baseball-data.com/player/bs/%E3%83%9E%E3%83%AC%E3%83%BC%E3%83%AD"                  
 [9] "http://baseball-data.com/player/bs/%E8%A5%BF%E9%87%8E%E3%80%80%E7%9C%9F%E5%BC%98"         
[10] "http://baseball-data.com/player/bs/%E9%A7%BF%E5%A4%AA"                      
……

こんな感じで「urlListForBsH」にはオリックスの野手が、「urlListForBsP」にはオリックスの投手のURL一覧が格納されている状態です。

そしてそれを

>urlListAllH = c (urlListForCH, urlListForGH,urlListForYbH,urlListForTH,urlListForSH,urlListForDH, urlListForFH, urlListForHH,urlListForMH, urlListForEH, urlListForLH,urlListForBsH)
>urlListAllP = c (urlListForCP, urlListForGP,urlListForYbP,urlListForTP,urlListForSP,urlListForDP, urlListForFP, urlListForHP,urlListForMP, urlListForEP, urlListForLP,urlListForBsP)

こんな感じで野手の全てのデータ・投手の全てのデータという形でそれぞれまとめています。

次に先ほど取得した各選手のURLにいき、2017年の年俸と2016年の成績データを取得してきます。

>Sys.sleep(2)

for 文中に出てくるこちらですが、実行を2秒止めるというコードになります。
通常Webスクレイピングする際に、このような処理をしないと短時間に対象サイトに多くのアクセスをすることになり、サーバー負荷が大きくかかってしまいます。
よって、次のリクエストまで最低1秒以上間隔を空けることでサーバーへの負担を軽くすることが必要になります。

この処理は必ず入れるようにしてください。
データを頂く訳ですから、サイト運営者へ配慮をするようにしましょう。

次に，取得したデータにヘッダーをつけます。

これをしないと，「0.111」という数字が打率なのか？出塁率なのか？など分からなくなります。

エクセルでデータをまとめる際に1行目にデータを表すヘッダーを以下のようにいれますよね？
それも自動でやってしまおう！という処理です。

f:id:Np-Ur:20170621140859p:plain

ここで、どんなデータを所得できているのか確認してみましょう。Rのコンソールで「PplayerInformationAndRecordAll」と入力すると

> PplayerInformationAndRecordAll 
                                        名前       生年月日     投打 年齢 ポジション プロ通算年   出身地  身長 血液型  体重               年俸                                   経歴
  1:       42 ジョンソン(くりす　じょんそん) 1984年10月14日 左投左打 32歳       投手      3年目 アメリカ 193cm   不明  93kg 31,000万円（推定）   ウィチタ州大 - パイレーツ - ツインズ
  2:         19 野村　祐輔(のむら　ゆうすけ)  1989年6月24日 右投右打 27歳       投手      6年目     岡山 177cm   AB型  82kg 10,000万円（推定）                        広陵高 - 明治大
  3:       58 ジャクソン(じぇい　じゃくそん) 1987年10月27日 右投右打 29歳       投手      2年目 アメリカ 185cm   不明  88kg  9,600万円（推定）                ファーマン大 - パドレス
  4:       21 中﨑　翔太(なかざき　しょうた)  1992年8月10日 右投右打 24歳       投手      7年目   鹿児島 186cm    B型 100kg  8,500万円（推定）                             日南学園高
  5: 66 ヘーゲンズ(ぶれいでぃん　へーげんず)  1989年5月12日 右投右打 28歳       投手      2年目 アメリカ 190cm   不明  95kg  7,700万円（推定）      マーセド大 - ダイヤモンドバックス
 ---                                                                                                                                                                                 
253:               36 髙木　伴(たかぎ　ばん)   1990年6月1日 右投右打 27歳       投手      3年目     埼玉 180cm    B型  80kg    750万円（推定） 市立川口高 - 東京農業大 - ＮＴＴ東日本
254:         69 大山　暁史(おおやま　さとし)  1988年10月6日 左投左打 28歳       投手      4年目     大分 168cm    A型  74kg    590万円（推定）     別府青山高 - 亜細亜大 - セガサミー
255:               68 鈴木　優(すずき　ゆう)   1997年2月5日 右投右打 20歳       投手      3年目     東京 181cm    B型  79kg    500万円（推定）                                 雪谷高
256:         48 齋藤　綱記(さいとう　こうき) 1996年12月18日 左投左打 20歳       投手      3年目   北海道 182cm    O型  89kg    500万円（推定）                                 北照高
257:        123 角屋　龍太(かどや　りゅうた) 1990年10月18日 右投右打 26歳       投手      2年目     岐阜 175cm    A型  75kg    440万円（推定）   富田高 - 名城大 - ジェイプロジェクト
        ドラフト                         タイトル 年度       球団 試合 勝利 敗北 セlブ 完投 完封勝 無四球 打者 投球回 被安打 被本塁打 与四球 与死球 奪三振 暴投 ボlク 失点 自責点 防御率 WHIP
  1:    <U+00A0>                 （防）15（沢）16 2016   広島東洋   26   15    7     0    3      2      0  736  180.1    154       11     49      3    141    3     0   50     43   2.15 1.13
  2: ドラフト1位 （率）16（勝）16（ベ）16（新）12 2016   広島東洋   25   16    3     0    1      1      0  633  152.2    139       11     37      6     91    3     0   50     46   2.71 1.15
  3:    <U+00A0>                                - 2016   広島東洋   67    5    4     0    0      0      0  271   68.1     46        4     23      0     89    1     0   15     13   1.71 1.01
  4: ドラフト6位                                - 2016   広島東洋   61    3    4    34    0      0      0  251   61.1     50        2     19      2     54    0     0   11      9   1.32 1.13
  5:    <U+00A0>                                - 2016   広島東洋   50    7    5     0    0      0      0  342   83.1     71        4     33      3     33    4     0   32     27   2.92 1.25
 ---                                                                                                                                                                                         
253: ドラフト4位                                - 2016 オリックス    4    0    0     0    0      0      0   25    6.0      5        2      2      0      5    0     1    4      4   6.00 1.17
254: ドラフト8位                                - 2016 オリックス    3    0    0     0    0      0      0   27    4.1      6        0      5      2      2    0     0    8      1   2.08 2.54
255: ドラフト9位                                - 2016 オリックス    1    0    0     0    0      0      0   12    1.1      5        0      3      0      3    0     0    5      5  33.75 6.00
256: ドラフト5位                                - 2016 オリックス    1    0    0     0    0      0      0   20    4.0      6        2      2      0      2    0     1    4      4   9.00 2.00
257: ドラフト8位                                - 2016 オリックス    2    0    0     0    0      0      0   16    2.1      5        1      4      0      1    0     0    3      3  11.57 3.86

このように選手の情報を取得することに成功しました！全部で257名の選手データです。

なお，取得したデータを保存しておきたい場合は以下を実行してcsvファイル形式で保存しておきましょう。

>write.csv(PplayerInformationAndRecordAll, "PplayerInformationAndRecordAll.csv")
>write.csv(HplayerInformationAndRecordAll, "HplayerInformationAndRecordAll.csv")

もし次にRを立ち上げた際にそのデータを使いたい場合は、次のコードを実行しましょう。

>HplayerInformationAndRecordAll = fread("HplayerInformationAndRecordAll.csv", encoding="UTF-8")
>PplayerInformationAndRecordAll = fread("PplayerInformationAndRecordAll.csv", encoding="UTF-8")

取得したデータの整備

次に先ほど取得したデータのうち分析に用いるデータを選択します。例えば、「生年月日」もスクレイピングして取得してあるのですが、正直生年月日が年俸と関係しているとは思えないので外します。

この辺りは何のデータを用いるか、いくつかパターンを変えて実行してみるのも面白いかもしれません。

さて、データを選び終えましたが、このままだと年俸の列に「○○万円」、年齢の列に「○○歳」と邪魔な文字列が入っています。
これを以下の処理で取り除きます。

これで分析ができそうなのですが……最後にもうひと手間必要です。
今回スクレイピングしたデータは基本的にRが数字ではなく文字列として認識してしまっています。

そこで、数値として扱いたいデータは数値に変換してあげる必要があります。
それが以下の部分です。

正直もっと綺麗にかける方法があるはずなのですが、こんな書き方で処理をしています。

まあ全体的に分かりやすさを重視して冗長なコードになっているので……
という言い訳をしておきます。

いよいよ線形回帰部分

最後に何を目的変数とするのか？何を説明変数とするのか？を選びます。
ここからは投手のデータを用いて分析することを考えます。（野手のデータの分析はみなさんやってみてください！）

ここでは、年俸を目的変数、年俸以外の変数を説明変数としています。
復習ですが、線形回帰をする場合、目的変数に正規分布を仮定していることを心に留めておいてください。

目的変数・説明変数については，以下を参照してください。
randpy.hatenablog.com

そしてその変数を用いて線形回帰を行います。
Rで線形回帰を行う関数は「lm」です。

このように

lm(目的変数 ~ 説明変数)

とするだけで線形回帰を実行してくれます。

結果を確認！

さて、長かったデータ処理も終わり、ついに線形回帰まで成功しました！

結果を確認するには「summary」という関数を用います。

>summary(resultP)

コンソール上で実行すると以下のような結果が出力されます。

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -21797.44    3148.70  -6.923 4.58e-11 ***
Xp年齢            825.56      90.75   9.097  < 2e-16 ***
Xp試合            33.83      28.55   1.185  0.23729    
Xp勝利            252.85     290.07   0.872  0.38431    
Xp敗北            -131.26     294.06  -0.446  0.65575    
Xpセlブ           375.93      72.68   5.173 5.10e-07 ***
Xp完投            1268.89     961.55   1.320  0.18830    
Xp完封勝          -153.32    1625.00  -0.094  0.92491    
Xp無四球          -373.31    1813.93  -0.206  0.83713    
Xp打者            -39.76     139.91  -0.284  0.77651    
Xp投球回          246.42     418.78   0.588  0.55684    
Xp被安打          -145.77     148.66  -0.981  0.32785    
Xp被本塁打        -64.76     201.50  -0.321  0.74820    
Xp与四球          -141.94     140.13  -1.013  0.31221    
Xp与死球          -331.74     257.56  -1.288  0.19907    
Xp奪三振          49.23      36.37   1.353  0.17731    
Xp暴投            -306.60     269.92  -1.136  0.25720    
Xpボlク           -1330.20     949.98  -1.400  0.16282    
Xp失点            611.53     209.35   2.921  0.00384 ** 
Xp自責点          -301.51     203.60  -1.481  0.14004    
Xp防御率          -63.62     180.79  -0.352  0.72523    
XpWHIP            1416.43    1277.07   1.109  0.26856    
横浜DeNA          -4461.17    1790.02  -2.492  0.01341 *  
広島東洋          -1675.88    1795.38  -0.933  0.35159    
阪神              -1867.45    1722.30  -1.084  0.27940    
埼玉西武          -1822.90    1724.58  -1.057  0.29164    
千葉ロッテ        -2588.56    1761.02  -1.470  0.14298    
東京ヤクルト      -2915.29    1772.39  -1.645  0.10140    
東北楽天          -2481.85    1796.36  -1.382  0.16847    
読売                737.68    1773.36   0.416  0.67782    
福岡ソフトバンク   4617.77    1803.17   2.561  0.01109 *  
北海道日本ハム     -962.28    1801.25  -0.534  0.59371    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5732 on 225 degrees of freedom
Multiple R-squared:  0.6211,    Adjusted R-squared:  0.5689 
F-statistic:  11.9 on 31 and 225 DF,  p-value: < 2.2e-16

左から

説明変数
推定された係数（1単位上がるとどれほど年俸に影響するか）
標準偏差（どれほど推定された値にばたつきがあるか）
t値（推定された値が有意かどうか、検定の際に用いるもの）
p値（推定された係数が0である確率）
有意性

を表しています。
推定された係数の見方として、例えば「セーブ」の変数では「375.93」という結果になっており、1セーブ増えごとに年棒は375万円増えるという意味になります。
また「敗北」の変数では「-131」という結果になっているので、1回負けるごとに年俸は131万円ずつ減る、という意味になります。

ただし、この推定された結果は偶然この結果になっただけで、もっと大きなデータを用いたり、用いるデータを変更すると結果が変わる可能性があります。
そこで、通常分析した際はその推定された係数が「実際は0（=年俸に影響はない）ではないか？」と検証する必要があります。

そこで必ず「p値」を確認しましょう！
例えば「年齢」の係数は「825」と一歳年を重ねるごとに825万円上がるという結果が出ており、p値は $2e-16 \fallingdotseq 0$ となっています。これは偶然この結果が得られる確率がとても低いということを意味しているため、「やっぱり偶然ではなく相関ありそう！」と結論付けられるという訳です。

ちなみに「じゃあp値はどのくらい低ければいいの？」という疑問が湧きそうですが、この基準として有意水準0.05,0.01というのが一般的に定められていて、この有意確率よりも低ければある程度信頼していいでしょう、と結論づけている論文が多いです。

一番右のマークはp値がどれほど小さいかを意味していて、

*** → p値が0.001以下
** → p値が0.01以下
* → p値が0.05以下
. → p値が0.1以下

となっています。それぞれの変数のp値を一つずつ確認するのは面倒な場合は、このマークをパッと確認することで有意な結果が得られた変数を見つけることができます。

最後に決定係数ですが、

Multiple R-squared:  0.6211,    Adjusted R-squared:  0.5689

と出力されています。つまり、今回の説明変数によって年俸の62%ぐらいの影響は推定できている、ということになります。

残差の分析

次に残差（実際の値と推定された値との差）を見てみます。

以下を実行することで、各選手と残差の一覧を取得できます。

さらに237行目で残差の昇順に並び替えを行っています。

取得できたら、コンソールで「orderedResidualsPWithName」と入力して結果をみてみましょう。

> orderedResidualsPWithName 
                                          名前 residualsP
  1:             16 東浜　巨(ひがしはま　なお) -13711.716
  2:           41 千賀　滉大(せんが　こうだい) -13002.159
  3:           66 柳瀬　明宏(やなせ　あきひろ) -11549.272
  4:           17 岩貞　祐太(いわさだ　ゆうた) -10035.113
  5:               21 武田　久(たけだ　ひさし)  -9976.915
 ---                                                     
253: 44 バンデンハーク(りっく　ばんでんはーく)  16948.043
254:       53 五十嵐　亮太(いがらし　りょうた)  19089.408
255:         18 松坂　大輔(まつざか　だいすけ)  20778.298
256:               50 攝津　正(せっつ　ただし)  24142.063
257:             19 金子　千尋(かねこ　ちひろ)  30679.20

このマイナスの値というのは、推定された年俸よりも実際の年俸が低い（本当はもっと貰っても良いのでは？）ということを意味しており、プラスの値は推定された年俸よりも実際の年俸が高い（成績のわりに貰いすぎでは？）ということを意味しています。

WBCで活躍され千賀選手などは残差がかなり低く、メジャー帰りの松坂選手は残差がかなり高いです。
この辺りは何となく直観的に合っていそうですね！

モデルを改良するには？

今回、年俸と目的変数、2016年の成績を説明変数にして線形回帰を行いましたが、あまり有意な結果は得られませんでした。

なのでモデルを改良する必要がありそうです。

パッと思いつく改良点としては、

実際には去年の成績のみが影響を受ける訳ではないので、2016年以前の成績も変数に加える
2016年・2015年など過去の年俸から比較して今年の年俸が出ていそうなので、変数に加える
年俸は果たして正規分布なのか？という疑問があるので別の確率分布を使用する

などがあります。
最後の項目に関しては、一般的に年俸は正規分布ではなく「対数正規分布」というのがフィットすると言われています。

そこで今回の目的である線形回帰からはトピックが外れますが、対数正規分布を仮定して分析した結果も最後に提示してみます。ちなみに正規分布以外の確率分布を柔軟に設定して分析する手法を、一般化線形モデルと呼びます。

一般化線形モデルについては、また別の機会（2週間後ぐらい？）に執筆するとしますが、とりあえずコードだけ紹介しておきます。

>resultPglm = glm(Yp ~ Xp+as.factor(dataForLinearRegressionP[,1]),family = gaussian(link = "log"))

これだけです。とても簡単ですねー。
分析結果を出力すると、

>summary(resultPglm)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.837291   0.456514   6.215 2.46e-09 ***
Xp年齢           0.166038   0.012075  13.751  < 2e-16 ***
Xp試合            0.016152   0.002947   5.482 1.13e-07 ***
Xp勝利           0.082096   0.026043   3.152 0.001840 ** 
Xp敗北           -0.040338   0.033686  -1.197 0.232389    
Xpセlブ          0.007125   0.005712   1.247 0.213553    
Xp完投           0.294978   0.066170   4.458 1.31e-05 ***
Xp完封勝         -0.056938   0.120703  -0.472 0.637581    
Xp無四球         0.225408   0.127353   1.770 0.078091 .  
Xp打者           -0.004413   0.012200  -0.362 0.717928    
Xp投球回         0.031268   0.037063   0.844 0.399761    
Xp被安打         -0.023875   0.012903  -1.850 0.065579 .  
Xp被本塁打       -0.084251   0.016032  -5.255 3.43e-07 ***
Xp与四球         -0.026424   0.012608  -2.096 0.037225 *  
Xp与死球         -0.064705   0.022068  -2.932 0.003714 ** 
Xp奪三振         -0.003705   0.002827  -1.310 0.191405    
Xp暴投           0.088690   0.025666   3.456 0.000656 ***
Xpボlク          -0.241395   0.127326  -1.896 0.059258 .  
Xp失点           0.120867   0.021176   5.708 3.59e-08 ***
Xp自責点         -0.053907   0.019227  -2.804 0.005492 ** 
Xp防御率         -0.028942   0.040935  -0.707 0.480283    
XpWHIP           0.249171   0.169645   1.469 0.143289    
横浜DeNA         -1.209370   0.304512  -3.971 9.61e-05 ***
広島東洋         -0.125393   0.211011  -0.594 0.552943    
阪神             -0.661155   0.190470  -3.471 0.000621 ***
埼玉西武         -0.274208   0.224541  -1.221 0.223291    
千葉ロッテ       -0.670380   0.240984  -2.782 0.005864 ** 
東京ヤクルト     -0.467428   0.213405  -2.190 0.029527 *  
東北楽天         -1.120519   0.347897  -3.221 0.001467 ** 
読売              0.020838   0.167500   0.124 0.901106    
福岡ソフトバンク  0.642199   0.155514   4.130 5.13e-05 ***
北海道日本ハム   -0.123035   0.181312  -0.679 0.498100