제 7강: 반복

데이터과학 입문

Author

Affiliation

원중호

서울대학교 통계학과

시작하기 전에

다음의 패키지가 설치되어 있지 않으면 설치한다.

# install.packages("tidyverse")
# install.packages("mdsr")
# install.packages("Lahman")
# install.packages("bench")
# install.packages("NHANES")
# install.packages("patchwork")
library(tidyverse)
library(mdsr)
library(Lahman)
library(bench)
library(NHANES)
library(patchwork)
sessionInfo()

R version 4.3.3 (2024-02-29)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Asia/Seoul
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] patchwork_1.2.0 NHANES_2.1.0    bench_1.1.3     Lahman_11.0-0  
 [5] mdsr_0.2.7      lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1  
 [9] dplyr_1.1.4     purrr_1.0.2     readr_2.1.5     tidyr_1.3.1    
[13] tibble_3.2.1    ggplot2_3.5.0   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.4      jsonlite_1.8.8    compiler_4.3.3    tidyselect_1.2.0 
 [5] scales_1.3.0      yaml_2.3.8        fastmap_1.1.1     R6_2.5.1         
 [9] generics_0.1.3    skimr_2.1.5       knitr_1.45        htmlwidgets_1.6.4
[13] munsell_0.5.0     pillar_1.9.0      tzdb_0.4.0        rlang_1.1.3      
[17] utf8_1.2.4        stringi_1.8.3     repr_1.1.6        xfun_0.42        
[21] timechange_0.3.0  cli_3.6.2         withr_3.0.0       magrittr_2.0.3   
[25] digest_0.6.34     grid_4.3.3        rstudioapi_0.15.0 base64enc_0.1-3  
[29] hms_1.1.3         lifecycle_1.0.4   vctrs_0.6.5       evaluate_0.23    
[33] glue_1.7.0        fansi_1.0.6       colorspace_2.1-0  rmarkdown_2.25   
[37] tools_4.3.3       pkgconfig_2.0.3   htmltools_0.5.7  

반복

• 계산기: 사람이 직접 손으로 산수 계산을 할 필요가 없도록 해 줌

• 프로그래밍 언어: 사람이 직접 코드 덩어리를 다시 실행하거나 한두 부분만 변경해서 매번 복사하여 붙여넣는 반복 계산을 수행하지 않도록 해 줌

• 어떤 야구팀의 시즌별 성과를 보기위한 코드를 작성하면, 다른 팀에게도 똑같이 적용할 수 있어야 함

• 이러한 작업을 수행하기 위해 코드를 다시 입력할 필요가 없어야 함

벡터화 연산

루프

names(Lahman::Teams)

 [1] "yearID"         "lgID"           "teamID"         "franchID"      
 [5] "divID"          "Rank"           "G"              "Ghome"         
 [9] "W"              "L"              "DivWin"         "WCWin"         
[13] "LgWin"          "WSWin"          "R"              "AB"            
[17] "H"              "X2B"            "X3B"            "HR"            
[21] "BB"             "SO"             "SB"             "CS"            
[25] "HBP"            "SF"             "RA"             "ER"            
[29] "ERA"            "CG"             "SHO"            "SV"            
[33] "IPouts"         "HA"             "HRA"            "BBA"           
[37] "SOA"            "E"              "DP"             "FP"            
[41] "name"           "park"           "attendance"     "BPF"           
[45] "PPF"            "teamIDBR"       "teamIDlahman45" "teamIDretro"   

glimpse(Teams)

• Teams 데이터프레임의 15번부터 40번 열까지가 각 팀이 그 시즌에 어떻게 성과를 냈는지에 대한 숫자 데이터를 포함

• 15번부터 40번까지의 26개 수치형 변수의 평균 계산

• 방법 1: 모든 열 이름을 일일이 입력하여 mean() 명령어를 26번 반복

• 방법 2: for() 루프

averages <- NULL
for (i in 15:40) {
  averages[i - 14] <- mean(Teams[, i], na.rm = TRUE)
}
names(averages) <- names(Teams)[15:40]
averages

          R          AB           H         X2B         X3B          HR 
 681.158872 5132.173466 1339.252405  229.027861   45.431841  106.542289 
         BB          SO          SB          CS         HBP          SF 
 473.614594  768.118373  109.094810   46.260531   46.182553   44.038670 
         RA          ER         ERA          CG         SHO          SV 
 681.157877  573.989055    3.841881   47.084245    9.634163   24.587396 
     IPouts          HA         HRA         BBA         SOA           E 
4016.164842 1339.034163  106.542289  473.959867  767.538972  179.839469 
         DP          FP 
 132.602653    0.966538 

• 불편한 점: 원하는 인덱스를 찾아 써야 함 (14, 15, 40; “매직넘버”)

R과 벡터

• R 프로그래머는 벡터의 각 요소에 연산을 적용하여 이러한 유형의 문제를 해결하는 것을 선호
  • 인덱스에 호소하지 않고 코드 한 줄로 해결!
• 왜 그런가?
  • R의 기본 자료형은 벡터
  • 문자열도 길이가 1인 벡터
  • 원자 자료형(atomic data type)이 없음

a <- "a string"
class(a)

[1] "character"

is.vector(a)

[1] TRUE

length(a)

[1] 1

---

>>> a = "a string"  # Python
>>> len(a)
8

• 그 결과, R은 벡터화된 연산에 고도로 최적화
• 루프는 본질적으로 이 최적화의 이점을 활용하지 못함

Note

for() 반복문은 최대한 피하라.

벡터화 함수

• 벡터를 입력으로 받아 같은 길이의 벡터로 출력

exp(1:3)

[1]  2.718282  7.389056 20.085537

mean(1:3)

[1] 2

벡터화하면 큰일나는 함수

if (c(TRUE, FALSE)) {
  cat("This is a great book!")
}

Error in if (c(TRUE, FALSE)) {: the condition has length > 1

벡터화 함수와 반복 함수

x <- 1:10000
bench::mark(
  exp(x),
  map_dbl(x, exp)   # we will study this later
)

# A tibble: 2 × 6
  expression           min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 exp(x)           27.87µs   99.3µs     9404.   156.3KB     29.0
2 map_dbl(x, exp)   5.76ms   6.66ms      147.    85.1KB     16.1

dplyr::across()

하드코딩 피하기

• 앞 예제에서, Teams 데이터프레임의 15번부터 40번 열이 수치형임을 확인하고 이를 매직넘버로 하드코딩
• tidyselect::where() 함수를 사용하여 수치형(is.numeric() == TRUE)인 열을 자동으로 식별
• dplyr::across() 함수를 사용, 지정된 열 모두에 mean() 함수 적용

Teams %>%
  summarize(across(where(is.numeric), mean, na.rm = TRUE))

    yearID     Rank        G    Ghome        W        L        R       AB
1 1959.492 4.028192 150.1254 78.08257 74.67463 74.67463 681.1589 5132.173
         H      X2B      X3B       HR       BB       SO       SB       CS
1 1339.252 229.0279 45.43184 106.5423 473.6146 768.1184 109.0948 46.26053
       HBP       SF       RA       ER      ERA       CG      SHO      SV
1 46.18255 44.03867 681.1579 573.9891 3.841881 47.08425 9.634163 24.5874
    IPouts       HA      HRA      BBA     SOA        E       DP       FP
1 4016.165 1339.034 106.5423 473.9599 767.539 179.8395 132.6027 0.966538
  attendance      BPF      PPF
1    1385100 100.1877 100.2109

• 결과에는 이전에 정의한 범위를 벗어나는 변수(yearID, attendance 등)도 포함됨에 유의

---

• 변수 이름으로도 지정 가능

Teams %>%
  summarize(across(c(yearID, R:SF, BPF), mean, na.rm = TRUE))

    yearID        R       AB        H      X2B      X3B       HR       BB
1 1959.492 681.1589 5132.173 1339.252 229.0279 45.43184 106.5423 473.6146
        SO       SB       CS      HBP       SF      BPF
1 768.1184 109.0948 46.26053 46.18255 44.03867 100.1877

purrr::map()

map() 및 그 변종

• across()보다 일반적으로 목록(list), 벡터, 또는 데이터프레임의 열의 각 항목에 함수를 적용

Teams %>% 
  select(15:40) %>%
  map_dbl(mean, na.rm = TRUE)

          R          AB           H         X2B         X3B          HR 
 681.158872 5132.173466 1339.252405  229.027861   45.431841  106.542289 
         BB          SO          SB          CS         HBP          SF 
 473.614594  768.118373  109.094810   46.260531   46.182553   44.038670 
         RA          ER         ERA          CG         SHO          SV 
 681.157877  573.989055    3.841881   47.084245    9.634163   24.587396 
     IPouts          HA         HRA         BBA         SOA           E 
4016.164842 1339.034163  106.542289  473.959867  767.538972  179.839469 
         DP          FP 
 132.602653    0.966538 

• 인수 1: 작업을 수행하고자 하는 대상(이 경우 데이터프레임)
• 인수 2: 함수의 이름을 지정 (mean())
• 인수 3: 인수 2의 함수에 들어갈 인수 (na.rm)

---

Teams %>% 
  select(teamID) %>% 
  map_dbl(mean, na.rm = TRUE)

teamID 
    NA 

• map_dbl(): double형 벡터를 반환
• map(): list를 반환

1차원 벡터에서의 반복

기존 함수를 반복 적용

• 현재 Los Angeles Angels of Anaheim으로 알려진 메이저리그 야구팀의 팀명 변천사

angels <- Teams %>%
  filter(franchID == "ANA") %>%
  group_by(teamID, name) %>%
  summarize(began = first(yearID), ended = last(yearID)) %>%
  arrange(began)
angels

# A tibble: 4 × 4
# Groups:   teamID [3]
  teamID name                          began ended
  <fct>  <chr>                         <int> <int>
1 LAA    Los Angeles Angels             1961  1964
2 CAL    California Angels              1965  1996
3 ANA    Anaheim Angels                 1997  2004
4 LAA    Los Angeles Angels of Anaheim  2005  2022

---

• 이제 각 이름의 길이(글자 수)를 구해 보자.

angels_names <- angels %>% pull(name)
nchar(angels_names[1])

[1] 18

nchar(angels_names[2])

[1] 17

nchar(angels_names[3])

[1] 14

nchar(angels_names[4])

[1] 29

• 이보다는 벡터 angels_names의 각 원소에 nchar()함수를 적용하는 것이 더 간단하고 효율적이며, 확장 가능

map_int(angels_names, nchar)   # returns an integer vector

[1] 18 17 14 29

---

• 실은 nchar()는 이미 벡터화 되어 있다.

nchar(angels_names)

[1] 18 17 14 29

사용자 정의 함수를 반복 적용

• 승수 기준 상위 5개 시즌 표시 함수

top5 <- function(data, team_name) {
  data %>%
    filter(name == team_name) %>%
    select(teamID, yearID, W, L, name) %>%
    arrange(desc(W)) %>%
    head(n = 5)
}

---

• Angels의 상위 5개 시즌

angels_names %>%
  map(top5, data = Teams)   # list

[[1]]
  teamID yearID  W  L               name
1    LAA   1962 86 76 Los Angeles Angels
2    LAA   1964 82 80 Los Angeles Angels
3    LAA   1961 70 91 Los Angeles Angels
4    LAA   1963 70 91 Los Angeles Angels

[[2]]
  teamID yearID  W  L              name
1    CAL   1982 93 69 California Angels
2    CAL   1986 92 70 California Angels
3    CAL   1989 91 71 California Angels
4    CAL   1985 90 72 California Angels
5    CAL   1979 88 74 California Angels

[[3]]
  teamID yearID  W  L           name
1    ANA   2002 99 63 Anaheim Angels
2    ANA   2004 92 70 Anaheim Angels
3    ANA   1998 85 77 Anaheim Angels
4    ANA   1997 84 78 Anaheim Angels
5    ANA   2000 82 80 Anaheim Angels

[[4]]
  teamID yearID   W  L                          name
1    LAA   2008 100 62 Los Angeles Angels of Anaheim
2    LAA   2014  98 64 Los Angeles Angels of Anaheim
3    LAA   2009  97 65 Los Angeles Angels of Anaheim
4    LAA   2005  95 67 Los Angeles Angels of Anaheim
5    LAA   2007  94 68 Los Angeles Angels of Anaheim

---

angels_names %>% 
  map_dfr(top5, data = Teams) # data frame

   teamID yearID   W  L                          name
1     LAA   1962  86 76            Los Angeles Angels
2     LAA   1964  82 80            Los Angeles Angels
3     LAA   1961  70 91            Los Angeles Angels
4     LAA   1963  70 91            Los Angeles Angels
5     CAL   1982  93 69             California Angels
6     CAL   1986  92 70             California Angels
7     CAL   1989  91 71             California Angels
8     CAL   1985  90 72             California Angels
9     CAL   1979  88 74             California Angels
10    ANA   2002  99 63                Anaheim Angels
11    ANA   2004  92 70                Anaheim Angels
12    ANA   1998  85 77                Anaheim Angels
13    ANA   1997  84 78                Anaheim Angels
14    ANA   2000  82 80                Anaheim Angels
15    LAA   2008 100 62 Los Angeles Angels of Anaheim
16    LAA   2014  98 64 Los Angeles Angels of Anaheim
17    LAA   2009  97 65 Los Angeles Angels of Anaheim
18    LAA   2005  95 67 Los Angeles Angels of Anaheim
19    LAA   2007  94 68 Los Angeles Angels of Anaheim

---

• 이름별 상위 5개 시즌의 평균 승수

angels_names %>% 
  map_dfr(top5, data = Teams) %>%
  group_by(teamID, name) %>%
  summarize(N = n(), mean_wins = mean(W)) %>%
  arrange(desc(mean_wins))

# A tibble: 4 × 4
# Groups:   teamID [3]
  teamID name                              N mean_wins
  <fct>  <chr>                         <int>     <dbl>
1 LAA    Los Angeles Angels of Anaheim     5      96.8
2 CAL    California Angels                 5      90.8
3 ANA    Anaheim Angels                    5      88.4
4 LAA    Los Angeles Angels                4      77  

• 최상위 5개 시즌 동안의 평균 성적을 비교할 때 “Los Angeles Angels of Anaheim” 시절이 가장 성공적인 것으로 보임

하위 그룹에 대한 반복

purrr::group_modify()

• 먼저 데이터프레임의 그룹을 정의한 다음 각 그룹에 함수를 적용
• map_dfr()와 유사.
  • map_dfr(): 벡터의 개별 원소를 사용
  • group_modify(): 데이터프레임에서 정의된 그룹을 사용

---

예제: 기대 승률

• Bill James의 공식(4.2절) \[\widehat{WPct} = \frac{RS^2}{RS^2 + RA^2} = \frac{1}{1 + (RA/RS)^2}\]
  • \(RS\): 팀 득점, \(RA\): 팀 실점

---

exp_wpct <- function(x) { 
  return(1/(1 + (1/x)^2))
}

TeamRuns <- Teams %>% 
  filter(yearID >= 1954) %>%
  rename(RS = R) %>% 
  mutate(WPct = W / (W + L), run_ratio = RS/RA) %>%
  select(yearID, teamID, lgID, WPct, run_ratio)

plt <- ggplot(data = TeamRuns, aes(x = run_ratio, y = WPct)) +
  geom_vline(xintercept = 1, color = "darkgray", linetype = 2) +
  geom_hline(yintercept = 0.5, color = "darkgray", linetype = 2) +
  geom_point(alpha = 0.2) + 
  stat_function(fun = exp_wpct, size = 2, color = "blue") + 
  xlab("Ratio of Runs Scored to Runs Allowed") + 
  ylab("Winning Percentage")

---

---

최적 지수?

• 실제 데이터와 모형함수 \(f(x) = 1/(1+(1/x)^2)\)의 그래프는 굉장히 잘 맞는 것을 알 수 있다.
• 그러나 이 모형의 지수 2는 Bill James가 제안한 값이고, 실제 경기 데이터에는 지수가 1.85에 가까운 모형이 더 잘 적합한다는 사실이 알려져 있다.

TeamRuns %>%
  nls(
    formula = WPct ~ 1/(1 + (1/run_ratio)^k), 
    start = list(k = 2)
  ) %>%
  coef()

       k 
1.831562 

---

• 또한 시대에 따라 모형의 최적 지수가 달라짐이 알려져 있다.
  
• 시대별 최적 지수를 구하기 위해 group_modify() 함수를 사용한다.

fit_k <- function(x) {
  mod <- nls(
    formula = WPct ~ 1/(1 + (1/run_ratio)^k), 
    data = x,
    start = list(k = 2)
  )
  return(tibble(k = coef(mod), n = nrow(x)))
}

---

TeamRuns %>% 
  mutate(decade = yearID %/% 10 * 10) %>%
  group_by(decade) %>% 
  group_modify(~fit_k(.))

# A tibble: 8 × 3
# Groups:   decade [8]
  decade     k     n
   <dbl> <dbl> <int>
1   1950  1.69    96
2   1960  1.90   198
3   1970  1.74   246
4   1980  1.93   260
5   1990  1.88   278
6   2000  1.94   300
7   2010  1.77   300
8   2020  1.77    90

• 지수가 매 10년마다 동일하지는 않지만 1.69에서 1.94 사이의 범위 안에 있음

---

시즌별 홈런 1위 팀

hr_leader <- function(x) {
# x is a subset of Teams for a single year and league
  x %>% 
    select(teamID, HR) %>% 
    arrange(desc(HR)) %>% 
    head(1)
}
Teams %>% 
  filter(yearID == 1961 & lgID == "AL") %>% 
  hr_leader()

  teamID  HR
1    NYA 240

• 1961년에는 뉴욕 양키스가 아메리칸 리그에서 홈런을 가장 많이 친 것을 알 수 있다.

---

• group_modify()를 사용하면 홈런 수에서 리그 1위를 차지한 모든 팀을 찾을 수 있다.
  • .keep 인수를 사용하여 그룹화 변수가 계산에 나타나도록 함

hr_leaders <- Teams %>% 
  group_by(yearID, lgID) %>% 
  group_modify(~hr_leader(.x), .keep = TRUE)

tail(hr_leaders, 4)

# A tibble: 4 × 4
# Groups:   yearID, lgID [4]
  yearID lgID  teamID    HR
   <int> <fct> <fct>  <int>
1   2021 AL    TOR      262
2   2021 NL    SFN      241
3   2022 AL    NYA      254
4   2022 NL    ATL      243

---

• 가장 많은 홈런을 친 팀의 리그별 평균

hr_leaders %>%
  group_by(lgID) %>%
  summarize(mean_hr = mean(HR))

# A tibble: 7 × 2
  lgID  mean_hr
  <fct>   <dbl>
1 AA       40.5
2 AL      158. 
3 FL       51  
4 NA       13.8
5 NL      131. 
6 PL       66  
7 UA       32  

---

• 1916년 전에는 AL과 NL 이외의 리그도 있었음.

hr_leaders %>%
  filter(yearID >= 1916) %>%
  group_by(lgID) %>%
  summarize(mean_hr = mean(HR))

# A tibble: 2 × 2
  lgID  mean_hr
  <fct>   <dbl>
1 AL       176.
2 NL       162.

---

• 시간에 따른 시즌별 최다 팀 홈런수의 변화

hr_leaders %>% 
  filter(yearID >= 1916) %>%
  ggplot(aes(x = yearID, y = HR, color = lgID)) + 
  geom_line() + 
  geom_point() + 
  geom_smooth(se = FALSE) + 
  geom_vline(xintercept = 1973) + 
  annotate(
    "text", x = 1974, y = 25, 
    label = "AL adopts DH", hjust = "left"
  ) +
  geom_vline(xintercept = 2022) + 
  annotate(
    "text", x = 1997, y = 40, 
    label = "NL adopts DH", hjust = "left"
  ) +
  labs(x = "Year", y = "Home runs", color = "League")

---

• 1973년 AL이 지명타자를 도입한 이후에는 AL이 우위를 점하고 있는 경향 (NL도 2022년에 도입)

시뮬레이션

반복과 분포

• 지금까지 벡터의 요소를 반복하면서 연산을 반복하는 방법을 배움
• 이 연산에 무작위성이 포함된 경우, 매번 같은 결과가 나오지 않으므로 무작위 연산이 생성하는 값의 분포를 이해하는 것이 유용

Bill James 기대승률 모형

• 1954년부터 2022년까지 69개 시즌 데이터를 사용하여 모형의 최적 지수가 1.84임을 계산

• 그러나 매 10년대별로 모형을 적합하면 지수가 1.69부터 1.94까지 변동성을 띔

• 10년대별 구분은 임의적

시즌별 최적 지수의 분포는 어떻게 되는가? 그 추정치인 1.84는 얼마나 신뢰할 수 있는가?

---

• group_modify()와 fit_k() 재활용

k_actual <- TeamRuns %>% 
  group_by(yearID) %>%         # by season
  group_modify(~fit_k(.))
k_actual %>%
  ungroup() %>%
  skim(k)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
k	69	0	1.84	0.18	1.31	1.69	1.88	1.96	2.33

---

ggplot(data = k_actual, aes(x = k)) + 
  geom_density() + 
  xlab("Best fit exponent for a single season")

---

• 표본 크기가 69개 뿐이므로, 평균의 표본분포를 더 잘 이해하기 위해 10000번의 복원추출을 통한 재표집을 사용 (부트스트랩)

• 수행하려는 반복 횟수를 n으로 정의하고, 단일 재표집 표본의 평균을 계산하는 식을 작성한 다음, map_dbl()을 1:n에 적용하여 반복을 수행

n <- 10000

bstrap <- 1:n %>%
  map_dbl(
    ~k_actual %>%
      pull(k) %>%
      sample(replace = TRUE) %>% 
      mean()
  )

civals <- bstrap %>%
  quantile(probs = c(0.025, .975))
civals

    2.5%    97.5% 
1.798625 1.883928 

---

ggplot(data = enframe(bstrap, value = "k"), aes(x = k)) + 
  geom_density() + 
  xlab("Distribution of resampled means") + 
  geom_vline(
    data = enframe(civals), aes(xintercept = value), 
    color = "red", linetype = 3
  ) +
  geom_vline(xintercept=1.84)

• 95%의 재표집된 지수가 1.7986247에서 1.8839279 사이에 있음
• 따라서, 초기 추정치인 1.84가 분포 중심 근방에 있음

산점도 프로그래밍

NHANES 데이터

• 체질량지수(body mass index, BMI)는 흔히 사용되는 비만 척도로서, 몸무게를 키의 제곱으로 나눈 비율

높은 BMI와 관련된 요인은 무엇인가?

• 미국 국립보건통계센터(National Center for Health Statistics)에서 수집한 국립보건 및 영양조사( National Health and Nutrition Examination Survey, NHANES) 데이터를 사용해서 알아본다.

str(NHANES::NHANES)

tibble [10,000 × 76] (S3: tbl_df/tbl/data.frame)
 $ ID              : int [1:10000] 51624 51624 51624 51625 51630 51638 51646 51647 51647 51647 ...
 $ SurveyYr        : Factor w/ 2 levels "2009_10","2011_12": 1 1 1 1 1 1 1 1 1 1 ...
 $ Gender          : Factor w/ 2 levels "female","male": 2 2 2 2 1 2 2 1 1 1 ...
 $ Age             : int [1:10000] 34 34 34 4 49 9 8 45 45 45 ...
 $ AgeDecade       : Factor w/ 8 levels " 0-9"," 10-19",..: 4 4 4 1 5 1 1 5 5 5 ...
 $ AgeMonths       : int [1:10000] 409 409 409 49 596 115 101 541 541 541 ...
 $ Race1           : Factor w/ 5 levels "Black","Hispanic",..: 4 4 4 5 4 4 4 4 4 4 ...
 $ Race3           : Factor w/ 6 levels "Asian","Black",..: NA NA NA NA NA NA NA NA NA NA ...
 $ Education       : Factor w/ 5 levels "8th Grade","9 - 11th Grade",..: 3 3 3 NA 4 NA NA 5 5 5 ...
 $ MaritalStatus   : Factor w/ 6 levels "Divorced","LivePartner",..: 3 3 3 NA 2 NA NA 3 3 3 ...
 $ HHIncome        : Factor w/ 12 levels " 0-4999"," 5000-9999",..: 6 6 6 5 7 11 9 11 11 11 ...
 $ HHIncomeMid     : int [1:10000] 30000 30000 30000 22500 40000 87500 60000 87500 87500 87500 ...
 $ Poverty         : num [1:10000] 1.36 1.36 1.36 1.07 1.91 1.84 2.33 5 5 5 ...
 $ HomeRooms       : int [1:10000] 6 6 6 9 5 6 7 6 6 6 ...
 $ HomeOwn         : Factor w/ 3 levels "Own","Rent","Other": 1 1 1 1 2 2 1 1 1 1 ...
 $ Work            : Factor w/ 3 levels "Looking","NotWorking",..: 2 2 2 NA 2 NA NA 3 3 3 ...
 $ Weight          : num [1:10000] 87.4 87.4 87.4 17 86.7 29.8 35.2 75.7 75.7 75.7 ...
 $ Length          : num [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ HeadCirc        : num [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ Height          : num [1:10000] 165 165 165 105 168 ...
 $ BMI             : num [1:10000] 32.2 32.2 32.2 15.3 30.6 ...
 $ BMICatUnder20yrs: Factor w/ 4 levels "UnderWeight",..: NA NA NA NA NA NA NA NA NA NA ...
 $ BMI_WHO         : Factor w/ 4 levels "12.0_18.5","18.5_to_24.9",..: 4 4 4 1 4 1 2 3 3 3 ...
 $ Pulse           : int [1:10000] 70 70 70 NA 86 82 72 62 62 62 ...
 $ BPSysAve        : int [1:10000] 113 113 113 NA 112 86 107 118 118 118 ...
 $ BPDiaAve        : int [1:10000] 85 85 85 NA 75 47 37 64 64 64 ...
 $ BPSys1          : int [1:10000] 114 114 114 NA 118 84 114 106 106 106 ...
 $ BPDia1          : int [1:10000] 88 88 88 NA 82 50 46 62 62 62 ...
 $ BPSys2          : int [1:10000] 114 114 114 NA 108 84 108 118 118 118 ...
 $ BPDia2          : int [1:10000] 88 88 88 NA 74 50 36 68 68 68 ...
 $ BPSys3          : int [1:10000] 112 112 112 NA 116 88 106 118 118 118 ...
 $ BPDia3          : int [1:10000] 82 82 82 NA 76 44 38 60 60 60 ...
 $ Testosterone    : num [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ DirectChol      : num [1:10000] 1.29 1.29 1.29 NA 1.16 1.34 1.55 2.12 2.12 2.12 ...
 $ TotChol         : num [1:10000] 3.49 3.49 3.49 NA 6.7 4.86 4.09 5.82 5.82 5.82 ...
 $ UrineVol1       : int [1:10000] 352 352 352 NA 77 123 238 106 106 106 ...
 $ UrineFlow1      : num [1:10000] NA NA NA NA 0.094 ...
 $ UrineVol2       : int [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ UrineFlow2      : num [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ Diabetes        : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ DiabetesAge     : int [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ HealthGen       : Factor w/ 5 levels "Excellent","Vgood",..: 3 3 3 NA 3 NA NA 2 2 2 ...
 $ DaysPhysHlthBad : int [1:10000] 0 0 0 NA 0 NA NA 0 0 0 ...
 $ DaysMentHlthBad : int [1:10000] 15 15 15 NA 10 NA NA 3 3 3 ...
 $ LittleInterest  : Factor w/ 3 levels "None","Several",..: 3 3 3 NA 2 NA NA 1 1 1 ...
 $ Depressed       : Factor w/ 3 levels "None","Several",..: 2 2 2 NA 2 NA NA 1 1 1 ...
 $ nPregnancies    : int [1:10000] NA NA NA NA 2 NA NA 1 1 1 ...
 $ nBabies         : int [1:10000] NA NA NA NA 2 NA NA NA NA NA ...
 $ Age1stBaby      : int [1:10000] NA NA NA NA 27 NA NA NA NA NA ...
 $ SleepHrsNight   : int [1:10000] 4 4 4 NA 8 NA NA 8 8 8 ...
 $ SleepTrouble    : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 1 1 1 ...
 $ PhysActive      : Factor w/ 2 levels "No","Yes": 1 1 1 NA 1 NA NA 2 2 2 ...
 $ PhysActiveDays  : int [1:10000] NA NA NA NA NA NA NA 5 5 5 ...
 $ TVHrsDay        : Factor w/ 7 levels "0_hrs","0_to_1_hr",..: NA NA NA NA NA NA NA NA NA NA ...
 $ CompHrsDay      : Factor w/ 7 levels "0_hrs","0_to_1_hr",..: NA NA NA NA NA NA NA NA NA NA ...
 $ TVHrsDayChild   : int [1:10000] NA NA NA 4 NA 5 1 NA NA NA ...
 $ CompHrsDayChild : int [1:10000] NA NA NA 1 NA 0 6 NA NA NA ...
 $ Alcohol12PlusYr : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 2 2 2 ...
 $ AlcoholDay      : int [1:10000] NA NA NA NA 2 NA NA 3 3 3 ...
 $ AlcoholYear     : int [1:10000] 0 0 0 NA 20 NA NA 52 52 52 ...
 $ SmokeNow        : Factor w/ 2 levels "No","Yes": 1 1 1 NA 2 NA NA NA NA NA ...
 $ Smoke100        : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 1 1 1 ...
 $ Smoke100n       : Factor w/ 2 levels "Non-Smoker","Smoker": 2 2 2 NA 2 NA NA 1 1 1 ...
 $ SmokeAge        : int [1:10000] 18 18 18 NA 38 NA NA NA NA NA ...
 $ Marijuana       : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 2 2 2 ...
 $ AgeFirstMarij   : int [1:10000] 17 17 17 NA 18 NA NA 13 13 13 ...
 $ RegularMarij    : Factor w/ 2 levels "No","Yes": 1 1 1 NA 1 NA NA 1 1 1 ...
 $ AgeRegMarij     : int [1:10000] NA NA NA NA NA NA NA NA NA NA ...
 $ HardDrugs       : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 1 1 1 ...
 $ SexEver         : Factor w/ 2 levels "No","Yes": 2 2 2 NA 2 NA NA 2 2 2 ...
 $ SexAge          : int [1:10000] 16 16 16 NA 12 NA NA 13 13 13 ...
 $ SexNumPartnLife : int [1:10000] 8 8 8 NA 10 NA NA 20 20 20 ...
 $ SexNumPartYear  : int [1:10000] 1 1 1 NA 1 NA NA 0 0 0 ...
 $ SameSex         : Factor w/ 2 levels "No","Yes": 1 1 1 NA 2 NA NA 2 2 2 ...
 $ SexOrientation  : Factor w/ 3 levels "Bisexual","Heterosexual",..: 2 2 2 NA 2 NA NA 1 1 1 ...
 $ PregnantNow     : Factor w/ 3 levels "Yes","No","Unknown": NA NA NA NA NA NA NA NA NA NA ...

• BMI 반응변수에 대해 75개의 잠재적 설명 변수가 있음

---

• BMI와 나이 사이의 산점도를 그리고 지역 회귀선을 추가하여 일반적인 추세를 보여주는 것으로 시작해 보자.

ggplot(NHANES, aes(x = Age, y = BMI)) + geom_point() +  geom_smooth()

---

나이 포함 NHANES의 모든 변수에 대해 비슷한 산점도를 그리도록 프로그램할 수 있을까?

• 절차
  1. 변수 이름을 입력으로 받고 산점도를 반환하는 함수를 작성
  2. 변수 집합을 정의하고 map()을 사용하여 1에서 작성한 함수를 반복 적용

bmi_plot <- function(.data, x_var) {
  ggplot(.data, aes(y = BMI)) +
    aes_string(x = x_var) + 
    geom_jitter(alpha = 0.3) + 
    geom_smooth() + 
    labs(
      title = paste("BMI by", x_var),
      subtitle = "NHANES",
      caption = "US National Center for Health Statistics (NCHS)"
    )
}

• x aesthetic에 x_var라는 이름의 개체가 아니라 x_var라는 이름의 변수를 대응시키려면 aes_string() 함수 사용이 필요

---

bmi_plot(NHANES, "Age")

---

• map()

c("Age", "HHIncomeMid", "PhysActiveDays", 
  "TVHrsDay", "AlcoholDay", "Pulse") %>%
  map(bmi_plot, .data = NHANES) %>%
  patchwork::wrap_plots(nrow = 2, ncol = 3)

---

과제

• names()를 사용하여 76개 산점도 그리기
• across()를 사용하여 특정 조건을 충족하는 변수와 BMI간의 산점도 그리기

참고 자료

범함수(functional)

• 함수를 입력으로 받아 벡터를 출력으로 반환하는 함수
  • https://adv-r.hadley.nz/functionals.html
  • 함수적 프로그래밍(functional programming)의 기본
    • Python: https://docs.python.org/3/howto/functional.html

깔끔한 평가(tidy evaluation)

https://dplyr.tidyverse.org/articles/programming.html

• Tidyverse 전체에서 사용되는 특수한 유형의 비표준 평가
  1. 데이터 마스킹: arrange(), count(), filter(), group_by(), mutate(), summarise(). 데이터 변수를 환경 변수처럼 사용(즉, df$my_variable가 아닌 my_variable로 쓸 수 있음).
  2. 깔끔한 선택(tidy selection): across(), relocate(), rename(), select(), pull(). 변수의 위치, 이름 또는 유형을 기준으로 쉽게 변수를 선택 (예: starts_with("x") 또는 is.numeric)

R 환경 (environments, MDSR Appendix C.4)

• 표현식(expression)을 평가(evaluate)할 때 R은 환경(environment)에서 개체를 검색
• Global environment
  • 가장 일반적인 환경
  • 내용은 RStudio의 Environment 탭 또는 ls() 명령을 통해 볼 수 있음
  • 전역 환경에서 찾을 수 없는 개체에 액세스하려고 하면 오류 발생

---

NHANES 데이터 예제

nhanes_small <- NHANES::NHANES %>%
  select(ID, SurveyYr, Gender, Age, AgeMonths, Race1, Poverty)
glimpse(nhanes_small)

Rows: 10,000
Columns: 7
$ ID        <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 51647, 5164…
$ SurveyYr  <fct> 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_1…
$ Gender    <fct> male, male, male, male, female, male, male, female, female, …
$ Age       <int> 34, 34, 34, 4, 49, 9, 8, 45, 45, 45, 66, 58, 54, 10, 58, 50,…
$ AgeMonths <int> 409, 409, 409, 49, 596, 115, 101, 541, 541, 541, 795, 707, 6…
$ Race1     <fct> White, White, White, Other, White, White, White, White, Whit…
$ Poverty   <dbl> 1.36, 1.36, 1.36, 1.07, 1.91, 1.84, 2.33, 5.00, 5.00, 5.00, …

---

R의 개체 식별 범위(scope)

"ID"

[1] "ID"

• 단일 문자열 “ID”를 원소로 하는 길이 1의 문자형 벡터 생성

ID

Error in eval(expr, envir, enclos): object 'ID' not found

• 전역 환경에서 존재하지 않는 ID라는 개체를 검색 — 오류

nhanes_small %>%
  pull(ID) %>%    # access within a data frame
  summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  51624   56904   62160   61945   67039   71915 

• 전역환경에서 접근할 수 있는 nhanes_small 데이터프레임 내의 ID 변수에 올바르게 접근

---

• nhanes_small의 많은 변수가 요인형
• 각 변수를 문자형으로 변환하고 싶을 때, mutate() 함수 사용

nhanes_small %>%
  mutate(SurveyYr = as.character(SurveyYr)) %>%
  select(ID, SurveyYr) %>% glimpse()

Rows: 10,000
Columns: 2
$ ID       <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 51647, 51647…
$ SurveyYr <chr> "2009_10", "2009_10", "2009_10", "2009_10", "2009_10", "2009_…

• 이 구조에서는 변환하려는 변수의 이름(예: SurveyYr)을 알고 이를 명시적으로 나열해야 함
• 우리의 목표는 데이터랭글링의 자동화!

---

• 열 이름(예: SurveyYr)을 변수 varname로 설정하고 해당 변수를 사용하여 이름 변경

varname <- "SurveyYr"
nhanes_small %>%
  mutate(varname = as.character(varname)) %>%
  select(ID, SurveyYr, varname) %>%
  glimpse()

Rows: 10,000
Columns: 3
$ ID       <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 51647, 51647…
$ SurveyYr <fct> 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_10…
$ varname  <chr> "SurveyYr", "SurveyYr", "SurveyYr", "SurveyYr", "SurveyYr", "…

• 왜 실패하는가?

---

• 이 경우 SurveyYr의 자료형을 변경하는 대신 SurveyYr의 값을 원소로 하는 문자 벡터인 varname이라는 새 변수를 만듦

• 이 동작은 비표준 평가(non-standard evaluation)라는 R 언어 기능의 결과

• dplyr 동사가 사용하는 깔끔한 평가(tidy evaluation)는 비표준 평가의 한 형태

  • 전역환경에 SurveyYr라는 개체가 없더라도 R이 이 개체를 찾을 수 있음
  • 앞 예제에서 mutate()는 nhanes_small 데이터프레임 내에서 SurveyYr를 찾아야 함을 알고 있음

---

해결책 1: dplyr 부사

• is.factor()를 across(), where()와 함께 사용하여 모든 요인형 변수를 찾은 다음 as.character()를 사용하여 문자로 변환

nhanes_small %>%
  mutate(across(where(is.factor), as.character))

# A tibble: 10,000 × 7
      ID SurveyYr Gender   Age AgeMonths Race1 Poverty
   <int> <chr>    <chr>  <int>     <int> <chr>   <dbl>
 1 51624 2009_10  male      34       409 White    1.36
 2 51624 2009_10  male      34       409 White    1.36
 3 51624 2009_10  male      34       409 White    1.36
 4 51625 2009_10  male       4        49 Other    1.07
 5 51630 2009_10  female    49       596 White    1.91
 6 51638 2009_10  male       9       115 White    1.84
 7 51646 2009_10  male       8       101 White    2.33
 8 51647 2009_10  female    45       541 White    5   
 9 51647 2009_10  female    45       541 White    5   
10 51647 2009_10  female    45       541 White    5   
# ℹ 9,990 more rows

---

해결책 2: base R

var_to_char_base <- function(data, varname) {
   data[[varname]] <- as.character(data[[varname]])
   data
}
var_to_char_base(nhanes_small, "SurveyYr")

# A tibble: 10,000 × 7
      ID SurveyYr Gender   Age AgeMonths Race1 Poverty
   <int> <chr>    <fct>  <int>     <int> <fct>   <dbl>
 1 51624 2009_10  male      34       409 White    1.36
 2 51624 2009_10  male      34       409 White    1.36
 3 51624 2009_10  male      34       409 White    1.36
 4 51625 2009_10  male       4        49 Other    1.07
 5 51630 2009_10  female    49       596 White    1.91
 6 51638 2009_10  male       9       115 White    1.84
 7 51646 2009_10  male       8       101 White    2.33
 8 51647 2009_10  female    45       541 White    5   
 9 51647 2009_10  female    45       541 White    5   
10 51647 2009_10  female    45       541 White    5   
# ℹ 9,990 more rows

• 주의사항
  1. [[ 같은 선택기를 사용하면 파이프라인이 끊어질 수 있음
  2. 변수 이름 주위의 따옴표가 번거로울 수 있음