聰明不如鈍筆
총명불여둔필
원래 tidyverse 코어(core) 패키지를 처음 블로그에 소개할 때는 '최대한 친절하게 쓴'을 붙이는 게 원칙. 하지만 purrr 패키지까지 찾아오실 정보면 이미 R에 친숙하신 분이라고 생각해 붙이지 않습니다.
원래 영어 낱말 'purr'(/pɜːr/) 또는 'purrr'는 고양이가 가르랑 소리를 낸다는 뜻입니다. 이것만 보면 purrr 패키지가 어떤 구실을 하는지 짐작하기가 쉽지 않습니다.
게다가 tidyverse 홈페이지조차 purrr 패키지를 소개하면서 '함수용 프로그램 도구 모음'(functional programming toolkit)라는 알쏭달쏭한 표현을 씁니다.
힌트는 해들리 위컴 박사가 쓴 그 유명한 책 'R를 활용한 데이터 과학'에 들어 있습니다. 이 책에서 purrr 패키지를 다루고 있는 꼭지 제목은 'Iteration' 그러니까 반복입니다.
원래 코딩 세계에서 반복을 상징하는 명령은 'for'라고 할 수 있습니다. 반복 작업에 purrr를 쓰면 for를 쓸 때와 어떤 차이가 있을까요? 지금부터 한 번 알아보겠습니다.
시작은 물론 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개 패키지를 동시에 탑재하는데 그 안에 purrr도 들어 있습니다.
자료도 불러와야겠죠? 이번에도 (이 블로그 전통에 따라) 한국 프로야구 타격 데이터를 불러오도록 하겠습니다.
위에 있는 CSV 파일을 내려 받으시고 아래 명령어로 불러오시면 됩니다.
kbo.csvkbo <- read.csv('kbo.csv') %>% as_tibble()원래 tidyverse에 있는 readr에는 CSV 파일을 읽을 때 쓰는 read_csv() 함수가 들어 있습니다. 애석하게도 이 함수를 쓰면 한글 인코딩에 문제가 생기는 일이 잦아서 저는 보통 저런 식으로 자료를 읽으라고 말씀드립니다.
노파심에 말씀드리면 여기 등장한 '%>%' 기호는 파이프(pipe)를 뜻합니다. 현실 세계에서 파이프가 한 곳에서 다른 곳으로 액체 또는 기체를 보낼 수 있게 도와주는 것처럼 이 녀석도 한 함수에서 다른 함수로 자료를 보내주는 구실을 합니다.
위에 쓴 코드는 아래처럼 나누어 쓸 수 있습니다.
kbo <- read.csv('kbo.csv')
kbo <- as_tibble(kbo)(혹시 파이프가 더 궁금하신 분이 계시면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.)
자료를 불러 왔으니 어떻게 생겼는지도 확인해야겠죠? 그냥 변수 이름만 치면 됩니다.
kbo## # A tibble: 303 x 8 ## 연도 X10년대 구단 타석당득점 타율 출루율 장타력 OPS ## <int> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.282 0.354 0.41 0.764 ## 2 1982 a 삼성 0.141 0.266 0.347 0.392 0.739 ## 3 1982 a OB 0.129 0.283 0.349 0.412 0.762 ## 4 1982 a 해태 0.125 0.261 0.328 0.408 0.736 ## 5 1982 a 롯데 0.115 0.256 0.344 0.373 0.717 ## 6 1982 a 삼미 0.102 0.24 0.304 0.345 0.648 ## 7 1983 a 삼성 0.116 0.263 0.333 0.393 0.726 ## 8 1983 a 해태 0.113 0.267 0.332 0.385 0.717 ## 9 1983 a MBC 0.109 0.256 0.323 0.354 0.677 ## 10 1983 a OB 0.111 0.259 0.324 0.362 0.687 ## # ... with 293 more rows
나머지는 아실 테고 'X10년대'라는 열에는 1982~1989년 'a', 1990~1999년 'b', 2000~2009년 'c', 2010~2019년 'd'를 각각 넣었습니다. 열 이름 앞에 X가 붙은 건 R는 숫자로 시작하는 변수 이름을 받아들이지 않기 때문에 R가 알아서 수정한 겁니다.
우리 목표는 이 자료를 가지고 아래 같은 결과를 출력하는 코드를 만드는 것. 아래는 타율, 출루율, 장타력, OPS(출루율+장타력) 변화로 타석당 득점 변화를 얼마나 설명할 수 있는지 나타내는 결정계수(R²)입니다.
## # A tibble: 1 x 4 ## 타율 출루율 장타력 OPS ## <dbl> <dbl> <dbl> <dbl> ## 1 0.776 0.770 0.876 0.916
R²를 구하는 첫 단계는 선형 회귀분석입니다. R에서는 lm() 함수를 써서 회귀분석을 진행할 수 있습니다. 예를 들어 타율과 타석당 득점을 가지고 회귀분석을 하는 코드는 아래처럼 쓸 수 있습니다.
lm(타석당득점~타율, kbo)## ## Call: ## lm(formula = 타석당득점 ~ 타율, data = kbo) ## ## Coefficients: ## (Intercept) 타율 ## -0.1259 0.9222
이 결과에는 우리가 원하는 R²가 없습니다. 이를 구하려면 summary() 함수를 한 번 더 써야 합니다.
lm(타석당득점~타율, kbo) %>%
summary()## ## Call: ## lm(formula = 타석당득점 ~ 타율, data = kbo) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.0257544 -0.0053930 0.0005455 0.0042343 0.0216233 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.125931 0.007615 -16.54 <2e-16 *** ## 타율 0.922209 0.028581 32.27 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.007795 on 301 degrees of freedom ## Multiple R-squared: 0.7757, Adjusted R-squared: 0.775 ## F-statistic: 1041 on 1 and 301 DF, p-value: < 2.2e-16
'(Multiple) R-squared' 보이시죠? 이 값만 뽑아내는 코드는 아래처럼 쓸 수 있습니다.
lm(타석당득점~타율, kbo) %>%
summary() %>%
.$r.squared## [1] 0.7757327
출루율, 장타력, OPS에 대해서도 같은 계산을 하려면 아래 같은 코드를 반복하면 됩니다.
lm(타석당득점~출루율, kbo) %>%
summary() %>%
.$r.squaredlm(타석당득점~장타력, kbo) %>%
summary() %>%
.$r.squaredlm(타석당득점~OPS, kbo) %>%
summary() %>%
.$r.squared물론 우리는 이렇게 하지 않을 겁니다. '반복문'을 공부하는 포스트니까요.
아래처럼 1부터 5까지 들어 있는 변수가 있다고 할 때 각 항목에 1을 반복해 더하려면 어떻게 하면 될까요?
a <- c(1, 2, 3, 4, 5)위에서 말씀드린 것처럼 제일 기본은 for() 함수를 쓰는 겁니다. 이렇게 말입니다.
for(i in 1:5){
print(a[i]+1)
}## [1] 2 ## [1] 3 ## [1] 4 ## [1] 5 ## [1] 6
잘 나왔습니다. 단, for()는 변수(what)에 집중하기 때문에 우리가 하려는 동작(how)을 소홀히 한다는 문제가 있습니다.
맞습니다. R에 조금 익숙하신 분이라면 apply() 계열 함수를 쓰면 된다고 생각하셨을 겁니다. apply() 계열 함수는 *apply(데이터, 함수) 형태로 씁니다. 이렇게 말입니다.
lapply(a, function(x){x+1})## [[1]] ## [1] 2 ## ## [[2]] ## [1] 3 ## ## [[3]] ## [1] 4 ## ## [[4]] ## [1] 5 ## ## [[5]] ## [1] 6
여기서 function(x){x+1} 은 어떤 입력값(x)이 들어오면 그 값에 1을 더하라는 뜻입니다. 여기서는 a에 들어 있는 각 항목이 순환하면서 반복적으로 x가 됩니다.
purrr에 들어 있는 map() 함수도 기본적으로 같은 방식으로 씁니다.
map(a, function(x){x+1})## [[1]] ## [1] 2 ## ## [[2]] ## [1] 3 ## ## [[3]] ## [1] 4 ## ## [[4]] ## [1] 5 ## ## [[5]] ## [1] 6
완전히 똑같은 결과가 나왔습니다.
apply() 계열 함수는 어떤 자료 형태를 받아 어떤 형태로 출력하는지에 따라 종류를 구분합니다. 위에서 본 lapply()는 리스트 형태로 자료를 출력하는 함수입니다.
sapply() 함수를 쓰면 벡터 형태로 결과를 나타냅니다.
sapply(a, function(x){x+1})## [1] 2 3 4 5 6
purrr에서는 이런 차이가 '_' 뒤에 붙습니다. sapply()와 같은 결과를 만들고 싶을 때는 map_dbl()을 쓰면 됩니다.
map_dbl(a, function(x){x+1})## [1] 2 3 4 5 6
함수를 아예 바깥에 두는 방법도 있습니다.
우리가 지금껏 쓴 함수를 'add_one'이라는 이름으로 저장해 보겠습니다.
add_one <- function(x){x+1}그리고 나서 함수가 들어가는 자리에 'add_one'만 쓰면 같은 결과가 나옵니다.
map_dbl(a, add_one)## [1] 2 3 4 5 6
네, 네, apply() 계열 함수도 완전히 같은 방식으로 쓸 수 있습니다.
sapply(a, add_one)## [1] 2 3 4 5 6
이렇게 보면 apply() 계열과 map() 계열이 별 차이가 없는 것 같기도 합니다.
아닙니다. purrr를 쓸 때 첫 번째 장점은 함수를 알아보기 쉽게 줄여 쓸 수 있다는 것. 만약 function(x){x+1} 을 그냥 x+1 처럼 쓰면 무슨 뜻인지 더 알아보기 쉽습니다.
애석하게도 purrr에서도 이렇게 간단하게 되는 않습니다. purrr에서는 '여기서부터 함수야'하고 알려주려는 목적으로 ~ 기호를 씁니다. 그리고 x는 .x 가 됩니다.
따라서 function(x){x+1} 은 ~.x+1 처럼 쓸 수 있습니다. 그래서 아래처럼 쓰면 위와 같은 결과가 나옵니다.
map_dbl(a, ~.x+1)## [1] 2 3 4 5 6
이 작은 차이가 얼마나 큰 차이를 만들어 내는지 이 포스트를 통해 차차 확인해 보겠습니다.
여기까지 따라오시느라 고생하셨습니다. 잠시 쉬었다 가셔도 좋습니다.
map() 계열 함수에 어떤 종류가 있고 어떤 기능을 하는지 예제를 가지고 활용해 보겠습니다.
아래처럼 b에 서로 (문자) 길이가 다른 숫자를 넣었습니다.
b <- c(1, 22, 333, 4444, 55555)각 숫자 길이가 알고 싶습니다. 일단 문자 길이는 str_length() 함수로 알 수 있습니다.
문자 길이는 정수(integer)로 나오겠죠? 정수 벡터로 결과물을 내놓고 싶을 때는 map_int() 함수를 쓰시면 됩니다. 아래처럼 말입니다.
map_int(b, str_length)## [1] 1 2 3 4 5
지금까지 제일 많이 썼던 map_dbl은 '배정 밀도(double precision) 부동 소수점' 방식으로 숫자를 표시합니다. 말이 어려운데 결국 같은 결과를 내놓습니다.
map_dbl(b, str_length)## [1] 1 2 3 4 5
논리적(locigal)으로 어떤 값이 참(True)인지 거짓인지(False) 알려주는 함수도 있겠죠? 네, map_lgl()입니다. b에 들어 있는 값이 전부 숫자가 맞는지(is.numeric) 물어보면 전부 그렇다고 답이 나옵니다.
map_lgl(b, is.numeric)## [1] TRUE TRUE TRUE TRUE TRUE
그렇다고 map() 계열에서 숫자만 다룰 수 있는 건 아닙니다. 문자(character)도 물론 가능합니다.
변수 c에 문자를 넣은 다음 공백 없이(paste0()) 'z'를 붙이라고 한 번 써보겠습니다. 문자를 다룰 때는 map_chr()를 씁니다.
c <- c('abc', 'def', 'ghi')map_chr(c, ~paste0(.x, 'z'))## [1] "abcz" "defz" "ghiz"
지금까지는 전부 변수 하나만 가지고 계산했습니다. 작업 대상이 변수 두 개일 때는 어떻게 해야 할까요? 정답은 map2_*()입니다.
d에 a와 반대 방향으로 숫자를 넣은 다음 이 둘을 더하는 코드를 써보겠습니다.
d <- c(5, 4, 3, 2, 1)map2(a, d, sum)## [[1]] ## [1] 6 ## ## [[2]] ## [1] 6 ## ## [[3]] ## [1] 6 ## ## [[4]] ## [1] 6 ## ## [[5]] ## [1] 6
map()만 썼을 때처럼 결과가 리스트로 나왔습니다. 벡터 형태로 나타내고 싶으시면 그냥 map2_dbl()을 쓰면 그만입니다.
map2_dbl(a, d, sum)## [1] 6 6 6 6 6
함수 안에 들어가는 변수 x를 purrr에서는 .x 로 쓴다는 것 기억하고 계시죠? 변수가 두 개일 때는 두 번째를 .y 로 쓰시면 그만입니다. 그래서 아래처럼 쓰면 위와 똑같은 결과가 나옵니다.
map2_dbl(a, d, ~.x+.y)## [1] 6 6 6 6 6
자료가 세 개 이상일 때는 pmap_*()을 쓰시면 되고 이때는 .x , .y , .z 가 아니라 ..1 , ..2 , ..3 처럼 숫자로 구분합니다.
그래서 아래는 각 자리별로 b-a+d 를 계산하는 코드가 되겠습니다.
pmap_dbl(list(a, b, d), ~..2-..1+..3)## [1] 5 24 333 4442 55551
데이터를 계산할 때뿐 아니라 데이터를 골라낼 때도 map() 계열 함수를 써야 할 일이 있습니다. 먼저 아래처럼 리스트로 만든 리스트를 준비합니다.
e <- list(
list(-1, x=1, y=c(1), z='a'),
list(-2, x=4, y=c(2, 3), z='b'),
list(-3, x=9, y=c(4, 5, 6))
)e## [[1]] ## [[1]][[1]] ## [1] -1 ## ## [[1]]$x ## [1] 1 ## ## [[1]]$y ## [1] 1 ## ## [[1]]$z ## [1] "a" ## ## ## [[2]] ## [[2]][[1]] ## [1] -2 ## ## [[2]]$x ## [1] 4 ## ## [[2]]$y ## [1] 2 3 ## ## [[2]]$z ## [1] "b" ## ## ## [[3]] ## [[3]][[1]] ## [1] -3 ## ## [[3]]$x ## [1] 9 ## ## [[3]]$y ## [1] 4 5 6
여기서 맨 처음에 있는 데이터 그러니까 -1, -2, -3을 골라내고 싶을 때는 예상하시는 것처럼 함수 대신 숫자 1을 지정하면 됩니다.
map_dbl(e, 1)## [1] -1 -2 -3
그러면 두 번째를 골라낼 때는 2가 되겠죠?
map_dbl(e, 2)## [1] 1 4 9
두 번째 항목에는 'x'라는 이름이 붙어 있습니다. 이 이름을 써서 데이터를 골라낼 수도 있습니다.
map_dbl(e, 'x')## [1] 1 4 9
세 번째 항목에도 'y'라는 데이터가 붙어 있지만 데이터 개수가 달라서 벡터로 출력하려고 하면 에러 메시지가 나옵니다.
map_dbl(e, 'y')## Error: Result 2 must be a single double, not a double vector of length 2
이 때는 map()을 써서 리스트로 받아야 합니다.
map(e, 'y')## [[1]] ## [1] 1 ## ## [[2]] ## [1] 2 3 ## ## [[3]] ## [1] 4 5 6
y 가운데 첫 번째 값(1, 2, 4)만 뽑아 내고 싶을 때는 아래처럼 쓰면 됩니다.
map(e, list('y', 1))## [[1]] ## [1] 1 ## ## [[2]] ## [1] 2 ## ## [[3]] ## [1] 4
'z' 역시 세 번째 리스트에는 아예 자료가 없습니다. 그래서 이런 식으로 명령을 입력하면 결과를 볼 수 없습니다.
map_chr(e, 'z')## Error: Result 3 must be a single string, not NULL of length 0
단, .default 항목에 기본값을 따로 지정하면 이때도 결과를 확인하는 게 가능합니다.
map_chr(e, 'z', .default=NA)## [1] "a" "b" NA
네, 또 고생 많으셨습니다. 기지개라도 한 번 펴고 가겠습니다.
기본 연습 차원에서 이런 저런 자료 형태를 소개했지만 아마 R를 쓰시는 분들 대부분 가장 많이 다루는 자료형은 데이터 프레임일 겁니다. tidyverse에서는 tibble이라는 형태로 데이터 프레임을 다룹니다.
한번 예시 프레임을 만들어 보겠습니다.
f <- tibble(a=c(17, 23, 4, 10, 11),
b=c(24, 5, 6, 12, 18),
c=c(1, 7, 13, 19, 25),
d=c(8, 14, 20, 21, 2),
e=c(15, 16, 22, 3, 9))이때 각 행(가로) 합을 구하는 건 쉽습니다. 그저 mutate() 함수로 더해주면 그만입니다.
f %>%
mutate(sum=a+b+c+d+e)## # A tibble: 5 x 6 ## a b c d e sum ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 17 24 1 8 15 65 ## 2 23 5 7 14 16 65 ## 3 4 6 13 20 22 65 ## 4 10 12 19 21 3 65 ## 5 11 18 25 2 9 65
아, 혹시 mutate() 함수가 나와서 당황한 분이 계시다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
각 행 가운데 제일 큰 값을 뽑아내고 싶을 때는 어떻게 하면 될까요?
이때는 rowwise() 함수 도움을 받으면 됩니다. 아래처럼 말입니다.
f %>%
rowwise() %>%
mutate(max=max(a, b, c, d, e))## Source: local data frame [5 x 6] ## Groups: <by row> ## ## # A tibble: 5 x 6 ## a b c d e max ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 17 24 1 8 15 24 ## 2 23 5 7 14 16 23 ## 3 4 6 13 20 22 22 ## 4 10 12 19 21 3 21 ## 5 11 18 25 2 9 25
그러면 열(세로)은 어떻게 하느냐? 네, 이번에는 map() 함수 도움을 받아야 합니다. 일단 기본형.
map(f, sum)## $a ## [1] 65 ## ## $b ## [1] 65 ## ## $c ## [1] 65 ## ## $d ## [1] 65 ## ## $e ## [1] 65
원래 자료가 데이터 프레임(Data Frame) 형태니까 같은 형태로 받으면 더 좋겠죠? 그럴 때는 map_df() 함수를 쓰면 됩니다.
map_df(f, sum)## # A tibble: 1 x 5 ## a b c d e ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 65 65 65 65 65
물론 각 원소별로 함수를 적용하는 것도 가능합니다. 각 자리 숫자에 1씩 더해볼까요?
map_df(f, ~.x+1)## # A tibble: 5 x 5 ## a b c d e ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 18 25 2 9 16 ## 2 24 6 8 15 17 ## 3 5 7 14 21 23 ## 4 11 13 20 22 4 ## 5 12 19 26 3 10
map() 함수 사촌 가운데 modify() 함수가 있습니다. 이 함수는 원래 자료형 그대로 결과를 노출합니다.
우리는 데이터 프레임을 다루고 있으니까 modify() 함수를 쓰면 데이터 프레임으로 출력하빈다.
modify(f, ~.x+1)## # A tibble: 5 x 5 ## a b c d e ## <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 18 25 2 9 16 ## 2 24 6 8 15 17 ## 3 5 7 14 21 23 ## 4 11 13 20 22 4 ## 5 12 19 26 3 10
앞부분을 잘 보셨다면 별로 어려울 게 없는 내용입니다. 그러면 지금부터 처음에 던진 질문을 해결해 보겠습니다.
일단 kbo에서 우리에게 필요한 열만 뽑겠습니다.
kbo %>%
select(타율, 출루율, 장타력, OPS)## # A tibble: 303 x 4 ## 타율 출루율 장타력 OPS ## <dbl> <dbl> <dbl> <dbl> ## 1 0.282 0.354 0.41 0.764 ## 2 0.266 0.347 0.392 0.739 ## 3 0.283 0.349 0.412 0.762 ## 4 0.261 0.328 0.408 0.736 ## 5 0.256 0.344 0.373 0.717 ## 6 0.24 0.304 0.345 0.648 ## 7 0.263 0.333 0.393 0.726 ## 8 0.267 0.332 0.385 0.717 ## 9 0.256 0.323 0.354 0.677 ## 10 0.259 0.324 0.362 0.687 ## # ... with 293 more rows
그리고 그냥 lm() 함수를 써서 회귀분석을 진행합니다.
kbo %>%
select(타율, 출루율, 장타력, OPS) %>%
map(~lm(kbo$타석당득점~.x))## $타율 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Coefficients: ## (Intercept) .x ## -0.1259 0.9222 ## ## ## $출루율 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Coefficients: ## (Intercept) .x ## -0.1647 0.8328 ## ## ## $장타력 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Coefficients: ## (Intercept) .x ## -0.03826 0.39921 ## ## ## $OPS ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Coefficients: ## (Intercept) .x ## -0.09679 0.29370
오, 뭔가 됩니다. lm() 다음에 summary()도 써야 한다고 했으니 똑같이 씁니다.
kbo %>%
select(타율, 출루율, 장타력, OPS) %>%
map(~lm(kbo$타석당득점~.x)) %>%
map(summary)## $타율 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.0257544 -0.0053930 0.0005455 0.0042343 0.0216233 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.125931 0.007615 -16.54 <2e-16 *** ## .x 0.922209 0.028581 32.27 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.007795 on 301 degrees of freedom ## Multiple R-squared: 0.7757, Adjusted R-squared: 0.775 ## F-statistic: 1041 on 1 and 301 DF, p-value: < 2.2e-16 ## ## ## $출루율 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.0249185 -0.0052528 -0.0002377 0.0047366 0.0232427 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.164745 0.008966 -18.38 <2e-16 *** ## .x 0.832831 0.026250 31.73 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.007897 on 301 degrees of freedom ## Multiple R-squared: 0.7698, Adjusted R-squared: 0.769 ## F-statistic: 1007 on 1 and 301 DF, p-value: < 2.2e-16 ## ## ## $장타력 ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.0136000 -0.0043552 -0.0002445 0.0037635 0.0227674 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.038256 0.003431 -11.15 <2e-16 *** ## .x 0.399205 0.008649 46.16 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.005791 on 301 degrees of freedom ## Multiple R-squared: 0.8762, Adjusted R-squared: 0.8758 ## F-statistic: 2130 on 1 and 301 DF, p-value: < 2.2e-16 ## ## ## $OPS ## ## Call: ## lm(formula = kbo$타석당득점 ~ .x) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.010735 -0.003197 -0.000217 0.003359 0.020748 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.096792 0.003787 -25.56 <2e-16 *** ## .x 0.293700 0.005132 57.23 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.004775 on 301 degrees of freedom ## Multiple R-squared: 0.9158, Adjusted R-squared: 0.9156 ## F-statistic: 3275 on 1 and 301 DF, p-value: < 2.2e-16
이번에도 잘 나왔습니다. R²만 뽑아내고 싶을 때는 어떻게 하면 된다고 말씀드렸었죠?
kbo %>%
select(타율, 출루율, 장타력, OPS) %>%
map(~lm(kbo$타석당득점~.x)) %>%
map(summary) %>%
map('r.squared')## $타율 ## [1] 0.7757327 ## ## $출루율 ## [1] 0.7698036 ## ## $장타력 ## [1] 0.8762027 ## ## $OPS ## [1] 0.9158364
네, 잘 나왔습니다. 이제 출력 형태만 데이터 프레임으로 바꾸면 그만입니다. 맨 마지막 map() 함수만 map_df()로 바꿨습니다.
kbo %>%
select(타율, 출루율, 장타력, OPS) %>%
map(~lm(kbo$타석당득점~.x)) %>%
map(summary) %>%
map_df('r.squared')## # A tibble: 1 x 4 ## 타율 출루율 장타력 OPS ## <dbl> <dbl> <dbl> <dbl> ## 1 0.776 0.770 0.876 0.916
사실 이렇게 일일이 map()을 쓰지 않고 한번에 써도 같은 결과를 얻을 수 있습니다. 최종 결과물을 데이터 프레임으로 받아 보고 싶은 거니까 map_df()로 나머지를 전부 묶으면 됩니다.
kbo %>%
select(타율, 출루율, 장타력, OPS) %>%
map_df(~lm(kbo$타석당득점~.x) %>%
summary() %>%
.$r.squared)## # A tibble: 1 x 4 ## 타율 출루율 장타력 OPS ## <dbl> <dbl> <dbl> <dbl> ## 1 0.776 0.770 0.876 0.916
아, 맨 마지막 r.squared를 처리하는 방식이 하나 하나 map()을 쓸 때와 한꺼번에 처리할 때 차이가 납니다. 괜히 위에서 데이터를 골라내는 연습을 한 게 아닙니다.
문제가 그것만 있을까요?
아닙니다.
여기서 우리는 'kbo$타석당득점'이라는 자료와 나머지 자료를 비교했습니다. 지금은 이렇게 써도 결과를 얻는 데 문제가 없지만 기준이 달라지면 문제가 될 수 있습니다.
새로운 접근법이 필요한 이유입니다.
데이터 프레임을 쉬운 말로 하면 '표'라고 할 수 있습니다.
표는 크게 두 가지 방식으로 그릴 수 있는데 이를 흔히 와이드(wide), 롱(long) 폼으로 구분합니다. 아래는 같은 자료를 와이드 폼(왼쪽)과 롱 폼으로 그린 결과물입니다.
사람이 보기에는 확실히 와이드 폼이 낫습니다. 그런데 컴퓨터는 롱 폼 처리에 더 익숙합니다.
그래서 tidyverse 생태계에는 이 두 가지 폼 사이를 오갈 수 있도록 tidyr 패키지가 기다리고 있습니다.
이 작업이 낯설게 느껴지시는 분은 '최대한 친절하게 쓴 R로 데이터 깔끔하게 만들기(feat. tidyr)' 포스트가 도움이 될 수 있습니다.
우리가 쓰고 있는 kbo도 와이드 폼에 가깝습니다. 이 자료를 (일부) 롱 폼으로 바꿔보겠습니다.
아래 있는 코드는 타율부터 OPS에 해당하는 열을 기록과 값이라는 열로 모아 길게 바꾸라는 내용입니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값')## # A tibble: 1,212 x 6 ## 연도 X10년대 구단 타석당득점 기록 값 ## <int> <fct> <fct> <dbl> <chr> <dbl> ## 1 1982 a MBC 0.137 타율 0.282 ## 2 1982 a MBC 0.137 출루율 0.354 ## 3 1982 a MBC 0.137 장타력 0.41 ## 4 1982 a MBC 0.137 OPS 0.764 ## 5 1982 a 삼성 0.141 타율 0.266 ## 6 1982 a 삼성 0.141 출루율 0.347 ## 7 1982 a 삼성 0.141 장타력 0.392 ## 8 1982 a 삼성 0.141 OPS 0.739 ## 9 1982 a OB 0.129 타율 0.283 ## 10 1982 a OB 0.129 출루율 0.349 ## # ... with 1,202 more rows
tidyr에 익숙하신 분 가운데도 pivot_longer() 함수는 낯설다는 분이 적지 않을 겁니다. tidyr를 업데이트하면서 따라 새로 생긴 함수라 그렇습니다.
옛날 방식으로 gather() 함수를 사용한다면 아래처럼 코드를 쓸 수 있습니다.
kbo %>%
gather(key=기록, value=값, 타율:OPS)## # A tibble: 1,212 x 6 ## 연도 X10년대 구단 타석당득점 기록 값 ## <int> <fct> <fct> <dbl> <chr> <dbl> ## 1 1982 a MBC 0.137 타율 0.282 ## 2 1982 a 삼성 0.141 타율 0.266 ## 3 1982 a OB 0.129 타율 0.283 ## 4 1982 a 해태 0.125 타율 0.261 ## 5 1982 a 롯데 0.115 타율 0.256 ## 6 1982 a 삼미 0.102 타율 0.24 ## 7 1983 a 삼성 0.116 타율 0.263 ## 8 1983 a 해태 0.113 타율 0.267 ## 9 1983 a MBC 0.109 타율 0.256 ## 10 1983 a OB 0.111 타율 0.259 ## # ... with 1,202 more rows
죄송하지만 아직 머리 아프다고 하시기에는 이릅니다.
이제 이 데이터 프레임을 기록 그러니까 타율, 출루율, 장타력, OPS를 기준으로 그룹을 만들겠습니다. 그리고 이 nest() 함수를 써서 이 내용을 중첩하겠습니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(기록) %>%
nest()## # A tibble: 4 x 2 ## # Groups: 기록 [4] ## 기록 data ## <chr> <list<df[,5]>> ## 1 타율 [303 x 5] ## 2 출루율 [303 x 5] ## 3 장타력 [303 x 5] ## 4 OPS [303 x 5]
데이터가 모두 어디론가 숨어 버렸습니다. 작업 편의상 이 결과물을 kbo2라는 변수에 넣겠습니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(기록) %>%
nest() -> kbo2데이터 프레임을 중첩하면서 생긴 data 열에는 어떤 내용이 들어 있을까요? 열어 보면 됩니다.
kbo2$data## [[1]] ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.282 ## 2 1982 a 삼성 0.141 0.266 ## 3 1982 a OB 0.129 0.283 ## 4 1982 a 해태 0.125 0.261 ## 5 1982 a 롯데 0.115 0.256 ## 6 1982 a 삼미 0.102 0.24 ## 7 1983 a 삼성 0.116 0.263 ## 8 1983 a 해태 0.113 0.267 ## 9 1983 a MBC 0.109 0.256 ## 10 1983 a OB 0.111 0.259 ## # ... with 293 more rows ## ## [[2]] ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.354 ## 2 1982 a 삼성 0.141 0.347 ## 3 1982 a OB 0.129 0.349 ## 4 1982 a 해태 0.125 0.328 ## 5 1982 a 롯데 0.115 0.344 ## 6 1982 a 삼미 0.102 0.304 ## 7 1983 a 삼성 0.116 0.333 ## 8 1983 a 해태 0.113 0.332 ## 9 1983 a MBC 0.109 0.323 ## 10 1983 a OB 0.111 0.324 ## # ... with 293 more rows ## ## [[3]] ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.41 ## 2 1982 a 삼성 0.141 0.392 ## 3 1982 a OB 0.129 0.412 ## 4 1982 a 해태 0.125 0.408 ## 5 1982 a 롯데 0.115 0.373 ## 6 1982 a 삼미 0.102 0.345 ## 7 1983 a 삼성 0.116 0.393 ## 8 1983 a 해태 0.113 0.385 ## 9 1983 a MBC 0.109 0.354 ## 10 1983 a OB 0.111 0.362 ## # ... with 293 more rows ## ## [[4]] ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.764 ## 2 1982 a 삼성 0.141 0.739 ## 3 1982 a OB 0.129 0.762 ## 4 1982 a 해태 0.125 0.736 ## 5 1982 a 롯데 0.115 0.717 ## 6 1982 a 삼미 0.102 0.648 ## 7 1983 a 삼성 0.116 0.726 ## 8 1983 a 해태 0.113 0.717 ## 9 1983 a MBC 0.109 0.677 ## 10 1983 a OB 0.111 0.687 ## # ... with 293 more rows
우리가 데이터 프레임을 묶을 때 쓴 각 공격 기록을 기준으로 데이터를 나눈 걸 알 수 있습니다. [[1]] , [[2]] , [[3]] , [[4]] 보다는 기록 이름으로 구분하는 게 보기 편하겠죠?
set_names() 함수를 써서 붙여 주겠습니다. 우리가 쓰려는 기록 이름은 kbo2$기록에 들어 있습니다.
kbo2$data %>%
set_names(., kbo2$기록)## $타율 ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.282 ## 2 1982 a 삼성 0.141 0.266 ## 3 1982 a OB 0.129 0.283 ## 4 1982 a 해태 0.125 0.261 ## 5 1982 a 롯데 0.115 0.256 ## 6 1982 a 삼미 0.102 0.24 ## 7 1983 a 삼성 0.116 0.263 ## 8 1983 a 해태 0.113 0.267 ## 9 1983 a MBC 0.109 0.256 ## 10 1983 a OB 0.111 0.259 ## # ... with 293 more rows ## ## $출루율 ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.354 ## 2 1982 a 삼성 0.141 0.347 ## 3 1982 a OB 0.129 0.349 ## 4 1982 a 해태 0.125 0.328 ## 5 1982 a 롯데 0.115 0.344 ## 6 1982 a 삼미 0.102 0.304 ## 7 1983 a 삼성 0.116 0.333 ## 8 1983 a 해태 0.113 0.332 ## 9 1983 a MBC 0.109 0.323 ## 10 1983 a OB 0.111 0.324 ## # ... with 293 more rows ## ## $장타력 ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.41 ## 2 1982 a 삼성 0.141 0.392 ## 3 1982 a OB 0.129 0.412 ## 4 1982 a 해태 0.125 0.408 ## 5 1982 a 롯데 0.115 0.373 ## 6 1982 a 삼미 0.102 0.345 ## 7 1983 a 삼성 0.116 0.393 ## 8 1983 a 해태 0.113 0.385 ## 9 1983 a MBC 0.109 0.354 ## 10 1983 a OB 0.111 0.362 ## # ... with 293 more rows ## ## $OPS ## # A tibble: 303 x 5 ## 연도 X10년대 구단 타석당득점 값 ## <int> <chr> <chr> <dbl> <dbl> ## 1 1982 a MBC 0.137 0.764 ## 2 1982 a 삼성 0.141 0.739 ## 3 1982 a OB 0.129 0.762 ## 4 1982 a 해태 0.125 0.736 ## 5 1982 a 롯데 0.115 0.717 ## 6 1982 a 삼미 0.102 0.648 ## 7 1983 a 삼성 0.116 0.726 ## 8 1983 a 해태 0.113 0.717 ## 9 1983 a MBC 0.109 0.677 ## 10 1983 a OB 0.111 0.687 ## # ... with 293 more rows
이렇게 자료가 깔끔하게 들어 있기 때문에 위에서 했던 것처럼 map_df()를 써서 반복 작업을 처리하면 그만입니다. 사실 이름을 지정한 건 데이터 프레임에 이름이 나오면 더 보기 좋다는 게 제일 큰 이유였습니다.
kbo2$data %>%
set_names(., kbo2$기록) %>%
map_df(~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)## # A tibble: 1 x 4 ## 타율 출루율 장타력 OPS ## <dbl> <dbl> <dbl> <dbl> ## 1 0.776 0.770 0.876 0.916
생각대로 잘 나왔습니다.
그런데 처음에 각 기록 이름이 나왔는데 굳이 다른 변수에 자료를 저장했다가 이름까지 붙이는 과정을 거치는 게 최선이었을까요?
사실은 중첩 데이터 프레임 안이 어떻게 생겼는지 보여드리려고 일부러 열어본 겁니다. 곧바로 작업을 하는 것도 가능합니다.
작업을 마치면 계산한 열을 새로 추가해야 하니까 mutate() 함수를 쓰면 됩니다.
아, 데이터 프레임을 묶었으면 풀 수도 있겠죠? 그럴 때는 unnest() 함수를 씁니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model)## # A tibble: 4 x 3 ## # Groups: 기록 [4] ## 기록 data model ## <chr> <list<df[,5]>> <dbl> ## 1 타율 [303 x 5] 0.776 ## 2 출루율 [303 x 5] 0.770 ## 3 장타력 [303 x 5] 0.876 ## 4 OPS [303 x 5] 0.916
짜잔, 한 번에 잘 나왔습니다. 이 결과에서 data 열만 사라지면 좀 더 그럴 듯하겠죠? 특정 열을 제외할 때는 이름에 빼기(-) 부호를 붙이기만 하면 됩니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model) %>%
select(-data)## # A tibble: 4 x 2 ## # Groups: 기록 [4] ## 기록 model ## <chr> <dbl> ## 1 타율 0.776 ## 2 출루율 0.770 ## 3 장타력 0.876 ## 4 OPS 0.916
이렇게 한 번 코드를 만들고 나면 기준을 바꿔서 반복 작업을 처리할 수 있습니다.
이번에는 10년대 구분까지 넣어 묶어 보겠습니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(X10년대, 기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model) %>%
select(-data)## # A tibble: 16 x 3 ## # Groups: X10년대, 기록 [16] ## X10년대 기록 model ## <fct> <chr> <dbl> ## 1 a 타율 0.673 ## 2 a 출루율 0.649 ## 3 a 장타력 0.797 ## 4 a OPS 0.829 ## 5 b 타율 0.828 ## 6 b 출루율 0.831 ## 7 b 장타력 0.883 ## 8 b OPS 0.929 ## 9 c 타율 0.626 ## 10 c 출루율 0.680 ## 11 c 장타력 0.800 ## 12 c OPS 0.872 ## 13 d 타율 0.815 ## 14 d 출루율 0.722 ## 15 d 장타력 0.896 ## 16 d OPS 0.936
이 결과는 롱 폼입니다. 와이드 폼으로 바꾸려면 맨 아래 pivot_wider()를 추가하면 됩니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(X10년대, 기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model) %>%
select(-data) %>%
pivot_wider(names_from=기록, values_from=model)## # A tibble: 4 x 5 ## # Groups: X10년대 [4] ## X10년대 타율 출루율 장타력 OPS ## <fct> <dbl> <dbl> <dbl> <dbl> ## 1 a 0.673 0.649 0.797 0.829 ## 2 b 0.828 0.831 0.883 0.929 ## 3 c 0.626 0.680 0.800 0.872 ## 4 d 0.815 0.722 0.896 0.936
pivot_longer()와 마찬가지로 이번에도 낯선 분이 계실 텐데 옛날 식으로 spread() 함수를 쓰면 이렇게 나타낼 수 있습니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(X10년대, 기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model) %>%
select(-data) %>%
spread(key=기록, value=model)## # A tibble: 4 x 5 ## # Groups: X10년대 [4] ## X10년대 OPS 장타력 출루율 타율 ## <fct> <dbl> <dbl> <dbl> <dbl> ## 1 a 0.829 0.797 0.649 0.673 ## 2 b 0.929 0.883 0.831 0.828 ## 3 c 0.872 0.800 0.680 0.626 ## 4 d 0.936 0.896 0.722 0.815
물론 이렇게 얻은 결과를 가지고 곧바로 그래프도 그릴 수 있습니다.
컴퓨터는 롱 폼을 좋아하니까 pivot_wider()는 빼고 가겠습니다.
아래는 각 10년대마다 각 타격 기록이 타석당 득점 기록을 얼마나 설명했는지 보여주는 그래프를 그리는 코드입니다.
kbo %>%
pivot_longer(cols=타율:OPS, names_to='기록', values_to='값') %>%
group_by(X10년대, 기록) %>%
nest() %>%
mutate(model=map(data, ~lm(타석당득점~값, data=.) %>%
summary() %>%
.$r.squared)) %>%
unnest(model) %>%
select(-data) %>%
ggplot(aes(X10년대, model, fill=기록)) +
geom_bar(stat='identity', position=position_dodge2(reverse=TRUE)) +
scale_fill_viridis_d()정말 고생 많으셨습니다. 여기까지 공부하셨으면 이제 컴퓨터를 닫고 (설마 지금 휴대전화로 보고 계신 건 아니시죠?) 좀 쉬셔도 됩니다.
아래 내용은 '방과후 교실' 같은 개념이니까 나중에 다시 오셔도 관계없습니다.
어쩌면 마음만 먹는다면 위에 나오는 12번 정도 반복 작업을 하는 건 그냥 일일이 수작업을 하는 게 더 쉬울지도 모릅니다.
400번은 어떨까요?
세계 보건 과학자 네트워크(NCD·RiSc) 홈페이지에 올라와 있는 국가별 평균 키 데이터를 불러서 한번 알아보겠습니다.
이건 영문 자료라 인코딩 문제가 없으니 read_csv() 함수를 쓰겠습니다.
height <- read_csv('http://ncdrisc.org/downloads/height/NCD_RisC_eLife_2016_height_age18_countries.csv')## Parsed with column specification: ## cols( ## Country = col_character(), ## ISO = col_character(), ## Sex = col_character(), ## `Year of birth` = col_double(), ## `Mean height (cm)` = col_double(), ## `Mean height lower 95% uncertainty interval (cm)` = col_double(), ## `Mean height upper 95% uncertainty interval (cm)` = col_double() ## )
이 자료가 어떻게 생겼는지도 한번 map() 계열 함수를 써서 확인해 볼까요?
height %>%
map_df(~(tibble(class=class(.x),
count=n_distinct(.x))),
.id='variable')## # A tibble: 7 x 3 ## variable class count ## <chr> <chr> <int> ## 1 Country character 200 ## 2 ISO character 200 ## 3 Sex character 2 ## 4 Year of birth numeric 101 ## 5 Mean height (cm) numeric 40393 ## 6 Mean height lower 95% uncertainty interval (cm) numeric 40398 ## 7 Mean height upper 95% uncertainty interval (cm) numeric 40397
그러니까 이 파일에는 101년 동안 태어난 200개 나라 남녀별 평균 키가 들어 있습니다. 태어난 연도와 평균 키 사이에 어떤 관계가 있는지 알아보려면 코드를 어떻게 짜면 될까요?
위에서 썼던 코드를 살짝만 수정하면 됩니다.
height %>%
group_by(Country, Sex) %>%
nest() %>%
mutate(model=map_dbl(data, ~lm(`Mean height (cm)`~`Year of birth`, data=.) %>%
summary() %>%
.$r.squared))## # A tibble: 400 x 4 ## # Groups: Country, Sex [400] ## Country Sex data model ## <chr> <chr> <list<df[,5]>> <dbl> ## 1 Afghanistan Men [101 x 5] 0.961 ## 2 Albania Men [101 x 5] 0.943 ## 3 Algeria Men [101 x 5] 0.792 ## 4 American Samoa Men [101 x 5] 0.958 ## 5 Andorra Men [101 x 5] 0.988 ## 6 Angola Men [101 x 5] 0.795 ## 7 Antigua and Barbuda Men [101 x 5] 0.930 ## 8 Argentina Men [101 x 5] 0.975 ## 9 Armenia Men [101 x 5] 0.925 ## 10 Australia Men [101 x 5] 0.978 ## # ... with 390 more rows
열이 400개니까 회귀분석을 400번 한 겁니다. system.time() 함수로 알아보면 이렇게 많은 작업을 하는 데 채 0.1초도 걸리지 않습니다.
height %>%
group_by(Country, Sex) %>%
nest() %>%
mutate(model=map_dbl(data, ~lm(`Mean height (cm)`~`Year of birth`, data=.) %>%
summary() %>%
.$r.squared)) %>%
system.time()## user system elapsed ## 0 0 0
사실 프로야구 자료를 다룰 때는 함수 쪽은 아예 건드리지 않았습니다. 그 부분도 이렇게 정리해서 쓸 수 있습니다. 먼저 함수 정의.
height_model <- function(x){
lm(`Mean height (cm)`~`Year of birth`, data=x) %>%
summary() %>%
.$r.squared
}
이렇게 정리를 하고 나면 add_one()을 썼을 때처럼 뒤에 함수 이름만 쓰면 됩니다.
height %>%
group_by(Country, Sex) %>%
nest() %>%
mutate(model=map_dbl(data, height_model))## # A tibble: 400 x 4 ## # Groups: Country, Sex [400] ## Country Sex data model ## <chr> <chr> <list<df[,5]>> <dbl> ## 1 Afghanistan Men [101 x 5] 0.961 ## 2 Albania Men [101 x 5] 0.943 ## 3 Algeria Men [101 x 5] 0.792 ## 4 American Samoa Men [101 x 5] 0.958 ## 5 Andorra Men [101 x 5] 0.988 ## 6 Angola Men [101 x 5] 0.795 ## 7 Antigua and Barbuda Men [101 x 5] 0.930 ## 8 Argentina Men [101 x 5] 0.975 ## 9 Armenia Men [101 x 5] 0.925 ## 10 Australia Men [101 x 5] 0.978 ## # ... with 390 more rows
여기까지 오고 나니 어떤 나라 어떤 성별이 R²가 제일 높을지 궁금하지 않으신가요? arrange() 함수를 쓰시면 바로 정렬할 수 있습니다.
height %>%
group_by(Country, Sex) %>%
nest() %>%
mutate(model=map_dbl(data, height_model)) %>%
arrange(-model)## # A tibble: 400 x 4 ## # Groups: Country, Sex [400] ## Country Sex data model ## <chr> <chr> <list<df[,5]>> <dbl> ## 1 Brazil Men [101 x 5] 0.999 ## 2 Taiwan Women [101 x 5] 0.999 ## 3 North Korea Men [101 x 5] 0.999 ## 4 Thailand Men [101 x 5] 0.996 ## 5 Italy Women [101 x 5] 0.996 ## 6 France Men [101 x 5] 0.996 ## 7 Tunisia Men [101 x 5] 0.995 ## 8 South Korea Men [101 x 5] 0.995 ## 9 Canada Men [101 x 5] 0.994 ## 10 Lebanon Women [101 x 5] 0.993 ## # ... with 390 more rows
한국과 북한 남성 모두 10위 안에 이름을 이름을 올리고 있습니다.
두 나라 남성들 평균 키가 어떻게 변했는지 그래프로 그려보겠습니다.
height %>%
filter(Country %in% c('North Korea', 'South Korea') & Sex=='Men') %>%
ggplot(aes(`Year of birth`, `Mean height (cm)`, color=Country)) +
geom_line()여기서 R²가 높다는 건 직선 모형에 가깝게 그러니까 예상 가능하게 평균 신장이 늘었다는 뜻일 뿐 키가 더 많이 커졌다는 뜻은 아닙니다.
이 직선 모형이 어떤 값을 예상했는지 한번 점을 찍어 볼까요?
이때는 modelr 패키지에서 제공하는 add_predictions() 함수를 쓰면 됩니다.
moderl 패키지도 tidyverse 생태계에 속하기는 하지만 자동으로 설치하거나 불러오지 않으니 설치하고 불러오는 과정을 거치셔야 합니다.
#install.packages('modelr')
library('modelr')함수 쓰기 연습을 한번 더 해볼까요?
이번에 만들 함수는 나라 이름을 받으면 그 나라 남성 평균 키는 실제로 어떻게 변했는지 모형이 예상한 키는 얼마였다고 알려주는 녀석입니다.
그러려면 먼저 height에서 입력 받은 나라 이름(x)과 이름이 같은 나라(Country) 남성 자료를 불러와야 합니다.
그리고 이 자료를 데이터 프레임(df)에 넣은 다음 lm() 함수를 써서 선형 회귀분석 모형을 만듭니다.
이어서 add_predictions() 함수로 이 모형이 예상한 값을 데이터 프레임에 붙여줍니다.
그다음 나라 이름, 출생연도(Year of birth), 평균 키(Mean height (cm))를 출력하도록 만듭니다.
이런 내용을 코드로 정리하면 아래처럼 쓸 수 있습니다.
height_model2 <- function(x){
df <- height %>%
filter(Country==x & Sex=='Men')
fit <- lm(`Mean height (cm)`~`Year of birth`, data=df)
df %>%
add_predictions(fit) %>%
select(Country, `Year of birth`, `Mean height (cm)`, pred)
}이 height_model2가 잘 작동하는지는 함수 적용 결과를 확인하면 알 수 있습니다.
map_df(c('South Korea', 'North Korea'), height_model2) %>%
map_chr(n_distinct)## Country Year of birth Mean height (cm) pred ## "2" "101" "202" "202"
나라는 두 개, 출생연도는 101개 그래서 자료는 202개가 잘 들어왔습니다.
그림을 그려봐야겠죠? 위에서 본 것처럼 실제 평균 키는 직선으로 그리고, 모형 예상값은 점으로 찍겠습니다.
map_df(c('South Korea', 'North Korea'), height_model2) %>%
ggplot(aes(`Year of birth`, `Mean height (cm)`, color=Country)) +
geom_line(lwd=1.5) +
geom_point(aes(y=pred), alpha=.5)예, 북한은 정말 직선에 가깝게 움직였습니다. 한국은 최근 들어 성장세가 한풀 꺾인 분위기입니다.
이렇게 함수를 한 번 만들고 나면 당연히 다른 나라에도 쓸 수 있습니다. 나라 숫자도 상관이 없습니다. 지금 그래프에 일본을 더해 볼까요?
map_df(c('South Korea', 'North Korea', 'Japan'), height_model2) %>%
ggplot(aes(`Year of birth`, `Mean height (cm)`, color=Country)) +
geom_line(lwd=1.5) +
geom_point(aes(y=pred), alpha=.5)예, 그래프를 보시면 아실 수 있는 것처럼 일본은 남자들 키가 가파르게 커지다가 아주 많이 꺾였습니다.
이렇게 함수를 만들면 확실히 반복 작업을 별로 어렵지 않게 처리할 수 있습니다.
한번 더 말씀드리지만 R²가 높다는 건 예측 가능한 방향으로 평균 신장이 늘었다는 뜻일 뿐 평균 키가 제일 많이 늘었다는 뜻은 아닙니다.
그러면 그건 어떤 나라 어떤 성별일까요?
dplyr 패키지에 들어 있는 주요 동사를 조합하면 어렵지 않게 뽑아낼 수 있습니다.
height %>%
group_by(Country, Sex) %>%
summarise(height_min=min(`Mean height (cm)`),
height_max=max(`Mean height (cm)`),
diff=height_max-height_min) %>%
arrange(-diff)## # A tibble: 400 x 5 ## # Groups: Country [200] ## Country Sex height_min height_max diff ## <chr> <chr> <dbl> <dbl> <dbl> ## 1 South Korea Women 142. 162. 20.2 ## 2 Iran Men 157. 174. 16.5 ## 3 Japan Women 142. 159. 16.2 ## 4 Croatia Women 150. 166. 15.8 ## 5 Czech Republic Women 153. 168. 15.7 ## 6 Serbia Women 152. 168. 15.7 ## 7 Greenland Men 159. 175. 15.4 ## 8 Japan Men 156. 171. 15.3 ## 9 South Korea Men 160. 175. 15.2 ## 10 Greece Men 163. 177. 14.9 ## # ... with 390 more rows
네, 한국 여성이 지난 101년 동안 평균 신장이 제일 크게 늘었습니다. 1896년생은 평균 키가 142㎝밖에 되지 않았는데 1996년생은 162㎝로 20㎝ 늘었습니다.
물론 위에 있는 코드는 갈수록 평균 키가 커진다는 정보를 바탕으로 평균 키가 가장 작았을 때와 클 때를 비교한 겁니다.
실제로 1896년생과 1996년생 평균 키를 비교하는 코드는 아래와 같습니다.
height %>%
group_by(Country, Sex) %>%
filter(`Year of birth` %in% c(1896, 1996)) %>%
select(Country, Sex, `Year of birth`, `Mean height (cm)`) %>%
pivot_wider(names_from=`Year of birth`, values_from=`Mean height (cm)`) %>%
mutate(diff=`1996`-`1896`) %>%
arrange(-diff)## # A tibble: 400 x 5 ## # Groups: Country, Sex [400] ## Country Sex `1896` `1996` diff ## <chr> <chr> <dbl> <dbl> <dbl> ## 1 South Korea Women 142. 162. 20.2 ## 2 Iran Men 157. 174. 16.5 ## 3 Japan Women 142. 158. 16.0 ## 4 Czech Republic Women 153. 168. 15.7 ## 5 Serbia Women 152. 168. 15.7 ## 6 Croatia Women 150. 166. 15.6 ## 7 Greenland Men 159. 175. 15.4 ## 8 South Korea Men 160. 175. 15.2 ## 9 Greece Men 163. 177. 14.8 ## 10 Japan Men 156. 171. 14.7 ## # ... with 390 more rows
나머지 순위는 소소하게 바뀌었지만 한국 여성이 평균 키가 제일 늘었다는 사실은 달라지지 않았습니다. 한국 남성도 여덟 번째로 평균 키가 늘었습니다. 그러니 여러분 키 자랑하셔도 됩니다.
올해는 2020년 그러니까 20을 두 번 반복하는 해입니다. 게다가 이 글 번호도 1199. 그래서 어쩐지 새해 첫 글은 반복문이 좋지 않을까 하는 생각에 정리해 봤습니다.
언제나 그렇듯 잘못된 또는 이해가 잘 가지 않는 내용이 있거나 같은 작업을 처리하는 더 좋은 방법을 알고 계시면 댓글로 남겨주시면 대단히 감사하겠습니다.
그러면 여러분 모두 Happy Tidyversing!
댓글,