聰明不如鈍筆
총명불여둔필
assignment KindeR

R로 깔끔하게 머신러닝(랜덤 포레스트) 요리하기(feat. tidymodels)


언제인가부터 머신러닝 또는 기계학습이라는 표현이 유행하기 시작했습니다.


이 블로그와 참 좋은 친구인 위키피디아는 머신러닝을 이렇게 설명합니다.


기계학습(機械學習) 또는 머신러닝(영어: machine learning)은 인공지능의 한 분야로, 컴퓨터가 학습할 수 있도록 하는 알고리즘과 기술을 개발하는 분야를 말한다. 가령, 기계학습을 통해서 수신한 e메일스팸인지 아닌지를 구분할 수 있도록 훈련할 수 있다.


기계 학습의 핵심은 표현(representation)과 일반화(generalization)에 있다. 표현이란 데이터의 평가이며, 일반화란 아직 알 수 없는 데이터에 대한 처리이다. 이는 전산 학습 이론 분야이기도 하다. 다양한 기계 학습의 응용이 존재한다. 문자 인식은 이를 이용한 가장 잘 알려진 사례이다.


물론 이 331자만 가지고 머신러닝이 무엇인지 정확하게 알 수는 없지만 알 듯 모를 듯 감은 잡을 수 있습니다.


'R로 깔끔하게 통계 모형 쓸어담기(feat. broom)' 포스트 마지막에 소개해 드린 K 평균 군집 분석이 바로 머신러닝 기법 가운데 하나입니다.


이렇게 R로는 (적어도 간단한) 머신러닝 작업을 어렵지 않게 진행할 수 있습니다.


아직 머신러닝 작업을 한번도 해보지 않으신 분은 막상 해보면 너무 간단해서 놀라실지도 모릅니다.


이번 포스트에서는 지난해 9월 29일 Sportugese에 '머신러닝으로 2019 사이영상 수상자 예상하기' 포스트를 쓰면서 했던 작업을 되살려 보겠습니다.


노파심에 말씀드리면 같은 작업을 새로 하는 거라 그때하고 100% 똑같은 결과가 나온다고 장담할 수는 없습니다.



랜덤 포레스트란 무엇인가?

'머신러닝으로 2019 사이영상 수상자 예상하기'에서 사용한 머신러닝 기법은 '랜덤 포레스트(Random Forest)'였습니다.


랜덤 포레스트를 이해하시려면 먼저 '의사결정 나무'가 무엇인지 알아봐야 합니다.


의사결정 나무는 "의사 결정 규칙과 그 결과들을 트리 구조로 도식화한 의사 결정 지원 도구의 일종"(위키피디아 한국어판)입니다. 


말로는 어려워도 아래 그림을 보시면 의사결정 나무가 뭔지 짐작하실 수 있을 겁니다.



이 나무는 Sportugese에서 프로배구 2018~2019 V리그 전반기 시청률을 다루면서 그렸던 것. 


이런 나무 하나로는 의사결정 과정에 '구멍'이 있을 수 있습니다.


이때 이런 나무를 여러 개 그리면 구멍을 최대한 줄일 수 있습니다. 나무를 여러 개 그려서 다수결 또는 평균으로 최종 결과물을 내놓는 방식입니다. 



이런 방식을 의사결정 나무가 랜덤하게 숲을 이룬다는 뜻에서 랜덤 포레스트라고 부릅니다.


개념을 이해하셨을 줄로 믿고 작업을 진행해 보겠습니다.



패키지, 자료 불러오기

언제나 그렇듯 제일 먼저 할 일은 패키지와 자료 불러오기이고, 언제나 그렇듯 제일 먼저 불러 올 패지키는 tidyverse입니다.

library('tidyverse')
## -- Attaching packages -------------------------------------------------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.2.1     √ purrr   0.3.3
## √ tibble  2.1.3     √ dplyr   0.8.3
## √ tidyr   1.0.0     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.4.0
## -- Conflicts ----------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()


우리는 tidy한 방법으로 머신러닝 모델을 만들 계획이니까 tidymodels도 불러옵니다.

library('tidymodels')
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## -- Attaching packages ------------------------------------------------------------------------------- tidymodels 0.0.3 --
## √ broom     0.5.3     √ recipes   0.1.9
## √ dials     0.0.4     √ rsample   0.0.5
## √ infer     0.5.1     √ yardstick 0.0.4
## √ parsnip   0.0.5
## -- Conflicts ---------------------------------------------------------------------------------- tidymodels_conflicts() --
## x scales::discard()   masks purrr::discard()
## x dplyr::filter()     masks stats::filter()
## x recipes::fixed()    masks stringr::fixed()
## x dplyr::lag()        masks stats::lag()
## x dials::margin()     masks ggplot2::margin()
## x yardstick::spec()   masks readr::spec()
## x recipes::step()     masks stats::step()
## x recipes::yj_trans() masks scales::yj_trans()


위에서 보신 것처럼 두 패키지 모두 한 번에 여러 패키지를 불러옵니다.


필요한 패키지만 불러오지 않고 이렇게 동시에 불러오는 건 실제로 작업을 하다 보면 각 패키지를 유기적으로 연결해 쓰는 게 능률을 높이는 데 도움이 되기 때문입니다.


자료도 불러와야겠죠? 우리가 이번 작업에 쓸 파일은 아래 CSV 파일에 들어 있습니다.


cya.csv


CSV 형식은 그냥 쉼표로 구분한 값(Comma-Seperated Values)이라는 뜻으로 각 열을 쉼표로 구분해 이런 이름이 붙었습니다.


기본적으로 R에는 CSV 파일을 읽을 때 쓰라고 read.csv() 함수가 들어 있습니다.


물론 이 함수를 써도 되지만 우리는 tidyverse(정확하게는 readr)에 들어 있는 read_csv() 함수를 쓰겠습니다.

cya <- read_csv('cya.csv')
## Parsed with column specification:
## cols(
##   season = col_double(),
##   name = col_character(),
##   team = col_character(),
##   lg = col_character(),
##   cy = col_character(),
##   vote = col_double(),
##   w = col_double(),
##   l = col_double(),
##   war = col_double(),
##   era = col_double(),
##   fip = col_double(),
##   ip = col_double(),
##   k = col_double(),
##   bb = col_double(),
##   hr = col_double(),
##   k9 = col_double(),
##   bb9 = col_double(),
##   kbb = col_double(),
##   babip = col_double()
## )


잘 들어왔는지 볼까요?

cya
## # A tibble: 807 x 19
##    season name  team  lg    cy     vote     w     l   war   era   fip    ip
##     <int> <fct> <fct> <fct> <fct> <int> <int> <int> <int> <int> <int> <int>
##  1   2019 Gerr~ Astr~ AL    X         0     2     2     1     1     1     3
##  2   2019 Lanc~ Rang~ AL    X         0     4    15     2     7     3     4
##  3   2019 Just~ Astr~ AL    X         0     1     3     3     2     4     1
##  4   2019 Char~ Rays  AL    X         0     4     3     4     3     2     9
##  5   2019 Shan~ Indi~ AL    X         0     7     8     5     4     5     2
##  6   2019 Luca~ Whit~ AL    X         0     9    10     6     5     6    15
##  7   2019 Jose~ Twins AL    X         0     9     8     7     8     7     8
##  8   2019 Mike~ Rang~ AL    X         0     9    14     8     6    11     4
##  9   2019 Edua~ Red ~ AL    X         0     3     3     9     9     8     6
## 10   2019 Marc~ Mari~ AL    X         0     4    19     9    13    10     7
## # ... with 797 more rows, and 7 more variables: k <int>, bb <int>, hr <int>,
## #   k9 <int>, bb9 <int>, kbb <int>, babip <int>


이 파일에는 2009년부터 지난해까지 11년 동안 규정이닝을 채운 투수들을 대상으로 △시즌(season) △이름(name) △팀(team) △리그(lg) △사이영상 수상 여부(cy) △사이영상 투표 득점(vote) △승리 순위(w) △패배 순위(l) △팬그래프스 대체 선수 대비 승리 기여도(fWAR) 순위(war) △평균자책점 순위(era) △수비 영향을 제거한 평균자책점(FIP) 순위(fip) △투구 이닝 순위(ip) △삼진 순위(k) △볼넷 순위(bb) △피홈런 순위(hr) △9이닝당 탈삼진(K/9) 순위(k9) △9이닝당 볼넷(BB/9) 순위(bb9) △삼진 대 볼넷 비율(K/BB) 순위(kbb) △인플레이타율(BABIP) 순위(babip)가 들어 있습니다.


이 가운데 2019년 데이터는 아직 사이영상 수상 여부도 없고, 투표 결과도 없기 때문에 일단 따로 빼겠습니다. 

cya %>% 
  filter(season==2019) -> cya_2019


'%>%'는 파이프(pipe) 기호입니다. 현실 세계에서 파이프가 액체나 기체 같은 유체를 한 곳에서 다른 곳으로 보내는 것처럼 코딩에서 파이프는 한 자료를 이 함수에서 저 함수로 보내는 구실을 합니다.


filter()dplyr 패키지에 들어 있는 함수(동사)로 문자 그대로 특정 기준에 따라 행을 골라내라는 뜻입니다.


이 코드는 season이 2019와 같은 열만 골라내라는 뜻이겠죠?


R에서 '='는 같다는 뜻이 아니라 변수에 데이터를 넣으라는 뜻이고, '=='가 같다는 의미입니다.


끝에 '->'라고 쓴 건 이렇게 정리한 결과를 특정 변수(여기서는 cya_2019)에 넣어두라는 뜻입니다.


그래서 저 코드는 사실 아래와 같은 뜻입니다.

cya_2019 <- filter(cya, season==2019)


이런 작업을 어떻게 진행하는지 잘 모르겠다고 하시는 분은 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.


다시 코딩으로 돌아와서 2009~2018년 데이터도 cya_pre라는 별도 변수에 넣어놓겠습니다.

cya %>% 
  filter(season!=2019) -> cya_pre


이것으로 데이터 불러오기도 일단 끝이 났습니다. 이제 본격적으로 머신러닝 작업을 시작할 차례입니다.



요리 재료 준비

머신러닝 첫 단계는 데이터를 학습용(또는 연습용)과 시험용(또는 검증용) 둘로 나누는 것.


이렇게 데이터를 둘로 나누는 건 만약 학습용 데이터로 시험을 보면 시험 문제를 다 알고 시험을 보는 것과 마찬가지 결과가 되기 때문입니다.


그래서 원래 데이터 특성을 (거의) 똑같이 유지하는 두 그룹으로 나누는 게 중요합니다.


tidymodels에서는 initial_split() 함수가 이 기닝을 담당합니다.    


아래는 cya_pre 데이터 가운데 70%(=0.7)를 학습용으로 (자연스레 나머지 30%를 시험용) 골라내 이를 cya_split라는 변수에 넣으라는 뜻입니다.

cya_pre %>%
  initial_split(prop=0.7) -> cya_split


어떻게 나눴는지 한 번 확인해 보겠습니다.

cya_split
## <526/225/751>


그냥 숫자가 나옵니다. 이건 학습용 데이터가 526개, 시험용 데이터가 225개고, 전체 데이터는 751개라는 뜻입니다.


어떤 데이터가 학습용인지 알아볼 때는 training() 함수를 쓰면 됩니다.

cya_split %>%
  training()
## # A tibble: 526 x 19
##    season name  team  lg    cy     vote     w     l   war   era   fip    ip
##     <int> <fct> <fct> <fct> <fct> <int> <int> <int> <int> <int> <int> <int>
##  1   2018 Marc~ Mari~ AL    X         0    12    13     9    17     8    19
##  2   2018 Dyla~ Orio~ AL    X         0    23    25    25    25    25     2
##  3   2018 Luis~ Yank~ AL    X         1     3     9     5     9     5    21
##  4   2018 Mike~ Mari~ AL    X         0    18    15    20    20    18     9
##  5   2018 Dall~ Astr~ AL    X         0    14    19    10    14    11    16
##  6   2018 Trev~ Indi~ AL    X        13    13     5     3     2     1    26
##  7   2018 Andr~ Ange~ AL    X         0    20    16    14    18    14    13
##  8   2018 Gerr~ Astr~ AL    X        26     9     3     2     4     2    25
##  9   2018 Rick~ Red ~ AL    X         0     6     7    17    19    15    12
## 10   2018 Kyle~ Twins AL    X         0    17    22    15    12    17    10
## # ... with 516 more rows, and 7 more variables: k <int>, bb <int>, hr <int>,
## #   k9 <int>, bb9 <int>, kbb <int>, babip <int>


같은 이치로 시험용 자료 확인에는 testing() 함수를 씁니다.

cya_split %>%
  testing()
## # A tibble: 225 x 19
##    season name  team  lg    cy     vote     w     l   war   era   fip    ip
##     <int> <fct> <fct> <fct> <fct> <int> <int> <int> <int> <int> <int> <int>
##  1   2018 Carl~ Indi~ AL    X         0     4    14     6     8     4    23
##  2   2018 Jako~ Roya~ AL    X         0    22    21    22    21    22     5
##  3   2018 Jake~ Twins AL    X         0    24    17    16    23    19     8
##  4   2018 Char~ Astr~ AL    X         0    10     1    13     7    10    17
##  5   2018 Davi~ Red ~ AL    X         0     8     8    18    11    16    11
##  6   2018 Jame~ Whit~ AL    X         0    26    26    24    24    24     3
##  7   2018 Blak~ Rays  AL    O       169     1     2     7     1     6    22
##  8   2018 Jon ~ Rock~ NL    X         0    16     9    15    29    18    12
##  9   2018 Gio ~ - - - NL    X         0    20    24    19    23    20    10
## 10   2018 Tann~ Nati~ NL    X         0    23    29    18    25    23     7
## # ... with 215 more rows, and 7 more variables: k <int>, bb <int>, hr <int>,
## #   k9 <int>, bb9 <int>, kbb <int>, babip <int>


이걸로 데이터 분류가 끝났습니다. 이제 요리를 시작할 시간입니다.



recipe 하나만 있으면 요리 끝!

tidymodels는 이 기본 모형을 만드는 작업을 요리법에 비유한 recipes 패키지를 포함하고 있습니다.


일단 학습용 데이터로 기본 모형을 만드는 코드는 이렇게 쓸 수 있습니다. 

cya_split %>% training() %>%
  recipe(vote~w+l+war+era+fip+ip+k+bb+hr+k9+bb9+kbb+babip)
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13


결과는 투표 결과 하나고 13개 변수를 예측에 쓴다는 뜻입니다.


모형을 만들 때는 이렇게 공식만 입력하는 게 아니라 결측값을 제거하고, 이상점(outlier)을 처리하고, 다중공선성이 나타나는 변수를 제거하고, 척도를 조정해 각 자료를 정규화하고, 가변수(dummy 변수)를 만드는 작업 등이 필요합니다.


tidymodels에서는 step_*() 함수로 이 작업을 진행합니다. 그리고 prep() 함수로 실제 준비를 마칩니다.


우리는 step_corr() 함수로 상관관계가 지나치게 큰 변수를 제거하고, step_center()step_scale() 함수를 써서 평균을 0으로 하는 척도를 만들어 보겠습니다.


이를 코드로 쓰면 아래처럼 나타납니다.

cya_split %>% training() %>%
  recipe(vote~w+l+war+era+fip+ip+k+bb+hr+k9+bb9+kbb+babip) %>%
  step_corr(all_predictors()) %>%
  step_center(all_predictors(), -all_outcomes()) %>%
  step_scale(all_predictors(), -all_outcomes()) %>%
  prep() -> cya_recipe


어떻게 생겼는지 열어봐야겠죠?

cya_recipe
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13
## 
## Training data contained 526 data points and no missing data.
## 
## Operations:
## 
## Correlation filter removed fip, bb9, k [trained]
## Centering for w, l, war, era, ip, bb, hr, k9, kbb, babip [trained]
## Scaling for w, l, war, era, ip, bb, hr, k9, kbb, babip [trained]


13개 변수 가운데 △수비 영향을 제거한 평균자책점(FIP) 순위(fip) △9이닝당 볼넷 순위(bb) △탈삼진 순위(k)를 다른 기록과 상관관계가 너무 높아서 제거해야 한다는 사실을 알 수 있습니다.


그리고 나머지 10개 변수에 대해서는 중심화와 척도화가 필요합니다.


이제 실제로 이 작업을 데이터에 적용할 차례.


receipes 패키지에서 시험용 데이터 처리는 bake() 함수가 담당합니다.

cya_recipe %>%
  bake(cya_split %>% testing()) -> cya_testing
cya_testing
## # A tibble: 225 x 11
##          w      l     war     era     ip      bb     hr      k9    kbb    babip
##      <dbl>  <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>   <dbl>  <dbl>    <dbl>
##  1 -1.32   -0.483 -1.13   -0.998   0.245 -1.21   -0.828 -1.24   -1.41   0.331  
##  2  0.185   0.101  0.174   0.0792 -1.26  -1.12    0.255 -0.163  -0.835  0.00219
##  3  0.352  -0.233 -0.315   0.245  -1.01   0.137  -0.995 -0.329   0.154 -0.244  
##  4 -0.815  -1.57  -0.559  -1.08   -0.255 -0.0305 -1.33  -1.16   -0.175 -0.573  
##  5 -0.982  -0.984 -0.152  -0.749  -0.755 -0.786  -0.412 -0.578  -0.588 -0.819  
##  6  0.519   0.518  0.337   0.328  -1.42   0.305   0.422  0.253   0.401 -1.31   
##  7 -1.57   -1.48  -1.05   -1.58    0.162 -0.114  -1.50  -1.33   -0.423 -1.56   
##  8 -0.315  -0.901 -0.396   0.742  -0.672 -0.618   0.505 -1.08   -0.752  0.577  
##  9  0.0185  0.351 -0.0705  0.245  -0.839  0.473  -1.08  -0.0795  0.731  0.331  
## 10  0.269   0.769 -0.152   0.411  -1.09  -0.954   0.172  0.170  -0.340  0.166  
## # ... with 215 more rows, and 1 more variable: vote <int>


위에서 예고했던 대로 자료를 처리한 걸 알 수 있습니다.


이어서 학습용 데이터 처리. 이번에 쓸 함수는 juice()입니다. 

cya_recipe %>%
  juice() -> cya_training
cya_training
## # A tibble: 526 x 11
##          w       l     war      era      ip     bb       hr       k9     kbb
##      <dbl>   <dbl>   <dbl>    <dbl>   <dbl>  <dbl>    <dbl>    <dbl>   <dbl>
##  1 -0.648  -0.567  -0.885  -0.252   -0.0884 -1.63  -1.41     0.00363 -1.25  
##  2  0.269   0.435   0.418   0.411   -1.51   -0.534  0.505   -0.910   -0.340 
##  3 -1.40   -0.901  -1.21   -0.915    0.0783 -0.954 -1.08    -1.08    -1.33  
##  4 -0.148  -0.400   0.0110 -0.00362 -0.922  -1.46  -0.662    0.502   -0.505 
##  5 -0.482  -0.0660 -0.804  -0.501   -0.338  -0.366 -1.25     0.336    0.0718
##  6 -0.565  -1.23   -1.37   -1.49     0.495  -0.450 -1.58    -1.41    -0.917 
##  7  0.0185 -0.316  -0.478  -0.169   -0.589  -1.04   0.0886  -0.495   -1.08  
##  8 -0.898  -1.40   -1.46   -1.33     0.412  -0.198 -1.16    -1.58    -1.16  
##  9 -1.15   -1.07   -0.233  -0.0865  -0.672  -0.870  0.00523 -0.412   -1.00  
## 10 -0.232   0.184  -0.396  -0.666   -0.839   0.389 -0.745   -0.0795   0.237 
## # ... with 516 more rows, and 2 more variables: babip <dbl>, vote <int>


이제 정말 자료 준비가 끝이 났습니다. 이제 진짜 랜덤포레스트 모형을 적용할 차례입니다.



엔진 시동 부릉부릉~

tidymodels 특징 가운데 하나는 여러 패키지를 같은 코드로 처리할 수 있도록 만든다는 것.


R에서 랜덤 로레스트 작업을 진행할 때는 randomForest 패키지를 써도 되고 ranger 패키지를 써도 됩니다.


둘은 그냥 패키지만 다른 게 아니라 속성을 쓰는 방식도 다릅니다.


예를 들어 나무 숫자를 지정할 때 randomForest에서는 ntree라는 속성을 써야 하는데 ranger에서는 num.trees입니다.


tidymodels에 들어 있는 parsnip 패키지를 쓰면 그냥 ntrees 속성으로 두 패키지에 전부 원하는 나무 숫자를 지정할 수 있도록 해줍니다.


그저 set_engine() 안에 어떤 패키지를 쓸 것인지만 알려주면 됩니다.


랜덤 포레스트 모형은 크게 분류(classification)와 회귀(regression) 모드로 나눌 수 있습니다.


이번에 우리가 하려는 작업에는 회귀 모드가 어울립니다.


그러면 randomForest 패키지를 쓰는 코드는 이렇게 쓸 수 있습니다.

rand_forest(trees=100, mode='regression') %>%
  set_engine('randomForest') %>%
  fit(vote~w+l+war+era+ip+bb+hr+k9+kbb+babip, data=cya_training) -> cya_rf
cya_rf
## parsnip model object
## 
## Fit time:  130ms 
## 
## Call:
##  randomForest(x = as.data.frame(x), y = y, ntree = ~100) 
##                Type of random forest: regression
##                      Number of trees: 100
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 241.2098
##                     % Var explained: 82.57


ranger 패키지를 써서 작업을 진행하고 싶을 때는 set_engine() 부분만 바꿔주면 그만입니다.

rand_forest(trees=100, mode='regression') %>%
  set_engine('ranger') %>%
  fit(vote~w+l+war+era+ip+bb+hr+k9+kbb+babip, data=cya_training) -> cya_rg
cya_rg
## parsnip model object
## 
## Fit time:  31ms 
## Ranger result
## 
## Call:
##  ranger::ranger(formula = formula, data = data, num.trees = ~100,      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1)) 
## 
## Type:                             Regression 
## Number of trees:                  100 
## Sample size:                      526 
## Number of independent variables:  10 
## Mtry:                             3 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       247.2441 
## R squared (OOB):                  0.8216349


그저 '머신러닝으로 2019 사이영상 수상자 예상하기' 포스트를 쓸 때 randomForest 패키지를 썼다는 이유 하나만으로 앞으로는 cya_rf를 기준으로 작업을 진행하겠습니다.



랜덤 포레스팅 결과 확인

이렇게 통계 모형을 만들고 결과를 확인하고 싶을 때는 predict() 함수를 쓰면 됩니다.


우리는 랜덤 포레스트 모형 cya_rf를 시험용 데이터 그러니까 cya_testing에 적용하고 싶은 거니까 이렇게 쓰면 됩니다.

cya_rf %>%
  predict(cya_testing)
## # A tibble: 225 x 1
##        .pred
##        <dbl>
##  1  1.81e+ 1
##  2 -2.29e-14
##  3  3.60e- 2
##  4  1.06e+ 1
##  5  3.42e- 1
##  6  2.00e- 2
##  7  1.10e+ 2
##  8  3.96e- 1
##  9 -2.41e-14
## 10 -2.41e-14
## # ... with 215 more rows


결과만 나와서 헷갈립니다. dplyr 패키지에 있는 bind_cols() 함수를 써서 합쳐 보겠습니다.

cya_rf %>%
  predict(cya_testing) %>%
  bind_cols(cya_testing)
## # A tibble: 225 x 12
##        .pred       w      l     war     era     ip      bb     hr      k9    kbb
##        <dbl>   <dbl>  <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>   <dbl>  <dbl>
##  1  1.81e+ 1 -1.32   -0.483 -1.13   -0.998   0.245 -1.21   -0.828 -1.24   -1.41 
##  2 -2.29e-14  0.185   0.101  0.174   0.0792 -1.26  -1.12    0.255 -0.163  -0.835
##  3  3.60e- 2  0.352  -0.233 -0.315   0.245  -1.01   0.137  -0.995 -0.329   0.154
##  4  1.06e+ 1 -0.815  -1.57  -0.559  -1.08   -0.255 -0.0305 -1.33  -1.16   -0.175
##  5  3.42e- 1 -0.982  -0.984 -0.152  -0.749  -0.755 -0.786  -0.412 -0.578  -0.588
##  6  2.00e- 2  0.519   0.518  0.337   0.328  -1.42   0.305   0.422  0.253   0.401
##  7  1.10e+ 2 -1.57   -1.48  -1.05   -1.58    0.162 -0.114  -1.50  -1.33   -0.423
##  8  3.96e- 1 -0.315  -0.901 -0.396   0.742  -0.672 -0.618   0.505 -1.08   -0.752
##  9 -2.41e-14  0.0185  0.351 -0.0705  0.245  -0.839  0.473  -1.08  -0.0795  0.731
## 10 -2.41e-14  0.269   0.769 -0.152   0.411  -1.09  -0.954   0.172  0.170  -0.340
## # ... with 215 more rows, and 2 more variables: babip <dbl>, vote <int>


모형을 만들었으면 성능이 어느 정도 되는지 알아봐야겠죠?


이때는 yardstick 패키지에 있는 metrics() 함수를 씁니다. 

cya_rf %>%
  predict(cya_testing) %>%
  bind_cols(cya_testing) %>%
  metrics(truth=vote, estimate=.pred)
## # A tibble: 3 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard      16.1  
## 2 rsq     standard       0.798
## 3 mae     standard       6.15


결정계수(R²) .798이면 썩 쓸 만한 모형이라고 할 수 있습니다.


사실 우리에게 필요한 게 딱 떨어지는 값이 아니라는 점에서 더욱 그렇습니다.



사이영상 수상자를 예측해 보자

그럼 한 번 이 모형을 2018년 이전 전체 데이터에 적용해 사이영상 수상자를 예상해 보겠습니다.


일단 레시피를 새로 만들고

cya_pre %>%
  recipe(vote~w+l+war+era+fip+ip+k+bb+hr+k9+bb9+kbb+babip) %>%
  step_corr(all_predictors()) %>%
  step_center(all_predictors(), -all_outcomes()) %>%
  step_scale(all_predictors(), -all_outcomes()) %>%
  prep() -> cya_recipe2
cya_recipe2
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13
## 
## Training data contained 751 data points and no missing data.
## 
## Operations:
## 
## Correlation filter removed fip, bb9, k [trained]
## Centering for w, l, war, era, ip, bb, hr, k9, kbb, babip [trained]
## Scaling for w, l, war, era, ip, bb, hr, k9, kbb, babip [trained]


시험용 데이터는 굽고

cya_recipe2 %>%
  bake(cya_pre) -> cya_testing_pre


학습용 데이터에서는 즙을 짭니다.

cya_recipe2 %>%
  juice() -> cya_training_pre


그리고 새 데이터로 랜덤 포레스트 모형을 만듭니다.

rand_forest(trees=100, mode='regression') %>%
  set_engine('randomForest', localImp=TRUE) %>%
  fit(vote~w+l+war+era+ip+bb+hr+k9+kbb+babip, data=cya_training_pre) -> cya_rf2


원래 데이터와 열을 합쳐 보면 이런 결과가 나옵니다.

cya_rf2 %>%
  predict(cya_testing_pre) %>%
  bind_cols(cya_pre)
## # A tibble: 751 x 20
##        .pred season name  team  lg    cy     vote     w     l   war   era
##        <dbl>  <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  5.59e- 1   2018 Marc~ Mari~ AL    X         0    12    13     9    17
##  2 -3.73e-14   2018 Dyla~ Orio~ AL    X         0    23    25    25    25
##  3  1.14e+ 1   2018 Carl~ Indi~ AL    X         0     4    14     6     8
##  4  9.89e+ 0   2018 Luis~ Yank~ AL    X         1     3     9     5     9
##  5 -3.90e-14   2018 Mike~ Mari~ AL    X         0    18    15    20    20
##  6 -2.54e-14   2018 Dall~ Astr~ AL    X         0    14    19    10    14
##  7 -4.03e-14   2018 Jako~ Roya~ AL    X         0    22    21    22    21
##  8  3.43e+ 1   2018 Trev~ Indi~ AL    X        13    13     5     3     2
##  9 -2.87e-14   2018 Andr~ Ange~ AL    X         0    20    16    14    18
## 10 -4.24e-14   2018 Jake~ Twins AL    X         0    24    17    16    23
## # ... with 741 more rows, and 9 more variables: fip <dbl>, ip <dbl>,
## #   k <dbl>, bb <dbl>, hr <dbl>, k9 <dbl>, bb9 <dbl>, kbb <dbl>,
## #   babip <dbl>


변수가 너무 많아서 잘 안 보입니다. 


사실 우리에게 이제 저 순위 변수는 필요가 없습니다. 그리고 선수도 이렇게 많이 필요하지 않습니다.


그저 이 모형이 최고점을 예상한 선수가 실제로 사이영상을 탔는지 아닌지만 알면 됩니다.


이를 알아보려고 먼저 시즌과 리그에 따라 모형 예상값(.pred) 순위를 매겼습니다.


그리고 예상값이 1등이거나 아니면 사이영상을 탄, 즉 cy 변수에 O가 들어 있는 행만 골라냈습니다.


코드로 쓰면 이렇습니다.

cya_rf2 %>%
  predict(cya_testing_pre) %>%
  bind_cols(cya_pre) %>%
  group_by(season, lg) %>%
  mutate(rank=rank(-.pred)) %>%
  select(season, lg, team, name, cy, rank) %>%
  filter(rank==1 | cy=='O') %>%
  arrange(season, lg, cy)
## # A tibble: 22 x 6
## # Groups:   season, lg [20]
##    season lg    team      name             cy     rank
##     <int> <fct> <fct>     <fct>            <fct> <dbl>
##  1   2009 AL    Royals    Zack Greinke     O         1
##  2   2009 NL    Giants    Tim Lincecum     O         1
##  3   2010 AL    Mariners  Felix Hernandez  O         1
##  4   2010 NL    Phillies  Roy Halladay     O         1
##  5   2011 AL    Tigers    Justin Verlander O         1
##  6   2011 NL    Dodgers   Clayton Kershaw  O         1
##  7   2012 AL    Rays      David Price      O         2
##  8   2012 AL    Tigers    Justin Verlander X         1
##  9   2012 NL    Mets      R.A. Dickey      O         1
## 10   2013 AL    Tigers    Max Scherzer     O         1
## 11   2013 NL    Dodgers   Clayton Kershaw  O         1
## 12   2014 AL    Indians   Corey Kluber     O         1
## 13   2014 NL    Dodgers   Clayton Kershaw  O         1
## 14   2015 AL    Astros    Dallas Keuchel   O         1
## 15   2015 NL    Cubs      Jake Arrieta     O         1
## 16   2016 AL    Red Sox   Rick Porcello    O         2
## 17   2016 AL    Tigers    Justin Verlander X         1
## 18   2016 NL    Nationals Max Scherzer     O         1
## 19   2017 AL    Indians   Corey Kluber     O         1
## 20   2017 NL    Nationals Max Scherzer     O         1
## 21   2018 AL    Rays      Blake Snell      O         1
## 22   2018 NL    Mets      Jacob deGrom     O         1


그 결과 2012년과 2016년 아메리칸리그 사이영상 수상자가 이 모형이 예상한 것과 다른 선수가 받았다는 사실을 알 수 있습니다.


공교룝게도 두 번 모두 예상 점수 1위는 저스틴 벌랜더였습니다.


20번 중 18번을 맞췄으니까 적중률 90%입니다.


같은 모형을 2019년 자료에도 적용해 봐야겠죠?


먼저 시험용 데이터를 굽고 

cya_recipe2 %>%
  bake(cya_2019) -> cya_testing_2019


바로 모형을 적용해 어떤 선수가 상위권에 포진했는지 알아봅니다.

cya_rf2 %>%
  predict(cya_testing_2019) %>%
  bind_cols(cya_2019) %>%
  arrange(-.pred) %>%
  select(name, lg, team, .pred)
## # A tibble: 56 x 4
##    name              lg    team      .pred
##    <fct>             <fct> <fct>     <dbl>
##  1 Gerrit Cole       AL    Astros    157. 
##  2 Justin Verlander  AL    Astros    150. 
##  3 Jacob deGrom      NL    Mets      107. 
##  4 Stephen Strasburg NL    Nationals  88.6
##  5 Hyun-Jin Ryu      NL    Dodgers    73.8
##  6 Charlie Morton    AL    Rays       72.7
##  7 Max Scherzer      NL    Nationals  34.2
##  8 Shane Bieber      AL    Indians    33.5
##  9 Lance Lynn        AL    Rangers    30.9
## 10 Mike Soroka       NL    Braves     22.8
## # ... with 46 more rows


그러면 아메리칸리그에서는 게릿 콜, 내셔널리그에서는 제이컵 디그롬이 제일 높은 예상 점수를 받았다는 걸 알 수 있습니다.


안타깝게도(?) 지난해 아메리칸리그 사이영상은 콜이 아니라 예상 점수 2위 벌랜더가 받아 갔습니다.


이제 모형 예상과 실제 수상 점수 2-1입니다.



사이영상 수상에는 어떤 기록이 중요할까?

'머신러닝으로 2019 사이영상 수상자 예상하기' 포스트를 읽으신 분은 각 변수별 중요도를 나타낸 그래프를 보셨을 겁니다.


randomForestExpaliner 패키지를 설치하면 이런 결과를 손쉽게 확인할 수 있습니다.


패키지가 등장했으니 설치하고 불러오는 과정을 거쳐야 합니다.

#install.packages('randomForestExplainer')
library('randomForestExplainer')
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2


cya_rf2 모델을 적합(fit)하는 과정에서 중요한 기록을 뽑아내는 거니까 이런 코드를 쓰면 됩니다.

measure_importance(cya_rf2$fit)
##    variable mean_min_depth no_of_nodes mse_increase node_purity_increase
## 1     babip           2.85        1780    80.639925             52023.43
## 2        bb           3.55        1659    21.054700             20840.12
## 3       era           1.41        2014   549.119101            313366.48
## 4        hr           3.49        1610     6.607821             19655.34
## 5        ip           2.98        1764    28.909621             36859.15
## 6        k9           2.80        1774    24.114781             36035.92
## 7       kbb           2.49        1648    53.136692             54522.94
## 8         l           2.96        1661    15.834681             41244.26
## 9         w           1.24        2002   373.215213            273889.47
## 10      war           1.63        1813   233.223120            184095.73
##    no_of_trees times_a_root      p_value
## 1          100            2 4.291676e-01
## 2          100            0 9.980043e-01
## 3          100           25 1.587096e-09
## 4          100            0 9.999823e-01
## 5          100            6 5.879063e-01
## 6          100            4 4.886822e-01
## 7          100           10 9.992111e-01
## 8          100            5 9.976564e-01
## 9          100           28 8.867385e-09
## 10         100           20 1.582939e-01


이 가운데 저는 'node purity increase' 항목으로 중요도를 따졌습니다.


기왕이면 100점 만점을 기준으로 중요도를 나타낸다면 더 이해하기 쉽겠죠?

measure_importance(cya_rf2$fit) %>%
  as_tibble() %>%
  mutate(imp=node_purity_increase*100/max(node_purity_increase)) %>%
  arrange(-imp) %>%
  select(variable, imp)
## # A tibble: 10 x 2
##    variable    imp
##    <fct>     <dbl>
##  1 era      100   
##  2 w         87.4 
##  3 war       58.7 
##  4 kbb       17.4 
##  5 babip     16.6 
##  6 l         13.2 
##  7 ip        11.8 
##  8 k9        11.5 
##  9 bb         6.65
## 10 hr         6.27


이를 토대로 사이영상 투표권을 가진 이들이 평균자책점과 승리에 제일 큰 가중치를 둔다는 사실을 짐작할 수 있습니다. 


어떤가요? '머신러닝 작덥을 직접 해보시면 너무 간단하게 놀라실지 모른다'던 말씀을 거짓말이었는지 몰라도 완전 외계어 수준은 아니지 않나요?


그럼 모두들 Happy Tidyversing, Happy Tidymodeling!


댓글,

KindeR | 카테고리 다른 글 더 보기