dplyrとpurrrの読み方
日本人同士の会話だと
- dplyr : ディープライヤー
- purrr : プルル
と呼ぶことが多いと思います。
日本人以外と話すときには何と読めばよいか気になったため調べました。
youtubeで動画を探すと以下のように発音していました。
dplyr : ディープライヤー
https://www.youtube.com/watch?v=p-QgJtiwvjg
purrr : プァー
https://www.youtube.com/watch?v=aP_Zt6P9PXE
purrrの発音が難しいです。特に変わった読み方はせず単純に英語の発音ルールでpurr(猫の鳴き声を表す動詞。rは2つ)を読んでいるように思います。
Rでロジスティック回帰（２値分類）をやってみた
Ｒでロジスティック回帰（2値分類）の練習をやってみました。
概要
- NBAのデータを使用する
- 各年の成績データを使用して、その年にオールスターに選ばれたかを予測する（2値分類）
準備
最初に、ライブラリの読み込みと、モデルの説明変数（targets）を設定しておきます。
library(caret) library(formattable) library(nbastats) library(tidyverse) ### settings #target columns targets <- c( "MP", # minutes played "PTS", # points "FT", # free throws "AST", # assists "STL", # steals "BLK", # blocks "TOV" # tunrovers )
モデル構築用データ作成
次にデータを読み込んでモデル構築用のデータを作っていきます。今回はモデル構築に2015-2016年シーズンのデータを使用し、テストに2016-2017年シーズンのデータを使用します。
NBAの登録選手約480名に対し、オールスターに選出されるのは25名ほどのため、サンプルを各シーズンのサンプルを分割するのではなく、前年の年のデータでモデルを構築し、翌年のデータで検証する方針としました。
まずは、2015-2016年シーズンのデータを成績データ(seasons_stats)、オールスターデータ(allstars)からそれぞれ抽出します。
### data tranformation # load data stats <- seasons_stats stars <- allstars # train data (2016) stats <- stats %>% filter(Year == "2016") stars <- stars %>% filter(year == "2016")
次に、seasons_statsデータで１選手１レコードにします。これは、以下の選手のようにシーズン中に移籍した場合、各チームでの成績と、トータルの成績両方がデータとして存在しているためです。このような場合はトータルのデータを使用することにしました。
# in case of duplicate rows in stats (traded during season) # use total stats(TOT) stats <- stats %>% group_by(Player) %>% mutate(n = row_number(), max_n = max(n)) %>% filter(max_n != 1 & Tm == "TOT" | max_n == 1) %>% ungroup()
次にseasons_statsとallstarsをマージします。マージキーは選手の名前を使用し、全て小文字にしてから使用します。
# merge stats with all stars stats <- stats %>% mutate(name_s = tolower(Player)) stars <- stars %>% mutate(mflg = 1, name_a = tolower(name)) stats_stars <- left_join(stats, stars, by = c("name_s" = "name_a")) %>% mutate(is_allstar = coalesce(mflg, 0)) # check merge nrow(stars) stats_stars %>% group_by(is_allstar) %>% summarise(n = n()) > # check merge > nrow(stars) [1] 26 > stats_stars %>% + group_by(is_allstar) %>% + summarise(n = n()) # A tibble: 2 x 2 is_allstar n <dbl> <int> 1 0 450 2 1 26
マージの結果、26選手のオールスター選出情報が全て成績データと紐づきました。
モデル構築
次にいよいよモデルを構築します。概要で説明したとおり、2015-2016年シーズンのデータを全てモデル構築に使用します。
### modeling x <- stats_stars %>% select(targets) %>% sapply(as.integer) y <- stats_stars$is_allstar # standarise x <- scale(x) # make model model <- lm(y ~ MP + PTS + FT + AST + STL + BLK + TOV, data = x)
モデルが出来たので中身を確認します。
summary(model) > summary(model) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -0.60042 -0.05272 0.00042 0.03274 0.76162 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.054622 0.007321 7.461 4.21e-13 *** xMP -0.245296 0.023193 -10.576 < 2e-16 *** xPTS 0.250727 0.030534 8.211 2.15e-15 *** xFT 0.057598 0.019642 2.932 0.003529 ** xAST 0.090719 0.018559 4.888 1.40e-06 *** xSTL 0.061315 0.015763 3.890 0.000115 *** xBLK 0.055393 0.009783 5.662 2.61e-08 *** xTOV -0.104087 0.026695 -3.899 0.000111 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1597 on 468 degrees of freedom Multiple R-squared: 0.5143, Adjusted R-squared: 0.507 F-statistic: 70.79 on 7 and 468 DF, p-value: < 2.2e-16
各変数の係数を見ると、全て1%有意(**)を満たしております。係数の符号を見てみると以下のようになっています。
| 変数 | 符号 | 解釈 |
| xMP | - | プレイ時間が短い方がオールスターに選ばれやすい |
| xPTS | + | ポイント数が多い方がオールスターに選ばれやすい |
| xFT | + | フリースロー数が多い方がオールスターに選ばれやすい |
| xAST | + | アシスト数が多い方がオールスターに選ばれやすい |
| xSTL | + | スティール数が多い方がオールスターに選ばれやすい |
| xBLK | + | ブロック数が多い方がオールスターに選ばれやすい |
| xTOV | - | ターンオーバー数が少ない方がオールスターに選ばれやすい |
符号についても納得性が高いです。プレイ時間については一瞬逆のように思えますが、「他の条件（ポイント数、フリースロー数、アシスト数etc.）をコントロールした上では、出場時間が短い方が評価が高い」、つまり、「短い出場時間で好成績を出せた選手の評価が高い」と解釈すれば納得できると思います。
次に、２値分類のaccuracyを確認します。
### predict & check accuracy p_hat <- predict(model, data.frame(x)) y_hat <- ifelse(p_hat > 0.5, 1, 0) confusionMatrix(as.factor(y_hat), as.factor(y)) > confusionMatrix(as.factor(y_hat), as.factor(y)) Confusion Matrix and Statistics Reference Prediction 0 1 0 448 10 1 2 16 Accuracy : 0.9748 95% CI : (0.9564, 0.9869) No Information Rate : 0.9454 P-Value [Acc > NIR] : 0.001436 Kappa : 0.7145 Mcnemar's Test P-Value : 0.043308 Sensitivity : 0.9956 Specificity : 0.6154 Pos Pred Value : 0.9782 Neg Pred Value : 0.8889 Prevalence : 0.9454 Detection Rate : 0.9412 Detection Prevalence : 0.9622 Balanced Accuracy : 0.8055 'Positive' Class : 0
Balanced Accuracyは0.8055です。
次に、2016-2017年シーズンのデータでテストを行います。これまでと同じ手順を2016-2017年シーズンのデータに適用して、テスト用データを作成し、作成したモデル（model）で予測を行い、Balanced Accuracyを確認します。
### test in 2017 data data(seasons_stats) data(allstars) # use 2017 data stats_test <- seasons_stats %>% filter(Year == "2017") stars_test <- allstars %>% filter(year == "2017") # in case of duplicate rows in stats (traded during season) # use total stats(TOT) stats_test <- stats_test %>% group_by(Player) %>% mutate(n = row_number(), max_n = max(n)) %>% filter(max_n != 1 & Tm == "TOT" | max_n == 1) %>% ungroup() # merge stats with all stars stats_test <- stats_test %>% mutate(name_s = tolower(Player)) stars_test <- stars_test %>% mutate(mflg = 1, name_a = tolower(name)) stats_stars_test <- left_join(stats_test, stars_test, by = c("name_s" = "name_a")) %>% mutate(is_allstar = coalesce(mflg, 0)) # check merge stats_stars_test %>% group_by(is_allstar) %>% summarise(n = n()) ### predict x_test <- stats_stars_test %>% select(targets) %>% sapply(as.integer) y_test <- stats_stars_test$is_allstar # standarise x_test <- data.frame(scale(x_test)) p_test_hat <- predict(model, newdata = x_test) y_test_hat <- ifelse(p_test_hat > 0.5, 1, 0) confusionMatrix(as.factor(y_test_hat), as.factor(y_test)) > confusionMatrix(as.factor(y_test_hat), as.factor(y_test)) Confusion Matrix and Statistics Reference Prediction 0 1 0 459 10 1 2 15 Accuracy : 0.9753 95% CI : (0.9573, 0.9872) No Information Rate : 0.9486 P-Value [Acc > NIR] : 0.002596 Kappa : 0.7019 Mcnemar's Test P-Value : 0.043308 Sensitivity : 0.9957 Specificity : 0.6000 Pos Pred Value : 0.9787 Neg Pred Value : 0.8824 Prevalence : 0.9486 Detection Rate : 0.9444 Detection Prevalence : 0.9650 Balanced Accuracy : 0.7978 'Positive' Class : 0
Balanced Accuracyは0.7978で訓練データの時とほとんど同じです。そのためモデルはある程度頑健といえると思われます。
備考
今回の目的はロジスティック回帰の手順確認のため、訓練データとテストデータについては吟味しておりません。より頑健性の高いモデルを作るには、例えば複数シーズンのデータをモデル構築に使用したり、シーズン内のデータからサンプリングを行ったりする必要があると思います。
Rでスクレイピングをしてみた part2
Ｒでロジスティック回帰の練習をしようと思ったのですが、自分がもっているデータの中に２値分類の問題がつくれそうなデータがなかったため、新たにデータを取得することにしました。
今回はwikipediaからNBAのオールスター出場選手のデータを取得します。
List of NBA All-Stars - Wikipedia
やり方は以前スクレイピングをやったときと同じです。今回はネットから情報を取得した後のデータ加工が少し大変でした。
データの取得
まずはrvestパッケージを使用してウェブからhtmlファイルを取得します。
library(rvest) library(tidyverse) ### get data from web # set url url <- "https://en.wikipedia.org/wiki/List_of_NBA_All-Stars" # get html file h <- read_html(url)
次にhtmlファイルの中からデータの入っているテーブルを探します。
# parse html table <- h %>% html_nodes("table") # check the table tag data table > # check the table tag data > table {xml_nodeset (14)} [1] <table class="wikitable"><tbody>\n<tr>\n<td>\n<b>#</b>\n</ ... [2] <table class="wikitable sortable" style="width:100%"><tbod ... [3] <table class="nowraplinks hlist collapsible autocollapse n ... [4] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [5] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [6] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [7] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [8] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [9] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [10] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [11] <table style="width:100%;padding:0"><tbody><tr>\n<td style ... [12] <table class="nowraplinks hlist collapsible autocollapse n ... [13] <table class="nowraplinks navbox-subgroup" style="border-s ... [14] <table class="nowraplinks navbox-subgroup" style="border-s ...
テーブルがたくさんあります。
wikipediaのソースを表示して確認すると(table class="wikitable sortable" style="width:100%")のタグがついたテーブルが今回欲しいテーブルのようです。
そのため2番目のテーブルデータを取得します。
データを取得したらView()で中身を確認します。
# pick up all-star players table allstars <- table[[2]] %>% html_table() # check the table View(allstars)
データはうまくデータフレームに入っています。
データの整形
後はデータをキレイにしていきます。
今回は以下のように、名前とAll Starに選出された年が入った2列のデータを目指します。
まずは、列名を変更します。
### transform data # change names names(allstars) <- c("Player", "counts", "Selections", "Notes", "Reference") View(allstars)
次に実際に使用するPlayer列とSelections列のデータをキレイにします。
Player列は不要な記号が末尾についているため、正規表現を使用して、名前の部分だけを取り除きます。
併せて、Selections列にenダッシュがはいっているため、これをハイフンに変換します。
enダッシュの変換方法については以下を参考にしました。
character - How can I add an en dash to a plot in R? - Stack Overflow
library(Cairo) # clean data allstars <- allstars %>% mutate(names = str_extract(Player, "[A-Za-z'. -]+"), years = gsub("\U2013", "-", Selections)) View(allstars)
これでPlayer列とSelection列をきれいにしたnames列とyears列が作成できました。
確認してみると一部Player列→names列の作成がうまくできていないデータがあります。
理由は選手の名前がEnglish alphabet以外のアルファベットで表記されているためです。
選手の指名の表記に応じて複数の変換ロジックを組むのが大変だったため、今回は表記がおかしくなる5名の選手名を直接修正することにしました。
# fix irregular(non-English alphabet) name players allstars[244, "names"] <- "Zydrunas Ilgauskas" allstars[324, "names"] <- "Goran Dragic" allstars[355, "names"] <- "Nikola Jokic" allstars[387, "names"] <- "Kristaps Porzingis" allstars[414, "names"] <- "Nikola Vucevic"
以降はnames列とyears列を使用して、目的のデータを作っていきます。
かなりムリヤリになってしまいました。
まずnamesとyearsをそれぞれベクトルとして抽出します。
これらを加工して得たベクトルを初期化したallstarsというデータフレームに格納していきます。
### create data table # separate names column and years column as single dataframe names <- allstars["names"] years <- allstars["years"] # initialise allstars allstars <- data.frame(name = NULL, year = NULL)
ここで、yearsの各行から年を抽出してベクトルを作成します。ポイントは連続して選出された場合2000-2005のようにハイフンでつながっていることです。そのため、ベクトル作成用の関数を作成します。
### utility selection_years <- function(str) { # split string str_array <- str_split(str, ";", simplify = TRUE) # split by ";" str_array <- gsub(" ", "", str_array) # remove space str_array <- gsub(",", "", str_array) # remove camma # number of elements of selection years ncol <- ncol(str_array) selection_years <- c() for(i in 1:ncol) { # single year(e.g 1997) or multiple years(e.g. 1997-2000) is_single <- !(str_detect(str_array[1, i], pattern = "-")) # swith the way of making vector by is_single if (is_single) { selection_years <- append(selection_years, as.integer(str_array[1, i])) } else { start <- as.integer(str_split(str_array[1, i], "-", simplify = TRUE)[1, 1]) end <- as.integer(str_split(str_array[1, i], "-", simplify = TRUE)[1, 2]) selection_years <- append(selection_years, start:end) } } selection_years }
ダラダラ長くなってしまいましたが、yearsに格納された要素を行列に変換して(1×N列)、各列が単一の年（2000）の場合と、複数年の場合(2000-2005)で処理を分岐させて1つのベクトルを作成しています。
この関数を使用して、yearsから選出年が1つずつ格納されたベクトルを作成して、選出された年数と同じ長さの名前ベクトルをnamesから作成してつなげます。
これを全ての選手について繰り返すことで目的のデータを得ます。
# pick up selection years from years # , and combine it with names n <- nrow(names) for(i in 1:n) { n <- names[i, 1] ys <- selection_years(years[i, 1]) l <- length(ys) ns <- rep(n, l) allstars <- rbind(allstars, data.frame(name = ns, year = ys)) }
これで完成です。
後は、以下の手順で自作のパッケージに今回作成したallstarsのデータを加えていつでも使えるようにしておきます。
hadadada00.hatenablog.com
Rで因子分析をやってみた
Ｒによるデータ分析の練習として因子分析をやってみました。
手順は以下の記事を参考にしました。
Rで因子分析やってみた - Qiita
概要
- NBAの各選手の成績データを使用する
- PT（ポイント数）やBLK（ブロック数）などの成績データから選手のパフォーマンスを説明する共通因子がないか確認する
データ
以下の記事で作成したデータを使用します。
hadadada00.hatenablog.com
データ準備
分析対象の変数はtargetsで指定した5つとします。
library(formattable) library(nbastats) library(tidyverse) ### data # stats data stats <- seasons_stats # filter 2017 records stats <- seasons_stats %>% filter(Year == 2017) # targets for analysis targets <- c( "AST", # assists "STL", # steals "BLK", # blocks "TRB", # total rebounds "PTS" # points ) # extract variables and change data type vars <- stats %>% select(targets) %>% sapply(as.integer)
scree plot
因子数を決めるためにスクリープロットを作成します。
### scree plot # correlation matrix cor <- cor(vars) # eigen value ev <- eigen(cor)$values # scree plot plot(ev, type = "b")
3つ目以降の固有値はほとんど説明力をもっていないため、因子数は2とします。
因子分析の実施
factanal関数を使用して因子分析を実施します。以下の記事を参考に、因子の抽出法は最尤法、回転法はプロマックス法を使用します。また、因子得点は回帰法で算出します。
因子分析の因子抽出方法について | Sunny side up!
因子分析における因子軸の回転法について | Sunny side up!
### factor analysis # factor extraction : maximum likelihood # rotation : promax # factor score : regression fa <- factanal(x = vars, factors = 2, rotation = "promax", scores = "regression") fa
結果は以下です。
>fa Call: factanal(x = vars, factors = 2, scores = "regression", rotation = "promax") Uniquenesses: AST STL BLK TRB PTS 0.137 0.225 0.300 0.047 0.192 Loadings: Factor1 Factor2 AST 1.036 -0.202 STL 0.753 0.190 BLK -0.151 0.919 TRB 0.923 PTS 0.683 0.304 Factor1 Factor2 SS loadings 2.138 1.866 Proportion Var 0.428 0.373 Cumulative Var 0.428 0.801 Factor Correlations: Factor1 Factor2 Factor1 1.000 -0.601 Factor2 -0.601 1.000 Test of the hypothesis that 2 factors are sufficient. The chi square statistic is 0.21 on 1 degree of freedom. The p-value is 0.643
Cumulative Varを見ると、2つの因子で全体の80%を説明できています。
Factor1とFacto2の内訳を見ていると、Factor1はAST（アシスト数）、STL（スティール数）、PTS（ポイント数）が高い因子負荷量を示しており、Factor2はBLK（ブロック数）とTRB（リバウンド数）が高い因子負荷量を示しています。
因子の名前としては以下のようにつけることにしました。
- 因子1 ： non-power-play ability
- 因子2 ： power-play ability
因子得点の確認
因子得点のトップ10を求めてどのような選手がそれぞれの因子得点が高いか確認します。まずはFactor1（non-power-play ability）の上位10選手を見てみます。
### check factor scores stats <- cbind(stats, fa$scores) # top10 non-power-play ability tbl1 <- stats %>% arrange(desc(Factor1)) %>% mutate(n = row_number()) %>% filter(n <= 10) %>% select(Player, Pos, Factor1) formattable(tbl1)
ランキングを見ると、ポイントガードのスター選手がランクインしています。Facto1はポイントガード能力といってもいいかもしれません。
次にFactor2（power-play ability）の上位10選手を見てみます。
# top10 power-play ability tbl2 <- stats %>% arrange(desc(Factor2)) %>% mutate(n = row_number()) %>% filter(n <= 10) %>% select(Player, Pos, Factor2) formattable(tbl2)
ランキングを見るとセンターのスター選手がランクインしています。Factor2はセンター能力といってもいいかもしれません。
（以上）
Rでスクレイピングしたデータをpackageにして公開してみた
以下の記事でウェブスクレイピングの練習をしました。
hadadada00.hatenablog.com
せっかくなので以下の記事で作成したpackageに追加して今後分析で簡単に使えるようにしてみました。
hadadada00.hatenablog.com
前回作成したpackageのローカル環境は以下のようになっています。
こちらにファイルを追加していきます。
まず、data-raw下にウェブスクレイピング用のRスクリプトファイルを作成します。
read_html.R
# TITLE : read_html.R # PURPOSE : get salary data from web # HISTORY : Ver1 2019-06-16 library(rvest) library(tidyverse) # url(domain) url_domain <- "https://hoopshype.com/salaries/players/" # get_html function # param season : NBA season get_html <- function(season) { url <- paste0(url_domain, season) # get html h <- read_html(url) # parse html tab <- h %>% html_nodes("table") # get salary table from tab table <- tab[[1]] %>% html_table() # create data frame table <- table %>% mutate(n = row_number(), season = season, name = X2, salary = X3, salary_adj = X4) %>% filter(n != 1) %>% select(season, name, salary, salary_adj) } # initialise data frame salaries <- data.frame() # repeat get_html function and add salary table into salaries salaries <- rbind(salaries, get_html("2017-2018"), get_html("2016-2017"), get_html("2015-2016"), get_html("2014-2015"), get_html("2013-2014"), get_html("2012-2013"), get_html("2011-2012"), get_html("2010-2011"), get_html("2009-2010"), get_html("2008-2009"), get_html("2007-2008"), get_html("2006-2007"), get_html("2005-2006"), get_html("2004-2005"), get_html("2003-2004"), get_html("2002-2003"), get_html("2001-2002"), get_html("2000-2001")) # save as .rda file save(salaries, file = "../data/salaries.rda")
これでdataディレクト下に新たにsalaries.rdaというデータファイルが追加されました。
次にRディレクト下にsalaries.Rという名前でデータの説明ファイルを作成します。
#' Datasets of NBA players' salaries #' #' A dataset of NBA players' salaries since 2000-2001 season to 2017-2018 season. #' #' @format A data frame with 8,691 rows and 4 variables: #' \describe{ #' \item{season}{season year} #' \item{name}{players' name} #' \item{salary}{salary, in US dollar} #' \item{salary_adj}{salary adjusted for inflation, in US dollar} #' } #' @docType data #' #' @usage data(salaries) #' #' @keywords datasets #' #' @source \url{https://hoopshype.com/salaries/players/2017-2018/} #' "salaries"
次にドキュメントファイル（.Rd）をアップデートします。
devtools::document()
R/man下にsalaries.Rdが追加されました。
次にパッケージのロードとチェックを行います。
チェックでエラーがでなければローカルでの編集作業は終了です。
devtools::load_all() devtools::check()
チェックで問題がなければgithubにアップロードします。
コマンドプロンプトを開いてpackageの環境に移動します。
cd /Users/[your name]/Desktop/Rwork/nbastats
あとはおなじみのコマンドです。
git add . git commit -m "[your comments]" git push -u origin master
最後に適当なRスクリプトファイルでパッケージのインストールからデータの読み込みまでを確認します。
devtools::install_github("hadadada00/nbastats") library(nbastats) head(salaries) > head(salaries) season name salary salary_adj 1 2017-2018 Stephen Curry $34,682,550 $35,678,476 2 2017-2018 LeBron James $33,285,709 $34,241,524 3 2017-2018 Paul Millsap $30,769,231 $31,652,784 4 2017-2018 Gordon Hayward $29,727,900 $30,581,550 5 2017-2018 Blake Griffin $29,512,900 $30,360,377 6 2017-2018 Kyle Lowry $28,903,704 $29,733,687
しっかり読み込めています。これでnbastatsパッケージににsalariesのデータを追加できました。
Rでスクレイピングをやってみた
以下の記事で分析を行う際に、ウェブページからNBA選手の年俸データを取得しました。
年毎にページが分かれているため、Rでウェブスクレイピングをして、データ取得を自動化しようと思います。
ウェブスクレイピングのやり方については以下のサイトを参考にしました。
Chapter 24 Web Scraping | Introduction to Data Science
ライブラリはrvestを使用します。まずは、データ取得元のURLをセットします。
library(rvest) library(tidyverse) # set url url_domain <- "https://hoopshype.com/salaries/players/" url_page <- "2017-2018" url <- paste0(url_domain,url_page)
次に、htmlファイルを取得し、tableタグの中身を抽出します。
# get html h <- read_html(url) # parse html tab <- h %>% html_nodes("table")
tableタグがついている要素を確認します。
tab > tab {xml_nodeset (1)} [1] <table class="hh-salaries-ranking-table hh-salaries-table- ...
tableタグがついている要素は1つだけのため、これを抽出してデータフレームにすればよさそうです。
データフレームへの変換後、データの中身を確認します。
# get salary table from tab table <- tab[[1]] %>% html_table() # check the table head(table) >head(table) X1 X2 X3 X4 1 NA Player 2017/18 2017/18(*) 2 1 Stephen Curry $34,682,550 $35,678,476 3 2 LeBron James $33,285,709 $34,241,524 4 3 Paul Millsap $30,769,231 $31,652,784 5 4 Gordon Hayward $29,727,900 $30,581,550 6 5 Blake Griffin $29,512,900 $30,360,377
データは取得できていますが、1行目がカラム名になっており、カラム名が未定義になっています。
また、X1列は不要です。
不要な列の削除、列名の設定、必要な列（年）の追加をしてデータフレームを整えます。
# create data frame table <- table %>% mutate(n = row_number(), year = "2017-2018", name = X2, salary = X3, salary_adj = X4) %>% filter(n != 1) %>% select(year, name, salary, salary_adj) head(table) >head(table) > head(table) year name salary salary_adj 1 2017-2018 Stephen Curry $34,682,550 $35,678,476 2 2017-2018 LeBron James $33,285,709 $34,241,524 3 2017-2018 Paul Millsap $30,769,231 $31,652,784 4 2017-2018 Gordon Hayward $29,727,900 $30,581,550 5 2017-2018 Blake Griffin $29,512,900 $30,360,377 6 2017-2018 Kyle Lowry $28,903,704 $29,733,687
これでデータスクレイピングができました。
追記：
取得元のサイトでは2000-2001シーズンから、2017-2018シーズンまでのデータが取得できます。
今後の分析で使用することを踏まえて、まとめて取得してpackage化しておこうと思います。
作業内容は別記事（以下）に記載しております。
hadadada00.hatenablog.com
以上
Rで主成分分析をやってみた
統計の勉強で主成分分析をやってみました。
分析概要
- NBA選手のstats（各成績のデータ）を使用して、選手の総合力を測る指標を作成する
- データは2016-2017年シーズンのものを使用
- 結果の大まかな検証として年収のデータと比較する
データ
選手の成績データは以下の記事で作成したものを使用します。
hadadada00.hatenablog.com
年収のデータは以下のリンクから入手したものをcsvファイルにして読み込んで使用します。
https://hoopshype.com/salaries/players/2016-2017/
データの読み込み
まずは成績のデータと年収のデータを読み込みます。
成績のデータは項目数が多いのですが、今回はtargetsで指定した4変数のみ使用します。
library(formattable) library(ggbiplot) library(nbastats) library(tidyverse) ### data # stats data stats <- seasons_stats # filter 2017 records stats <- seasons_stats %>% filter(Year == 2017) # targets for analysis targets <- c( "AST", # assists "STL", # steals "BLK", # blocks "PTS" # points ) # change data type temp1 <- stats %>% select(-targets) temp2 <- stats %>% select(targets) %>% sapply(as.integer) stats <- cbind(temp1, temp2) # salary data # https://hoopshype.com/salaries/players/2016-2017/ salaries <- read_csv("./salary_2016_2017.csv")
成績データと年収データのマージ
次に読み込んだ成績データと年収データをマージします。
stats <- stats %>% mutate(flg_st = "1", name = tolower(Player)) salaries <- salaries %>% mutate(flg_sa = "1", name = tolower(name)) stats_salaries <- left_join(stats, salaries, by = c("name"))
マージ結果を確認します。
stats_salaries %>% group_by(flg_st, flg_sa) %>% dplyr::summarise(n = n()) # A tibble: 2 x 3 # Groups: flg_st [1] flg_st flg_sa n <chr> <chr> <int> 1 1 NA 44 2 1 1 551
一部くっついていないデータがありますが、選手の母集団が違う可能性があるため、今回はこれでよしとしようと思います。
年収が取得できなかったレコードは除外します。
stats_salaries <- stats_salaries %>% filter(!is.na(flg_sa))
少しデータを見てみようと思います。成績の相関係数と年収のTOP10です。
### check the data #correlation matrix cor <- data.frame(cor(stats_salaries[targets])) formattable(cor) # salary top5 bottom5 top5_bottom5 <- stats_salaries %>% arrange(desc(salary)) %>% dplyr::mutate(n = row_number(), salary = comma(salary)) %>% filter(n <= 5 | n >= (nrow(stats_salaries) - 4)) %>% select(n, Player, Pos, Age, Tm, salary) formattable(top5_bottom5)
相関係数行列を見ると4つの変数間に相関は概ねありそうですが、BLK（ブロック数）だけ、やや相関が低いといえそうです。
年収（USドル）のランキングを見ると、1位のLebron Jamesは3,000万ドル（約30億円）もらっているのに対して、一番低いDahntay Jonesは24,000ドル（約240万円）です。
主成分分析の実施
以下で主成分分析を実行します。
### principal component analysis # pick up target variables pca_vars <- stats_salaries[, targets] # pca pca_results <- prcomp(pca_vars, scale = T)
結果を確認します。
summary(pca_results) pca_results > summary(pca_results) Importance of components: PC1 PC2 PC3 PC4 Standard deviation 1.6865 0.8834 0.44539 0.4205 Proportion of Variance 0.7111 0.1951 0.04959 0.0442 Cumulative Proportion 0.7111 0.9062 0.95580 1.0000 > pca_results Standard deviations (1, .., p=4): [1] 1.6865263 0.8834317 0.4453887 0.4204833 Rotation (n x k) = (4 x 4): PC1 PC2 PC3 PC4 AST 0.5096261 -0.46438991 0.09538605 -0.7180005 STL 0.5480113 -0.13136200 -0.73544490 0.3762293 BLK 0.3675924 0.87550923 -0.04149933 -0.3108652 PTS 0.5521238 -0.02386692 0.66955217 0.4962757
累積寄与率（Cumulatrive Proportion）を見ると、第1主成分で71%ほど、第2主成分までで全体の91%の変動が説明できています。
結果をグラフでも確認してみます。
ggbiplot(pca_results, obs.scale = 1, var.scale = 1, alpha=0.5)
グラフを見ると第1主成分は概ねPTS（総得点数）、STL（スティール数）、AST（アシスト数）と相関が高く、第2主成分はBLK（ブロック数）と相関が高いです。
第1主成分はいずれの成績（PTS, STL, AST, BLK）とも正の相関のため、選手の能力の総合指標として使えそうです。
主成分スコア（PC1）と年収の比較
作成した主成分スコアが選手の能力を表わす総合指標として使えそうか、主観的な評価と年収との比較で確認します。
まずは、PC1のトップ10選手を確認します。
### compare PC1 scores with Salaries stats_salaries_scores <- cbind(stats_salaries, pca_results$x) # make rank variables stats_salaries_scores <- stats_salaries_scores %>% mutate(salary_rank = rank(desc(salary)), score1_rank = rank(desc(PC1))) # top10 score1 players stats_salaries_scores %>% arrange(score1_rank) %>% filter(score1_rank <= 10) %>% select(Player, Pos, Age, Tm, score1_rank, salary_rank) %>% formattable()
ランクインしている選手を見ると、salary_rankでも上位のスター選手がいる一方、salary_rankがそれほど高くない選手もいます。
以下に挙げたように年報と選手の評価やパフォーマンスは必ずしも一致しないため注意が必要ですが、PC1は選手のパフォーマンスを表わす指標として的外れではなさそうです。
- 2016-2017年にブレイクした選手は、その評価が年俸に反映されていない。
- 若い選手は契約の制約でそもそも年俸が上がりづらい
主観的な感想としては、ここにあがっている選手は現役のスター選手ばかりのためPC1がパフォーマンス指標としてある程度使えそうだと思われます。
次にPC1とsalaryの相関を見ます。
# correlation matrix stats_salaries_scores %>% select(score1_rank, salary_rank) %>% cor() > score1_rank 1.0000000 0.5784745 salary_rank 0.5784745 1.0000000 # scatter plot of score1_rank and salary_rank stats_salaries_scores %>% ggplot(aes(score1_rank, salary_rank)) + geom_point()
相関係数は0.58程度のため、そこまで高くはありません。散布図をみてもかなりバラけています。
総合的に見ると、主観的にはPC1は選手の総合的な能力を表わす指標として使えそうです。しかし、年俸といった選手の現時点での金銭的な価値とはかならずしも相関が高くないといえそうです。
補足
今回は主成分分析を用いて、主成分スコアを選手の総合能力を表わす指標として使えないか考えて見ました。しかし、わざわざ主成分スコアを計算しなくとも、各成績の平均値ではだめなのか気になりました。
ためしに、標準化してスケールを調整した上での成績の平均値を計算し、主成分スコアと比較してみました。
### supplement stats_salaries_scores %>% mutate(AST_std = scale(AST), STL_std = scale(STL), BLK_std = scale(BLK), PTS_std = scale(PTS), std_sum = (AST_std + STL_std + BLK_std + PTS_std), std_sum_rank = rank(desc(std_sum))) %>% ggplot(aes(score1_rank, std_sum_rank)) + geom_point()
結果を見ると、主成分1スコアと各成績の平均値（標準化後）はほぼおなじ指標であることがわかりました。そのため、今回の例では主成分をわざわざしなくとも成績指標を平均するだけで簡易的なパフォーマンス指標として使えると思います。
