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])

f:id:joure:20170610003031p:plain

因子分析

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))

f:id:joure:20170611010017p:plain

「長打力」「打撃技術」「走力」別の因子得点が上位の選手

 さて因子得点を算出したとしても、選手名や成績の年度、球団名といったカラムがありません。なので、色々方法があると思うのですが、ここではマトリックス形式の因子得点をデータフレーム形式にして元のデータフレームと結合(結合キーも作成する)します。そして結合後のデータフレームから年度、チーム名、選手名、因子得点のカラムを抽出して、長打力」「打撃技術」「走力」それぞれの因子得点が上位の選手を表示したいと思います。

# データフレーム形式にする
> 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