聰明不如鈍筆
총명불여둔필
영어 낱말 'broom'은 아이콘에서 보시는 것처럼 '빗자루'라는 뜻. 빗자루는 지저분한 무엇인가를 쓸어담는 도구입니다.
그렇다면 이 패키지가 쓸어담는 건 무엇일까요? 네, 지저분한 통계 분석 결과입니다.
broom은 통계 분석 결과를 깔끔하게 데이터 테이블에 담아서 이후 작업을 원활하게 진행할 수 있도록 도와줍니다.
이번 포스트에서는 (이 블로그 전통에 따라) 프로야구 타격 기록을 가지고 broom 패키지를 써서 통계 모형을 깔끔하게(tidy) 정리하는 방법을 알아보도록 하겠습니다.
제일 먼저 할 일은 물론 tidyverse 패키지를 (설치하고) 불러오기.
#install.packages('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()
tidyverse를 불러 오면 패키지를 총 8개 패키지를 불러오는데 그 안에 broom은 없습니다. 그래서 따로 (설치하고) 불러오셔야 합니다.
#install.packages('broom')
library('broom')broom 패키지는 tidyverse가 아니라 tidymodels 안에 들어 있습니다. 원하신다면 이를 설치하고 불러오셔도 됩니다.
#install.packages('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 파일에 우리가 이번 포스트에서 쓸 프로야구 타격 데이터가 들어 있습니다.
tidyverse에는 원래 CSV 파일을 읽을 때 쓰는 read_csv() 함수가 들어 있습니다. 문제는 이 함수를 쓰면 한글 인코딩이 말썽을 일으킬 때가 있다는 것. 그래서 저는 주로 이런 방법으로 CSV 파일을 불러옵니다.
kbo.csvkbo <- read.csv('kbo.csv') %>% as_tibble()kbo## # A tibble: 303 x 7 ## 연도 구단 타석당득점 출루율 장타력 OPS wOBA ## <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1982 MBC 0.137 0.354 0.41 0.764 0.346 ## 2 1982 OB 0.129 0.349 0.412 0.762 0.343 ## 3 1982 롯데 0.115 0.344 0.373 0.717 0.329 ## 4 1982 삼미 0.102 0.304 0.345 0.648 0.298 ## 5 1982 삼성 0.141 0.347 0.392 0.739 0.338 ## 6 1982 해태 0.125 0.328 0.408 0.736 0.332 ## 7 1983 MBC 0.109 0.323 0.354 0.677 0.312 ## 8 1983 OB 0.111 0.324 0.362 0.687 0.314 ## 9 1983 롯데 0.0989 0.314 0.366 0.68 0.313 ## 10 1983 삼미 0.0923 0.312 0.344 0.656 0.302 ## # ... with 293 more rows
'%>%'라는 기호가 낯선 분도 혹시 계실지 모르겠습니다. 이 기호는 R에서 파이프(pipe)를 뜻합니다.
현실 세계에서 파이프가 기체나 액체 같은 유체를 한 곳에서 다른 곳으로 보내는 것처럼 코딩에서 파이프는 한 자료를 이 함수에서 저 함수로 보내는 구실을 합니다.
그러니까 저 코드는 kbo.csv 파일을 불러와서 이를 tidyverse에서 쓰는 데이터 프레임 형식인 tibble로 바꾼 다음 이를 kbo라는 변수에 저장하라는 뜻입니다.
이제 통계 모형을 만들 차례. OPS(출루율+장타력)가 타석당득점과 어떤 관계가 있는지 알려주는 (선형) 회귀분석 모형을 만들겠습니다.
R에서 회귀분석을 진행할 때는 lm() 함수를 쓰면 됩니다. 코드는 이렇게 씁니다.
ops_fit <- lm(타석당득점~OPS, kbo)ops_fit## ## Call: ## lm(formula = 타석당득점 ~ OPS, data = kbo) ## ## Coefficients: ## (Intercept) OPS ## -0.09663 0.29345
이렇게 변수 내용만 열어보는 걸로는 우리가 원하는 내용을 확인하기 어려울 때가 많습니다.
그래서 보통 summary() 함수를 써서 내용을 확인하는 일이 많습니다. 이렇게 말입니다.
ops_fit %>%
summary()## ## Call: ## lm(formula = 타석당득점 ~ OPS, data = kbo) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.0105880 -0.0033245 -0.0001155 0.0034264 0.0207466 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.096626 0.003787 -25.51 <2e-16 *** ## OPS 0.293449 0.005133 57.17 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.004776 on 301 degrees of freedom ## Multiple R-squared: 0.9157, Adjusted R-squared: 0.9154 ## F-statistic: 3268 on 1 and 301 DF, p-value: < 2.2e-16
확실히 아까보다는 많은 정보를 알 수 있게 됐습니다.
그런데 결과가 이렇게 나오면 그다음 작업을 진행할 때 손이 많이 갑니다. 만약 이런 결과가 데이터 프레임 안에 들어 있다면 후반 작업을 하기가 더 편할 겁니다.
broom이 하는 일이 바로 이런 결과를 데이터 프레임, 조금 더 정확하게는 tibble에 넣어주는 겁니다.
broom은 크게 세 가지 동사를 써서 이런 작업을 처리합니다.
먼저 모형 작업 결과를 정리하는 tidy(),
ops_fit %>%
tidy()## # A tibble: 2 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) -0.0966 0.00379 -25.5 3.10e- 77 ## 2 OPS 0.293 0.00513 57.2 1.11e-163
이어서 모형 성능을 정리해서 알려주는 glance(),
ops_fit %>%
glance()## # A tibble: 1 x 11 ## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC ## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> ## 1 0.916 0.915 0.00478 3268. 1.11e-163 2 1190. -2375. -2364. ## # ... with 2 more variables: deviance <dbl>, df.residual <int>
마지막으로 모형을 각 데이터에 적용하는 augment()입니다.
ops_fit %>%
augment()## # A tibble: 303 x 9 ## 타석당득점 OPS .fitted .se.fit .resid .hat .sigma .cooksd .std.resid ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 0.137 0.764 0.128 0.000310 0.00931 0.00421 0.00475 8.08e-3 1.95 ## 2 0.129 0.762 0.127 0.000305 0.00181 0.00409 0.00478 2.96e-4 0.380 ## 3 0.115 0.717 0.114 0.000291 0.00151 0.00371 0.00478 1.86e-4 0.316 ## 4 0.102 0.648 0.0935 0.000528 0.00871 0.0122 0.00476 2.08e-2 1.83 ## 5 0.141 0.739 0.120 0.000275 0.0207 0.00331 0.00463 3.14e-2 4.35 ## 6 0.125 0.736 0.119 0.000274 0.00573 0.00330 0.00477 2.39e-3 1.20 ## 7 0.109 0.677 0.102 0.000408 0.00698 0.00731 0.00477 7.92e-3 1.47 ## 8 0.111 0.687 0.105 0.000372 0.00602 0.00606 0.00477 4.87e-3 1.26 ## 9 0.0989 0.68 0.103 0.000397 -0.00399 0.00691 0.00478 2.44e-3 -0.838 ## 10 0.0923 0.656 0.0959 0.000493 -0.00358 0.0107 0.00478 3.07e-3 -0.754 ## # ... with 293 more rows
사실 모형이 하나일 때는 이런 동사가 무슨 도움이 될까 싶습니다. 모형이 두 개만 되어도 사정이 달라집니다.
한번 OPS와 가중출루율(wOBA) 가운데 어떤 기록이 타석당 득점을 더 잘 설명하는지 알아볼까요?
지금 우리가 쓰고 있는 데이터는 와이드 폼(wide form)입니다. 데이터 분석에는 롱 폼(long form)이 더 유리할 때가 많습니다.
이렇게 데이터 프레임 유형을 바꿀 때는 tidyr를 활용하면 됩니다.
자세한 내용이 궁금하신 분이 계시면 '최대한 친절하게 쓴 R로 데이터 깔끔하게 만들기(feat. tidyr)' 포스트가 도움이 될 수 있습니다.
현재 데이터에서 OPS와 wOBA만 길게 뽑아 보겠습니다.
kbo %>%
pivot_longer(cols=OPS:wOBA, names_to='기록', values_to='값')## # A tibble: 606 x 7 ## 연도 구단 타석당득점 출루율 장타력 기록 값 ## <int> <fct> <dbl> <dbl> <dbl> <chr> <dbl> ## 1 1982 MBC 0.137 0.354 0.41 OPS 0.764 ## 2 1982 MBC 0.137 0.354 0.41 wOBA 0.346 ## 3 1982 OB 0.129 0.349 0.412 OPS 0.762 ## 4 1982 OB 0.129 0.349 0.412 wOBA 0.343 ## 5 1982 롯데 0.115 0.344 0.373 OPS 0.717 ## 6 1982 롯데 0.115 0.344 0.373 wOBA 0.329 ## 7 1982 삼미 0.102 0.304 0.345 OPS 0.648 ## 8 1982 삼미 0.102 0.304 0.345 wOBA 0.298 ## 9 1982 삼성 0.141 0.347 0.392 OPS 0.739 ## 10 1982 삼성 0.141 0.347 0.392 wOBA 0.338 ## # ... with 596 more rows
pivot_longer()는 tidyr 업데이트로 새로 생긴 함수. 예전에는 gather() 함수로 같은 작업을 처리했습니다. 이렇게 말입니다.
kbo %>%
gather(key=기록, value=값, OPS:wOBA)## # A tibble: 606 x 7 ## 연도 구단 타석당득점 출루율 장타력 기록 값 ## <int> <fct> <dbl> <dbl> <dbl> <chr> <dbl> ## 1 1982 MBC 0.137 0.354 0.41 OPS 0.764 ## 2 1982 OB 0.129 0.349 0.412 OPS 0.762 ## 3 1982 롯데 0.115 0.344 0.373 OPS 0.717 ## 4 1982 삼미 0.102 0.304 0.345 OPS 0.648 ## 5 1982 삼성 0.141 0.347 0.392 OPS 0.739 ## 6 1982 해태 0.125 0.328 0.408 OPS 0.736 ## 7 1983 MBC 0.109 0.323 0.354 OPS 0.677 ## 8 1983 OB 0.111 0.324 0.362 OPS 0.687 ## 9 1983 롯데 0.0989 0.314 0.366 OPS 0.68 ## 10 1983 삼미 0.0923 0.312 0.344 OPS 0.656 ## # ... with 596 more rows
이렇게 롱 폼으로 바꾸고 나면 기록 그러니까 OPS와 wOBA를 기준으로 그룹을 나눠 분석할 수 있게 됩니다.
데이터를 그룹을 나눌 때는 dplyr 패키지게 들어 있는 group_by() 함수를 쓰면 됩니다.
'이게 무슨 소리냐?' 하시는 분께는 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
자 그럼 이제 cor.test() 함수를 써서 두 지표가 타석당 득점과 어떤 피어슨 상관계수를 나타내는지 분석해 보겠습니다.
kbo %>%
pivot_longer(cols=OPS:wOBA, names_to='기록', values_to='값') %>%
group_by(기록) %>%
summarise(cor.test(타석당득점, 값, method='pearson') %>%
tidy())## # A tibble: 2 x 9 ## # Groups: 기록 [2] ## 기록 estimate statistic p.value parameter conf.low conf.high method ## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr> ## 1 OPS 0.957 57.2 1.11e-163 301 0.946 0.965 Pears~ ## 2 wOBA 0.927 42.9 3.62e-130 301 0.909 0.941 Pears~ ## # ... with 1 more variable: alternative <chr>
어떤가요? 한번에 깔끔하게 결과가 잘 나왔습니다.
OPS(.957)가 wOBA(.927)보다 타석당 득점과 피어슨 상관계수가 더 높게 나타납니다. 양쪽 모두 (요즘 말 많고 탈 많은) 유의확률(p-value)도 유의미한 수준으로 낮습니다.
원래 하던 대로 회귀분석도 해볼까요?
이번에는 glance()로 열어서 R²를 비교하겠습니다.
kbo %>%
pivot_longer(cols=OPS:wOBA, names_to='기록', values_to='값') %>%
group_by(기록) %>%
summarise(lm(타석당득점~값) %>%
glance())## # A tibble: 2 x 12 ## # Groups: 기록 [2] ## 기록 r.squared adj.r.squared sigma statistic p.value df logLik AIC ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> ## 1 OPS 0.916 0.915 0.00478 3268. 1.11e-163 2 1190. -2375. ## 2 wOBA 0.859 0.859 0.00617 1837. 3.62e-130 2 1113. -2219. ## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>
이번에도 마찬가지로 OPS(.916)가 wOBA(.859)보다 타석당 득점을 더 잘 설명하는 것으로 나타났습니다
(개인적으로 제가 Sportugese나 기사에 OPS만 주구장창 쓰는 이유입니다.)
'R에서 반복 작업 깔끔하게 처리하기(feat. purrr)' 포스트를 읽어 보신 분은 여기서 배신감을 느끼실지도 모르겠습니다. 같은 작업을 훨씬 간단하게 처리했으니까요.
이게 바로 broom 패키지 힘입니다.
이렇게 OPS는 그저 출루율과 장타력을 더한 값일 뿐인데 이렇게 정확하게 타자 능력을 보여줍니다.
OPS가 타석당 득점을 얼마나 잘 설명하는지 그림으로 그려보겠습니다.
아래는 x축에 OPS, y축에 타석당 득점을 두고 산점도를 그리는 코드입니다.
kbo %>%
ggplot(aes(OPS, 타석당득점)) +
geom_point()이제 이 그래프에 모형 예측 결과를 넣을 차례.
가장 기본적인 방법은 predict() 함수를 쓰는 겁니다. 위에서 만든 ops_fit 예측 결과를 쓰면 되겠죠?
이 결과는 선 그래프로 표시합니다.
kbo %>%
ggplot(aes(OPS, 타석당득점)) +
geom_point() +
geom_line(aes(y=predict(ops_fit)))ggplot2에 익숙하신 분은 굳이 이렇게 따로 쓰지 않아도 된다는 걸 알고 계실 겁니다. stat_smooth() 함수가 있으니까요.
우리는 선형 모델을 기준으로 하고 있으니까 method='lm'을 지정하겠습니다.
kbo %>%
ggplot(aes(OPS, 타석당득점)) +
geom_point() +
stat_smooth(method='lm')stat_smooth() 함수를 쓰면 모형 예측값뿐 아니라 신뢰도 95%를 기준으로 표준 오차(standard error)를 같이 보여줍니다.
표준 오차가 존재한다는 건 이 그래프 기울기가 맞지 않을 수도 있다는 뜻. 이 차이를 줄일 수 있는 방법은 없을까요?
여기서 우리는 부트스트랩(bootstrap)을 활용할 겁니다. (이 블로그와 좋은 친구인) 위키피디아는 부트스트랩을 이렇게 설명합니다.
부트스트랩(bootstrapping)은 무작위 표본 추출에 의존하는 어떤 시험이나 계측이다. 부트스트랩은 표본 추정치들의 (편향, 분포, 신뢰 구간, 오차 예측 또는 기타 추정치들로 정의 되는) 정확도를 할당할 수 있도록 한다.
뭔가 알쏭달쏭합니다. 조금 더 자세하게 알아볼까요?
원래 부트스트랩은 부츠 달린 끈을 뜻하는 낱말.
여기서 "pull one's own by one's bootstrap"이라는 표현이 생겼습니다. 이 말은 '자기 혼자 노력으로 어떤 일을 해낸다'는 뜻입니다.
통계학에서 이 표현을 쓰는 건 표본을 통해 모집단 성질을 추정할 수 있는 것처럼 표본을 재표집(resampling)하면 표본 성질도 추정할 수 있기 때문입니다.
그러니까 샘플이 있을 때 그 샘플을 여러 번 (복원) 추출해서 평균이나 분산이 어떤 분포를 나타내는지 알아보는 과정이 부트스트랩입니다.
tidymodels 생태계에서는 rsample 패키지에 들어 있는 bootstraps() 함수가 부트스트랩을 담당합니다.
위에서 broom 패키지만 불러 오셨다면 아래처럼 패키지를 설치하고 불러오는 과정이 필요합니다.
#install.packages('rsample')
library('rsample')부트스트랩을 진행하는 과정은 간단합니다.
아래는 kbo 자료에서 100번 재표집을 진행하는 코드입니다.
kbo %>%
bootstraps(100)## # Bootstrap sampling ## # A tibble: 100 x 2 ## splits id ## <list> <chr> ## 1 <split [303/121]> Bootstrap001 ## 2 <split [303/116]> Bootstrap002 ## 3 <split [303/107]> Bootstrap003 ## 4 <split [303/109]> Bootstrap004 ## 5 <split [303/107]> Bootstrap005 ## 6 <split [303/110]> Bootstrap006 ## 7 <split [303/115]> Bootstrap007 ## 8 <split [303/112]> Bootstrap008 ## 9 <split [303/113]> Bootstrap009 ## 10 <split [303/108]> Bootstrap010 ## # ... with 90 more rows
그러면 이렇게 100번에 걸쳐 총 303개 표본 가운데 몇 개를 보였는지 알려주는 결과가 뜹니다. 오른쪽은 각 재표집 결과에 이름을 붙인 것.
지금은 그냥 샘플을 어떻게 다시 추출하겠다는 결과만 나와있을 뿐 분석을 진행한 건 아닙니다.
여러 재표집 결과에 분석을 진행하려면 반복 작업이 필요할 겁니다.
tidyverse 생태계에서는 purrr 패키지에 들어 있는 map() 계열 함수가 반복 작업을 담당합니다.
한번 간단하게 반복 작업을 진행해 보겠습니다.
먼저 1부터 5까지 숫자를 담고 있는 변수를 하나 만듭니다.
n <- c(1, 2, 3, 4, 5)n## [1] 1 2 3 4 5
map() 함수로 각 원소에 전부 1을 더하려면 이렇게 쓰면 됩니다.
map(n, function(x){x+1})## [[1]] ## [1] 2 ## ## [[2]] ## [1] 3 ## ## [[3]] ## [1] 4 ## ## [[4]] ## [1] 5 ## ## [[5]] ## [1] 6
결과가 리스트로 나오는 게 마음에 들지 않으시다면 map_dbl() 함수가 기다리고 있습니다.
map_dbl(n, function(x){x+1})## [1] 2 3 4 5 6
물론 함수에 이름을 지정하고 이를 활용하는 것도 가능합니다.
아래는 지금까지 우리가 쓴 함수를 add_one에 저장하고 이를 적용하는 코드입니다.
add_one <- function(x){x+1}map_dbl(n, add_one)## [1] 2 3 4 5 6
어떻게 하는지 감이 잡히시죠?
이런 반복 작업에 대해 좀 더 자세히 알고 싶으시다면 'R에서 반복 작업 깔끔하게 처리하기(feat. purrr)' 포스트가 도움이 될 수 있습니다.
우리가 하려는 건 회귀분석입니다.
먼저 OPS와 타석당 득점 사이에 회귀분석을 진행하는 함수를 만들어 lm_on_boots라고 이름을 붙이겠습니다.
lm_on_boots <- function(x) {
lm(타석당득점~OPS, data=x)
}그다음 재표집한 각 샘플(splits)에 이 함수를 적용하고 그 결과를 별도 열에 저장하는 코드를 쓰겠습니다.
이렇게 계산 결과를 데이터 프레임이 추가할 때는 dplyr 패키지에 있는 mutate() 함수를 쓰면 됩니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots))## # Bootstrap sampling ## # A tibble: 100 x 3 ## splits id model ## * <list> <chr> <list> ## 1 <split [303/111]> Bootstrap001 <lm> ## 2 <split [303/111]> Bootstrap002 <lm> ## 3 <split [303/112]> Bootstrap003 <lm> ## 4 <split [303/114]> Bootstrap004 <lm> ## 5 <split [303/108]> Bootstrap005 <lm> ## 6 <split [303/115]> Bootstrap006 <lm> ## 7 <split [303/115]> Bootstrap007 <lm> ## 8 <split [303/102]> Bootstrap008 <lm> ## 9 <split [303/117]> Bootstrap009 <lm> ## 10 <split [303/118]> Bootstrap010 <lm> ## # ... with 90 more rows
예, 오른쪽에 잘 들어왔습니다.
이렇게 mutate() 함수를 쓰는 게 어떤 개념인지 모르는 분이 계시다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
그러면 도대체 저 model이라는 열에는 어떤 데이터가 들어 있을까요?
한번 열어보겠습니다.
먼저 할 일은 위에서 우리가 본 결과는 kbo_boots라는 변수에 담는 것.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots)) -> kbo_boots그러고 나서 첫 번째 행을 열어보면 이렇게 lm() 함수를 적용한 걸 알 수 있습니다.
kbo_boots$model[1]## [[1]] ## ## Call: ## lm(formula = 타석당득점 ~ OPS, data = x) ## ## Coefficients: ## (Intercept) OPS ## -0.09218 0.28768
당연히 Bootstrap001부터 Bootstrap100까지 같은 작업을 진행했을 겁니다.
그러면 이 작업 결과를 원래 자료에 붙이려면 어떻게 해야 할까요? 네, augment() 함수를 쓰면 됩니다.
여기까지 작업한 걸 다시 kbo_boots에 넣어 보겠습니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots),
augmented=map(model, augment)) -> kbo_bootskbo_boots## # Bootstrap sampling ## # A tibble: 100 x 4 ## splits id model augmented ## * <list> <chr> <list> <list> ## 1 <split [303/106]> Bootstrap001 <lm> <tibble [303 x 9]> ## 2 <split [303/107]> Bootstrap002 <lm> <tibble [303 x 9]> ## 3 <split [303/108]> Bootstrap003 <lm> <tibble [303 x 9]> ## 4 <split [303/112]> Bootstrap004 <lm> <tibble [303 x 9]> ## 5 <split [303/101]> Bootstrap005 <lm> <tibble [303 x 9]> ## 6 <split [303/105]> Bootstrap006 <lm> <tibble [303 x 9]> ## 7 <split [303/116]> Bootstrap007 <lm> <tibble [303 x 9]> ## 8 <split [303/114]> Bootstrap008 <lm> <tibble [303 x 9]> ## 9 <split [303/112]> Bootstrap009 <lm> <tibble [303 x 9]> ## 10 <split [303/120]> Bootstrap010 <lm> <tibble [303 x 9]> ## # ... with 90 more rows
augmented 열도 열어보지 않을 수가 없겠죠?
kbo_boots$augmented[1]## [[1]] ## # A tibble: 303 x 9 ## 타석당득점 OPS .fitted .se.fit .resid .hat .sigma .cooksd .std.resid ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 0.0995 0.67 0.100 0.000440 -7.46e-4 0.00795 0.00495 9.20e-5 -0.152 ## 2 0.0912 0.649 0.0942 0.000534 -3.00e-3 0.0117 0.00494 2.21e-3 -0.611 ## 3 0.126 0.773 0.130 0.000364 -3.65e-3 0.00543 0.00494 1.50e-3 -0.742 ## 4 0.129 0.771 0.129 0.000357 -3.34e-4 0.00523 0.00495 1.21e-5 -0.0679 ## 5 0.104 0.694 0.107 0.000350 -2.90e-3 0.00503 0.00494 8.77e-4 -0.589 ## 6 0.117 0.738 0.120 0.000286 -3.07e-3 0.00335 0.00494 6.53e-4 -0.623 ## 7 0.0755 0.595 0.0788 0.000800 -3.30e-3 0.0262 0.00494 6.17e-3 -0.677 ## 8 0.117 0.744 0.121 0.000292 -4.33e-3 0.00349 0.00494 1.35e-3 -0.879 ## 9 0.0904 0.639 0.0914 0.000581 -9.99e-4 0.0138 0.00495 2.91e-4 -0.204 ## 10 0.113 0.701 0.109 0.000329 3.47e-3 0.00444 0.00494 1.10e-3 0.704 ## # ... with 293 more rows
지금은 공부 차원에서 '속살'이 어떻게 생겼는지 알아보려고 일부러 열어보고 있습니다.
위처럼 데이터 프레임 열 안에 리스트 또는 데이터 프레임이 들어 있는 형태를 '중첩(nested) 데이터 프레임'이라고 합니다.
이 프레임은 unnest() 함수로 열어서 내용을 확인하고 작업을 진행할 수 있습니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots),
augmented=map(model, augment)) %>%
unnest(augmented)## # A tibble: 30,300 x 12 ## splits id model 타석당득점 OPS .fitted .se.fit .resid .hat .sigma ## &llt;list> <chr> <lis> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 <spli~ Boot~ <lm> 0.125 0.768 0.128 3.26e-4 -3.76e-3 0.00494 0.00464 ## 2 <spli~ Boot~ <lm> 0.101 0.683 0.103 3.73e-4 -2.05e-3 0.00647 0.00464 ## 3 <spli~ Boot~ <lm> 0.146 0.829 0.147 5.75e-4 -1.02e-3 0.0154 0.00464 ## 4 <spli~ Boot~ <lm> 0.135 0.782 0.133 3.73e-4 2.66e-3 0.00648 0.00464 ## 5 <spli~ Boot~ <lm> 0.117 0.738 0.120 2.68e-4 -2.88e-3 0.00334 0.00464 ## 6 <spli~ Boot~ <lm> 0.117 0.744 0.121 2.73e-4 -4.21e-3 0.00347 0.00464 ## 7 <spli~ Boot~ <lm> 0.136 0.829 0.147 5.75e-4 -1.05e-2 0.0154 0.00460 ## 8 <spli~ Boot~ <lm> 0.107 0.711 0.111 2.89e-4 -4.33e-3 0.00390 0.00464 ## 9 <spli~ Boot~ <lm> 0.135 0.787 0.134 3.92e-4 5.04e-4 0.00716 0.00464 ## 10 <spli~ Boot~ <lm> 0.120 0.738 0.120 2.68e-4 6.28e-5 0.00334 0.00464 ## # ... with 30,290 more rows, and 2 more variables: .cooksd <dbl>, ## # .std.resid <dbl>
여기까지 왔으면 이제 그림을 그리는 건 어렵지 않습니다.
위에서 prediction() 함수로 예상한 부분에 해당하는 결과는 샘플별 .fitted 열에 들어 있습니다. 선이 100개나 되니까 투명도(alpha)를 조절해 그리겠습니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots),
augmented=map(model, augment)) %>%
unnest(augmented) %>%
ggplot(aes(OPS, 타석당득점)) +
geom_point() +
geom_line(aes(y=.fitted, group=id), alpha=.05)예, 아주 잘 나왔습니다.
위에서 보신 것처럼 이 선은 재표집한 샘플 100개를 전부 그린 겁니다. 우리는 세련되게(?) 신뢰 수준 95%만 남겨놓겠습니다.
이를 알아볼 때는 tidy() 함수를 써서 추정치(estimate)를 알아내고 quantile() 함수로 신뢰도 95%에 해당하는 부분만 남겨 놓으면 됩니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots),
tidy=map(model, tidy)) %>%
unnest(tidy) %>%
group_by(term) %>%
summarise(low=quantile(estimate, .05/2),
median=median(estimate),
high=quantile(estimate, 1-.05/2))## # A tibble: 2 x 4 ## term low median high ## <chr> <dbl> <dbl> <dbl> ## 1 (Intercept) -0.103 -0.0966 -0.0897 ## 2 OPS 0.284 0.293 0.302
이러면 다 끝났습니다. 이제 이 범위를 골라서 그래프를 그리면 됩니다.
kbo %>%
bootstraps(100) %>%
mutate(model=map(splits, lm_on_boots),
tidy=map(model, tidy),
augmented=map(model, augment)) %>%
unnest(tidy) %>%
group_by(term) %>%
mutate(low=quantile(estimate, .05/2),
median=median(estimate),
high=quantile(estimate, 1-.05/2)) %>%
filter(estimate>=low & estimate<=high) %>%
unnest(augmented) %>%
ggplot(aes(OPS, 타석당득점)) +
geom_point() +
geom_line(aes(y=.fitted, group=id), alpha=.05)타석당 득점~OPS 모델은 워낙 잘 떨어져서 이렇게 모델을 비교하는 게 별로 재미있지는 않습니다.
그렇다고 부트스트랩을 통해 신뢰 구간을 구하는 게 의미가 없는 작업은 아니었으리라고 생각합니다.
'이런 작업이 필요해서 이 포스트를 찾아오신 분이 계시다'에 500원 겁니다.
이제 마지막으로 출루율, 장타력 데이터를 가지고 k-평균 군집 분석을 해보겠습니다.
역시 이 블로그와 아주 좋은 친구 위키피디아에서 인용하면:
k-평균 알고리즘(K-means algorithm)은 주어진 데이터를 k개의 클러스터로 묶는 알고리즘으로, 각 클러스터와 거리 차이의 분산을 최소화하는 방식으로 동작한다.
kbo %>% select(출루율, 장타력) -> obpslg그다음 작업은 k를 1~9개 그룹으로 나누는 k-평균 군집 분석 모형을 만들고 glance() 함수로 모형 성능을 알아보겠습니다.
map(1:9, ~kmeans(obpslg, centers=.x)) %>%
map(glance)## [[1]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.539 -1.11e-16 1 ## ## [[2]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.182 0.357 1 ## ## [[3]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.111 0.428 3 ## ## [[4]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0752 0.464 2 ## ## [[5]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0603 0.479 3 ## ## [[6]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0496 0.489 5 ## ## [[7]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0439 0.495 3 ## ## [[8]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0407 0.498 3 ## ## [[9]] ## # A tibble: 1 x 4 ## totss tot.withinss betweenss iter ## <dbl> <dbl> <dbl> <int> ## 1 0.539 0.0368 0.502 4
여기서 우리가 눈여겨 봐야 하는 숫자는 'tot.withnss'입니다. 한번 이 지표만 따로 뽑아 볼까요?
map(1:9, ~kmeans(obpslg, centers=.x)) %>%
map(glance) %>%
map('tot.withinss')## [[1]] ## [1] 0.5388278 ## ## [[2]] ## [1] 0.1823147 ## ## [[3]] ## [1] 0.1112682 ## ## [[4]] ## [1] 0.07515769 ## ## [[5]] ## [1] 0.05996553 ## ## [[6]] ## [1] 0.0494076 ## ## [[7]] ## [1] 0.0456507 ## ## [[8]] ## [1] 0.04048165 ## ## [[9]] ## [1] 0.03730518
여기서 withnss는 '응집도'를 나타냅니다. 응집도는 각 그룹 중심에서 각 원소가 얼마나 떨어져 있는지를 알려주는 값입니다. 당연히 이 숫자가 적을수록 군집화에 성공했다고 할 수 있습니다.
'tot.withnss'는 그룹별 응집도를 모두 더한 값입니다. 그리고 k값에 따라 이 응집도 합이 어떻게 변하는지 살펴 보면 k를 몇 개로 정해야 하는지 알 수 있습니다.
x축에 그룹 숫자를 넣고 y축에 tot.withnss를 넣고 선 그래프를 그려 보면 아래 같은 모양이 나타납니다.
map(1:9, ~kmeans(obpslg, centers=.x)) %>%
map(glance) %>%
set_names(1:9) %>%
map_df('tot.withinss') %>%
pivot_longer(everything()) %>%
ggplot(aes(as.numeric(name), value)) +
geom_line()이렇게 보면 k가 3개를 넘어가면 응집도 합이 별 변화가 없는 걸 알 수 있습니다. 그러면 3개가 적당한 k 개수입니다. 이런 식으로 개수를 정하는 걸 흔히 '팔꿈치 기법'이라고 부릅니다.
k를 3개로 놓고 각 원소가 어떤 그룹에 속하는지 augment() 함수까지 파이프로 연결하면 아래 같은 결과를 얻을 수 있습니다.
obpslg %>%
kmeans(., centers=3) %>%
augment(obpslg)## # A tibble: 303 x 3 ## 출루율 장타력 .cluster ## <dbl> <dbl> <fct> ## 1 0.354 0.41 1 ## 2 0.349 0.412 1 ## 3 0.344 0.373 1 ## 4 0.304 0.345 2 ## 5 0.347 0.392 1 ## 6 0.328 0.408 1 ## 7 0.323 0.354 2 ## 8 0.324 0.362 2 ## 9 0.314 0.366 2 ## 10 0.312 0.344 2 ## # ... with 293 more rows
그림으로 잘 나왔는지도 확인해야겠죠? 전혀 어려울 게 없습니다.
obpslg %>%
kmeans(., centers=3) %>%
augment(obpslg) %>%
ggplot(aes(출루율, 장타력, color=.cluster)) +
geom_point()냉정하게 말해서 broom는 tidyverse 그리고 tidymodels 생태계에서만 필요합니다.
그런데 이 생태계가 곧 R 생태계라고 할 수 있을 만큼 번식력(?)이 왕성하기 때문에 broom을 쓰면 통계 모형 정리가 아주 손쉬워지는 것도 사실입니다.
이리 저리 모형 분석 결과를 정리하느라 헤매지 마시고 broom을 써보시기를 추천합니다.
언제나 그렇듯 잘못된 또는 이해가 잘 가지 않는 내용이 있거나 같은 작업을 처리하는 더 좋은 방법을 알고 계시면 댓글로 남겨주시면 대단히 감사하겠습니다.
그러면 여러분 모두 Happy Tidyversing, Happy Tidymodeling!
댓글,
