【因子分析, R言語, スクレイピング】プロ野球選手の打撃成績を「長打力」「打撃技術」「走力」に要約すると上位の選手は誰になるのか！？
因子分析の練習としてどのような分析を行おうかと考えていたところ、以下のブログに面白そうな記事を発見しました。
nijyester.blog.fc2.com
nijyester.blog.fc2.com
この記事では過去10年間におけるパリーグの打者（規定打席到達に限る）の打撃成績から「長打力」「打撃技術」「走力」といった3つの共通因子を抽出し、それぞれの因子得点が上位の打者を紹介していました。残念ながらコードが記載されていませんでしたので、本記事では2011年から20016年の6年間におけるセパ両リーグの打者（規定打席の1/2以上）の打撃成績にデータを拡張して上記ブログの分析を再現してみます。
使用するパッケージ
library(rvest) library(dplyr) library(stringr) library(magrittr) library(ggplot2) library(psych) library(GPArotation)
データセットの作成
分析に利用するデータは、プロ野球データFreakさん（http://baseball-data.com/）のウェブサイトからスクレイピングして取得します。
スクレイピングのコードは以下のとおりです。
> session <- html_session("http://baseball-data.com/") # 2016年の12球団の打撃成績データの取得 > carp16 <- session %>% jump_to("16/stats/hitter2-c/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > giants16 <- session %>% jump_to("16/stats/hitter2-g/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > dena16 <- session %>% jump_to("16/stats/hitter2-yb/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > tigers16 <- session %>% jump_to("16/stats/hitter2-t/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > swallows16 <- session %>% jump_to("16/stats/hitter2-s/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > dragons16 <- session %>% jump_to("16/stats/hitter2-d/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > fighters16 <- session %>% jump_to("16/stats/hitter2-f/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > hawks16 <- session %>% jump_to("16/stats/hitter2-h/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > marines16 <- session %>% jump_to("16/stats/hitter2-m/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > lions16 <- session %>% jump_to("16/stats/hitter2-l/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > eagles16 <- session %>% jump_to("16/stats/hitter2-e/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > buffaloes16 <- session %>% jump_to("16/stats/hitter2-bs/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) # 2015年の12球団の打撃成績データの取得 > carp15 <- session %>% jump_to("15/stats/hitter2-c/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > giants15 <- session %>% jump_to("15/stats/hitter2-g/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > dena15 <- session %>% jump_to("15/stats/hitter2-yb/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > tigers15 <- session %>% jump_to("15/stats/hitter2-t/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > swallows15 <- session %>% jump_to("15/stats/hitter2-s/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > dragons15 <- session %>% jump_to("15/stats/hitter2-d/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > fighters15 <- session %>% jump_to("15/stats/hitter2-f/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > hawks15 <- session %>% jump_to("15/stats/hitter2-h/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > marines15 <- session %>% jump_to("15/stats/hitter2-m/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > lions15 <- session %>% jump_to("15/stats/hitter2-l/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > eagles15 <- session %>% jump_to("15/stats/hitter2-e/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) > buffaloes15 <- session %>% jump_to("15/stats/hitter2-bs/tpa-4.html") %>% html_nodes(xpath = '//*[@id="tbl"]') %>% > html_table() %>% extract2(1) #2011年から2014年の上記コードは省略
次に取得したデータに球団名と年度のカラムを追加するため以下のようにします。
#球団名カラムの追加 #2016年 carp16$team <- "広島" giants16$team <- "巨人" dena16$team <- "横浜" tigers16$team <- "阪神" swallows16$team <- "ヤクルト" dragons16$team <- "中日" fighters16$team <- "日本ハム" hawks16$team <- "ソフトバンク" marines16$team <- "ロッテ" lions16$team <- "西武" eagles16$team <- "楽天" buffaloes16$team <- "オリックス" #2015年 carp15$team <- "広島" giants15$team <- "巨人" dena15$team <- "横浜" tigers15$team <- "阪神" swallows15$team <- "ヤクルト" dragons15$team <- "中日" fighters15$team <- "日本ハム" hawks15$team <- "ソフトバンク" marines15$team <- "ロッテ" lions15$team <- "西武" eagles15$team <- "楽天" buffaloes15$team <- "オリックス" #2011年から2014年は省略 #年度カラムを追加 #2016年 carp16$year <- 2016 giants16$year <- 2016 dena16$year <- 2016 tigers16$year <- 2016 swallows16$year <- 2016 dragons16$year <- 2016 fighters16$year <- 2016 hawks16$year <- 2016 marines16$year <- 2016 lions16$year <- 2016 eagles16$year <- 2016 buffaloes16$year <- 2016 #2015年 carp15$year <- 2015 giants15$year <- 2015 dena15$year <- 2015 tigers15$year <- 2015 swallows15$year <- 2015 dragons15$year <- 2015 fighters15$year <- 2015 hawks15$year <- 2015 marines15$year <- 2015 lions15$year <- 2015 eagles15$year <- 2015 buffaloes15$year <- 2015 #2011年から2014年は省略
次に各年度のデータフレームを結合します。
#2011年から2016年の各球団の打撃成績データフレームを行結合する > hitter <- rbind(carp16,giants16,dena16,tigers16,swallows16,dragons16,fighters16,hawks16,marines16,lions16,eagles16,buffaloes16, > carp15,giants15,dena15,tigers15,swallows15,dragons15,fighters15,hawks15,marines15,lions15,eagles15,buffaloes15, > carp14,giants14,dena14,tigers14,swallows14,dragons14,fighters14,hawks14,marines14,lions14,eagles14,buffaloes14, > carp13,giants13,dena13,tigers13,swallows13,dragons13,fighters13,hawks13,marines13,lions13,eagles13,buffaloes13, > carp12,giants12,dena12,tigers12,swallows12,dragons12,fighters12,hawks12,marines12,lions12,eagles12,buffaloes12, > carp11,giants11,dena11,tigers11,swallows11,dragons11,fighters11,hawks11,marines11,lions11,eagles11,buffaloes11) > glimpse(hitter) Observations: 726 Variables: 26 $ 背番号 <chr> "2", "9", "33", "51", "25", "55", "60", "31", "44", "5", "背番号",... $ 選手名 <chr> "田中 広輔", "丸 佳浩", "菊池 涼介", "鈴木 誠也", "新井 貴浩", "エルドレッド", "安部 友裕", ... $ 打率 <chr> ".265", ".291", ".315", ".335", ".300", ".294", ".282", ".202"... $ 試合 <chr> "143", "143", "141", "129", "132", "95", "115", "106", "103", ... $ 打席数 <chr> "679", "652", "640", "528", "513", "354", "292", "289", "275", ... $ 打数 <chr> "581", "557", "574", "466", "454", "316", "259", "243", "254",... $ 得点 <chr> "102", "98", "92", "76", "66", "42", "23", "19", "34", "35", "... $ 安打 <chr> "154", "162", "181", "156", "136", "93", "73", "49", "74", "66... $ 二塁打 <chr> "17", "30", "22", "26", "23", "14", "12", "7", "10", "8", "二塁打"... $ 三塁打 <chr> "3", "8", "3", "8", "2", "0", "4", "0", "2", "1", "三塁打", "4", "... $ 本塁打 <chr> "13", "20", "13", "29", "19", "21", "6", "0", "10", "5", "本塁打",... $ 塁打 <chr> "216", "268", "248", "285", "220", "170", "111", "56", "118", ... $ 打点 <chr> "39", "90", "56", "95", "101", "53", "33", "17", "41", "34", "... $ 盗塁 <chr> "28", "23", "13", "16", "0", "1", "7", "4", "0", "6", "盗塁", "8... $ 盗塁刺 <chr> "19", "9", "5", "11", "1", "0", "0", "1", "0", "0", "盗塁刺", "2",... $ 犠打 <chr> "3", "1", "23", "3", "0", "0", "7", "12", "0", "0", "犠打", "1",... $ 犠飛 <chr> "1", "3", "3", "3", "4", "3", "5", "1", "1", "2", "犠飛", "3", "... $ 四球 <chr> "77", "84", "40", "53", "54", "31", "19", "29", "20", "22", "四... $ 敬遠 <chr> "1", "1", "0", "1", "1", "0", "2", "5", "1", "0", "敬遠", "0", "... $ 死球 <chr> "17", "7", "0", "3", "1", "4", "1", "4", "0", "1", "死球", "5", ... $ 三振 <chr> "119", "107", "106", "79", "101", "86", "64", "54", "29", "47"... $ 併殺打 <chr> "1", "9", "3", "10", "12", "8", "3", "8", "4", "5", "併殺打", "7",... $ 長打率 <chr> ".372", ".481", ".432", ".612", ".485", ".538", ".429", ".230",... $ 出塁率 <chr> ".367", ".389", ".358", ".404", ".372", ".362", ".327", ".296",... $ team <chr> "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", ... $ year <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, ...
データフレームhitterの各カラムの型がchr型になっているので、これをnumeric型に変換します。
> hitter %<>% mutate_each(funs(as.numeric),打率,試合,打席数,打数,得点,安打,二塁打,三塁打,本塁打, > 塁打,打点,盗塁,盗塁刺,犠打,犠飛,四球,敬遠,死球,三振,併殺打,長打率,出塁率) > glimpse(hitter) Observations: 726 Variables: 26 $ 背番号 <chr> "2", "9", "33", "51", "25", "55", "60", "31", "44", "5", "背番号",... $ 選手名 <chr> "田中 広輔", "丸 佳浩", "菊池 涼介", "鈴木 誠也", "新井 貴浩", "エルドレッド", "安部 友裕", ... $ 打率 <dbl> 0.265, 0.291, 0.315, 0.335, 0.300, 0.294, 0.282, 0.202, 0.291,... $ 試合 <dbl> 143, 143, 141, 129, 132, 95, 115, 106, 103, 67, NA, 143, 143, ... $ 打席数 <dbl> 679, 652, 640, 528, 513, 354, 292, 289, 275, 268, NA, 618, 576,... $ 打数 <dbl> 581, 557, 574, 466, 454, 316, 259, 243, 254, 243, NA, 576, 529... $ 得点 <dbl> 102, 98, 92, 76, 66, 42, 23, 19, 34, 35, NA, 58, 58, 96, 48, 2... $ 安打 <dbl> 154, 162, 181, 156, 136, 93, 73, 49, 74, 66, NA, 163, 160, 168... $ 二塁打 <dbl> 17, 30, 22, 26, 23, 14, 12, 7, 10, 8, NA, 28, 32, 28, 22, 12, 1... $ 三塁打 <dbl> 3, 8, 3, 8, 2, 0, 4, 0, 2, 1, NA, 4, 0, 3, 1, 1, 0, 0, 1, 0, NA... $ 本塁打 <dbl> 13, 20, 13, 29, 19, 21, 6, 0, 10, 5, NA, 11, 25, 23, 24, 4, 12,... $ 塁打 <dbl> 216, 268, 248, 285, 220, 170, 111, 56, 118, 91, NA, 232, 267, ... $ 打点 <dbl> 39, 90, 56, 95, 101, 53, 33, 17, 41, 34, NA, 42, 81, 75, 68, 3... $ 盗塁 <dbl> 28, 23, 13, 16, 0, 1, 7, 4, 0, 6, NA, 8, 1, 13, 0, 2, 0, 0, 7,... $ 盗塁刺 <dbl> 19, 9, 5, 11, 1, 0, 0, 1, 0, 0, NA, 2, 3, 3, 1, 1, 1, 2, 3, 2, ... $ 犠打 <dbl> 3, 1, 23, 3, 0, 0, 7, 12, 0, 0, NA, 1, 2, 1, 0, 19, 1, 0, 20, ... $ 犠飛 <dbl> 1, 3, 3, 3, 4, 3, 5, 1, 1, 2, NA, 3, 2, 6, 0, 1, 3, 3, 0, 0, N... $ 四球 <dbl> 77, 84, 40, 53, 54, 31, 19, 29, 20, 22, NA, 33, 38, 81, 39, 36... $ 敬遠 <dbl> 1, 1, 0, 1, 1, 0, 2, 5, 1, 0, NA, 0, 0, 2, 0, 5, 1, 0, 0, 0, N... $ 死球 <dbl> 17, 7, 0, 3, 1, 4, 1, 4, 0, 1, NA, 5, 5, 0, 4, 4, 4, 2, 0, 0, ... $ 三振 <dbl> 119, 107, 106, 79, 101, 86, 64, 54, 29, 47, NA, 78, 83, 67, 10... $ 併殺打 <dbl> 1, 9, 3, 10, 12, 8, 3, 8, 4, 5, NA, 7, 21, 6, 10, 10, 8, 7, 4, ... $ 長打率 <dbl> 0.372, 0.481, 0.432, 0.612, 0.485, 0.538, 0.429, 0.230, 0.465, ... $ 出塁率 <dbl> 0.367, 0.389, 0.358, 0.404, 0.372, 0.362, 0.327, 0.296, 0.342, ... $ team <chr> "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", "広島", ... $ year <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, ...
ここで、分析に利用する変数を以下のように選択します。
> hitter2 <- hitter %>% dplyr::select(year,team,選手名,打率,安打,二塁打,三塁打,本塁打,打点,得点,盗塁,盗塁刺, > 犠打,犠飛,四球,死球,併殺打,三振)
これらの変数がすべて欠測していないサンプルを利用して因子分析を行います。
> hitter3 <- na.omit(hitter2[,1:18]) > dim(hitter3) 654 18
本分析で利用するすべての変数が欠測していないサンプルサイズは654となります。
平行分析による因子数の決定
因子数の決定方法には様々な方法がありますが、ここでは平行分析を利用します。ざっくりとした説明をすると、乱数から作成した相関行列の固有値と分析に利用する変数の相関行列の固有値のプロットをみて因子数を決定する方法を平行分析といいます。Rではpsychパッケージのfa.parallel()関数を使用して実行ができます。
> fa.parallel(hitter3[,4:18])
因子分析
psych パッケージ のfa()関数を使用する。因子数を3とし、初期値解の算出に最尤法を使用し、因子の回転バリマックス法（直行回転法）を適用します。
> result <- fa(hitter3[,4:18],nfactors=3,rotate="varimax",fm="ml",scores=T) > print(result,sort=T) Factor Analysis using method = ml Call: fa(r = hitter3[, 4:18], nfactors = 3, rotate = "varimax", scores = T, fm = "ml") Standardized loadings (pattern matrix) based upon correlation matrix item ML1 ML2 ML3 h2 u2 com 安打 2 0.88 0.30 0.36 1.00 0.005 1.6 二塁打 3 0.72 0.40 0.11 0.70 0.304 1.6 打率 1 0.62 0.15 0.19 0.44 0.563 1.3 得点 7 0.60 0.53 0.49 0.89 0.106 2.9 併殺打 14 0.51 0.30 -0.26 0.42 0.581 2.2 犠飛 11 0.46 0.28 -0.16 0.31 0.686 1.9 本塁打 5 0.29 0.88 -0.20 0.90 0.097 1.3 打点 6 0.59 0.76 -0.11 0.94 0.062 1.9 三振 15 0.28 0.65 0.17 0.53 0.473 1.5 四球 12 0.44 0.54 0.25 0.54 0.456 2.4 犠打 10 -0.12 -0.36 0.35 0.27 0.731 2.2 死球 13 0.25 0.26 0.12 0.14 0.858 2.4 盗塁 8 0.02 0.08 0.86 0.75 0.249 1.0 盗塁刺 9 0.05 0.02 0.83 0.70 0.303 1.0 三塁打 4 0.20 -0.03 0.57 0.36 0.639 1.2 ML1 ML2 ML3 SS loadings 3.35 2.99 2.55 Proportion Var 0.22 0.20 0.17 Cumulative Var 0.22 0.42 0.59 Proportion Explained 0.38 0.34 0.29 Cumulative Proportion 0.38 0.71 1.00 （以下略）
print()関数の出力結果を説明しておきます。表頭のML1、ML2、ML3は共通因子名を表しており、その要素に因子負荷量を表示しています。次にh2、u2、comは共通性、独自性、複雑性を表しています。ML1の因子負荷量をみると安打、二塁打、打率といった打撃成績が高い値を有しています。そのためML1は「打撃技術」として解釈しても問題なさそうです。次に、ML2の因子負荷量をみると本塁打、打点、三振が高い値を有しており「長打力」として解釈できそうです。最後にML3は盗塁、盗塁刺が高い値を有しており「走力」示す因子として解釈できそうです。共通因子ごとの因子負荷量を視覚化すると以下のようになります。
> par(mfrow = c(3, 1)) > barplot(result$loadings[,1], main="ML1", ylim=c(0,1)) > barplot(result$loadings[,2], main="ML2", ylim=c(0,1)) > barplot(result$loadings[,3], main="ML3", ylim=c(0,1))
「長打力」「打撃技術」「走力」別の因子得点が上位の選手
さて因子得点を算出したとしても、選手名や成績の年度、球団名といったカラムがありません。なので、色々方法があると思うのですが、ここではマトリックス形式の因子得点をデータフレーム形式にして元のデータフレームと結合(結合キーも作成する)します。そして結合後のデータフレームから年度、チーム名、選手名、因子得点のカラムを抽出して、長打力」「打撃技術」「走力」それぞれの因子得点が上位の選手を表示したいと思います。
# データフレーム形式にする > fc_score <- as.data.frame(result$scores) # 結合キーの作成 > fc_score$id <- 1:654 > hitter3$id <- 1:654 # hitter3と内部結合 > con_vm <- dplyr::inner_join(hitter3,fc_score, by="id") # 年度、チーム名、選手名、因子得点のカラムを抽出 > fscores2 <- con_vm %>% dplyr::select(year,team,選手名,ML1,ML2,ML3) # 因子得点に名前を付ける > fscores %<>% rename(打撃技術=ML1, 長打力=ML2, 走力=ML3)
打撃技術
# 打撃技術 > bc <- fscores %>% dplyr::arrange(desc(打撃技術)) %>% head(5) > print(bc) year team 選手名 打撃技術 長打力 走力 1 2015 ヤクルト 川端 慎吾 2.939996 -1.06102638 0.3141611 2 2014 阪神 マートン 2.621416 -0.26765939 -0.6557168 3 2013 阪神 マートン 2.570096 0.03829577 -0.9561742 4 2013 ソフトバンク 内川 聖一 2.567277 0.29047422 -0.9587933 5 2013 ロッテ 今江 敏晃 2.472706 -0.81613314 -0.9445243
長打力
# 長打力 > pw <- fscores %>% dplyr::arrange(desc(長打力)) %>% head(5) > print(pw) year team 選手名 打撃技術 長打力 走力 1 2013 ヤクルト バレンティン -0.03436098 4.908928 -1.0891188 2 2011 西武 中村 剛也 -0.06175341 4.178670 -0.6219383 3 2016 ヤクルト 山田 哲人 -0.32908468 3.599799 0.9251124 4 2015 西武 中村 剛也 0.47319423 3.386818 -0.9606652 5 2016 横浜 筒香 嘉智 0.44410461 3.351923 -0.4769025
走力
# 走力 > run <- fscores %>% dplyr::arrange(desc(走力)) %>% head(5) > print(run) year team 選手名 打撃技術 長打力 走力 1 2011 ソフトバンク 本多 雄一 -0.1886137 0.09714669 4.458103 2 2016 西武 金子 侑司 -1.0446395 0.27608928 3.666863 3 2016 広島 田中 広輔 -0.3700709 1.12992915 3.596656 4 2012 中日 大島 洋平 0.9233752 -1.08544104 3.559731 5 2014 日本ハム 西川 遥輝 -0.5454467 1.17349420 3.545054
