PART IV Model

CHAPTER18 _남은정

library(modelr)
library(tidyverse)
options(na.action=na.warn)
sim1
ggplot(sim1,aes(x,y))+ geom_point()

models <- tibble(
  a1 = runif(250, -20, 40),
  a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) + geom_point() + geom_abline(aes(intercept=a1,slope=a2),data=models,alpha=1/4)
model1 <- function(a, data) { a[1] + data$x * a[2] }  # input: beta, data / output: predicted y
model1(c(7, 1.5), sim1)                               # predicted y 
measure_distance <- function(mod, data) {
  diff <- data$y - model1(mod, data)      # y-predicted y : distance
  sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
sim1_dist <- function(a1, a2) {         # input : beta1, beta2
  measure_distance(c(a1, a2), sim1)     # output: distance 
}
models <- models %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))   # dist 변수 생성 
models
ggplot(sim1, aes(x,y))+ geom_point(size=2,color="grey30")+ 
  geom_abline(aes(intercept=a1,slope=a2,color=-dist),data=filter(models,rank(dist)<=10))
                                           #  dist : dist가 작을수록 색이 어두움
                                           # -dist : dist가 작을수록 색이 밝음
ggplot(models, aes(a1, a2)) +
  geom_point(data = filter(models, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist))
grid <- expand.grid(                 # expand.grid : vector or factor 조합으로 data frame 생성
  a1 = seq(-5, 20, length = 25),
  a2 = seq(1, 3, length = 25)
  ) %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
grid %>% 
  ggplot(aes(a1, a2)) +
  geom_point(data = filter(grid, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist)) 
ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(grid, rank(dist) <= 10)
  )
(best <- optim(c(0, 0), measure_distance, data = sim1))
best$par
ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(intercept = best$par[1], slope = best$par[2])
(sim1_mod <- lm(y ~ x, data = sim1))
coef(sim1_mod)                  # optim 결과와 거의 일치
## Visualizing Models
sim1 
(grid <-sim1 %>%
 data_grid(x))      # sim1에는 20개의 obs 존재. data_grid(x) : unique한 x를 출력, 10개 출력
sim1_mod
(grid<-grid %>%
 add_predictions(sim1_mod))     # lm model의 predicted y 생성
ggplot(sim1, aes(x)) +
  geom_point(aes(y = y)) +
  geom_line(aes(y = pred), data = grid, colour = "red", size = 1)
(sim1 <- sim1 %>% 
  add_residuals(sim1_mod))
ggplot(sim1, aes(resid)) + 
  geom_freqpoly(binwidth = 0.5)
ggplot(sim1, aes(x, resid)) + 
  geom_ref_line(h = 0) +
  geom_point() 
df <- tribble(
  ~y, ~x1, ~x2,
  4, 2, 5,
  5, 1, 6
)
model_matrix(df, y ~ x1)        # Model Matrix 생성 (X) 
model_matrix(df, y ~ x1 - 1)        # intercept (1) 제외한 Matrix
model_matrix(df, y ~ x1 + x2)
(df <- tribble(
  ~ sex, ~ response,
  "male", 1,
  "female", 2,
  "male", 1
))
model_matrix(df, response ~ sex)    # Categorical variable은 1,0으로 표시
ggplot(sim2) + 
  geom_point(aes(x, y))
mod2 <- lm(y ~ x, data = sim2)
grid <- sim2 %>% 
  data_grid(x) %>% 
  add_predictions(mod2)
grid
ggplot(sim2, aes(x)) + 
  geom_point(aes(y = y)) +
  geom_point(data = grid, aes(y = pred), colour = "red", size = 4)
grid
mod2
tibble(x = "e") %>% 
  add_predictions(mod2)         # obs에 없는 level "e"에 대해서는 predict 불가능
sim3
ggplot(sim3, aes(x1, y)) + 
  geom_point(aes(colour = x2))
(mod1 <- lm(y ~ x1 + x2, data = sim3))
(mod2 <- lm(y ~ x1 * x2, data = sim3))  # x1 * x2 : interaction까지 포함
grid <- sim3 %>%                        # 120 rowS
          data_grid(x1,x2) %>%          # 40  rows  (unique한 x1, x2)
        gather_predictions(mod1,mod2)       # 80  rows  (mod1에 대해서 40개, mod2에 대해서 40개의 predictions)
grid
ggplot(sim3, aes(x1,y,color=x2)) + geom_point() + geom_line(data=grid, aes(y=pred))+ facet_wrap(~model)
sim3 <- sim3 %>% 
  gather_residuals(mod1, mod2)      # mod1, mod2 각각에 대해 residual 생성
ggplot(sim3, aes(x1, resid, colour = x2)) + 
  geom_point() + 
  facet_grid(model ~ x2)
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)
grid <- sim4 %>% 
  data_grid(
    x1 = seq_range(x1, 5), 
    x2 = seq_range(x2, 5) 
  ) %>% 
  gather_predictions(mod1, mod2)
grid
seq_range(c(0.0123, 0.923423), n = 5)
seq_range(c(0.0123, 0.923423), n = 5, pretty = TRUE)        # pretty : 반올림해서 보기 좋게 해줌 (table생성시 유용)
x1 <- rcauchy(100)
seq_range(x1, n = 5)                # seq_range  :min-max 구간을 5등분한 값
seq_range(x1, n = 5, trim = 0.10)       # trim: head, tail 각 5%를 제외한 후 시행
seq_range(x1, n = 5, trim = 0.25)       
seq_range(x1, n = 5, trim = 0.50)
x2 <- c(0, 1)
seq_range(x2, n = 5)
seq_range(x2, n = 5, expand = 0.10)     # expand : 양 끝을 5%씩 확장함
seq_range(x2, n = 5, expand = 0.25)
seq_range(x2, n = 5, expand = 0.50)
ggplot(grid, aes(x1, x2)) + 
  geom_tile(aes(fill = pred)) + 
  facet_wrap(~ model)
ggplot(grid, aes(x1, pred, colour = x2, group = x2)) + 
  geom_line() +
  facet_wrap(~ model)
ggplot(grid, aes(x2, pred, colour = x1, group = x1)) + 
  geom_line() +
  facet_wrap(~ model)

## Transformation
df <- tribble(
  ~y, ~x,
   1,  1,
   2,  2, 
   3,  3
)
df
model_matrix(df, y ~ x^2 + x)
model_matrix(df, y ~ I(x^2) + x)    # x^2 항 추가
model_matrix(df, y ~ poly(x, 2))
library(splines)
model_matrix(df, y ~ ns(x, 2))  # Generate the B-spline basis matrix for a natural cubic spline
sim5 <- tibble(
  x = seq(0, 3.5 * pi, length = 50),
  y = 4 * sin(x) + rnorm(length(x))
)
ggplot(sim5, aes(x, y)) +
  geom_point()
mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)

grid <- sim5 %>% 
  data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>% 
  gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")
grid                                            # model 5개 각각에 대해서 x의 seq_range 값에 대해 predictions 생성
ggplot(sim5, aes(x, y)) + 
  geom_point() +
  geom_line(data = grid, colour = "red") +
  facet_wrap(~ model)
df <- tribble(
  ~x, ~y,
  1, 2.2,
  2, NA,
  3, 3.5,
  4, 8.3,
  NA, 10
)
mod <- lm(y ~ x, data = df)
mod <- lm(y ~ x, data = df, na.action = na.exclude)
nobs(mod)   # number of observations (NA 제외하고 3개)

CHAPTER19_문지수

##############19. Model Building############
library(tidyverse)
library(modelr)
options(na.action = na.warn)
library(nycflights13)
library(lubridate)
## Why are low quantity diamonds more expensive?
ggplot(diamonds, aes(cut, price)) + geom_boxplot()
ggplot(diamonds, aes(color, price)) + geom_boxplot() #J색일수록 색이 좋지 않음
ggplot(diamonds, aes(clarity, price)) + geom_boxplot() #l1일수록 선명도가 좋지 않음
## Price and carat
#무게(carat)가 가격을 결정하는 중요한 요인이며, 낮은 품질의 다이아몬드가 무거운 경향이 있음.
install.packages("hexbin")
library(hexbin)
ggplot(diamonds, aes(carat, price)) + 
  geom_hex(bins = 50) #bins: 수직, 수평방향으로 bins의 개수
# carat의 영향을 분리하는 모델을 적합시켜봄으로써 
#다이아몬드의 다른 요소들이 상대적인 가격에 어떻게 영향을 미치는지 살펴볼 수 있음.
##데이터 작업 2가지를 먼저 하고 분석을 하겠다.->diamonds2
#Focus on diamonds smaller than 2.5 carats (99.7% of the data)
#Log-transform the carat and price variables. -> 선형관계로 만들 수 있음.
diamonds2 <- diamonds %>% 
  filter(carat <= 2.5) %>% 
  mutate(lprice = log2(price), lcarat = log2(carat))  
ggplot(diamonds2, aes(lcarat, lprice)) + 
  geom_hex(bins = 50)  
mod_diamond <- lm(lprice ~ lcarat, data = diamonds2)  
diamonds2 %>% 
  data_grid(carat = seq_range(carat, 20))  
diamonds2 %>% 
  data_grid(carat = seq_range(carat, 20)) %>% 
  mutate(lcarat = log2(carat))  
diamonds2 %>% 
  data_grid(carat = seq_range(carat, 20)) %>% 
  mutate(lcarat = log2(carat)) %>% 
  add_predictions(mod_diamond, "lprice")
grid <- diamonds2 %>% 
  data_grid(carat = seq_range(carat, 20)) %>% 
  mutate(lcarat = log2(carat)) %>% 
  add_predictions(mod_diamond, "lprice") %>% 
  mutate(price = 2 ^ lprice)
ggplot(diamonds2, aes(carat, price)) + 
  geom_hex(bins = 50) + 
  geom_line(data = grid, colour = "red", size = 1)
#적합시킨 모형에 따르면, 무거운 다이아몬드는 예상보다 비싸지 않다는 것을 볼 수 있음.
#(이 데이터셋에는 19000달러보다 비싼 다이아몬드는 없기 때문)
diamonds2 <- diamonds2 %>% 
  add_residuals(mod_diamond, "lresid")
ggplot(diamonds2, aes(lcarat, lresid)) + 
  geom_hex(bins = 50)
#잔차가 특별한 패턴을 보이지 않아 적절함을 알 수 있음.
##이제 가격(price)대신에 잔차를 사용하여 플롯을 다시 그려보자
ggplot(diamonds2, aes(cut, lresid)) + geom_boxplot()
ggplot(diamonds2, aes(color, lresid)) + geom_boxplot()
ggplot(diamonds2, aes(clarity, lresid)) + geom_boxplot()
#이를 통해 다이아몬드의 품질이 좋을수록 상대적인 가격도 높아짐을 볼 수 있음.
#잔차가 -1인 경우, lprice가 다이아몬드 무게에 근거한 예측치보다 1단위 낮다는 의미임.
## A more complicated model
mod_diamond2 <- lm(lprice ~ lcarat + color + cut + clarity, data = diamonds2)
diamonds2 %>% 
  data_grid(cut, .model = mod_diamond2) 
grid <- diamonds2 %>% 
  data_grid(cut, .model = mod_diamond2) %>% 
  add_predictions(mod_diamond2)
grid
ggplot(grid, aes(cut, pred)) + 
  geom_point()
diamonds2 <- diamonds2 %>% 
  add_residuals(mod_diamond2, "lresid2")
ggplot(diamonds2, aes(lcarat, lresid2)) + 
  geom_hex(bins = 50)
#몇몇의 다이아몬드에서 잔차가 크게 나타나고 있음을 볼 수 있음.
diamonds2 %>% 
  filter(abs(lresid2) > 1) %>% 
  add_predictions(mod_diamond2) %>% 
  mutate(pred = round(2 ^ pred)) %>% 
  select(price, pred, carat:table, x:z) %>% 
  arrange(price)  

## Exerciese
#4
diamonds2 %>% 
  add_predictions(mod_diamond2) %>%
  add_residuals(mod_diamond2) %>%
  summarise(sq_err = sqrt(mean(resid^2)),
            abs_err = mean(abs(resid)),
            p975_err = quantile(resid, 0.975),
            p025_err = quantile(resid, 0.025))

## What affects the number of  daily flights?
daily <- flights %>% 
  mutate(date = make_date(year, month, day)) %>% 
  group_by(date) %>% 
  summarise(n = n())
daily  ##the number of flights per day
ggplot(daily, aes(date, n)) + 
  geom_line()
## Day of week
daily <- daily %>% 
  mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wday, n)) + 
  geom_boxplot()
#주중보다 주말에 비행이 적은 경향을 보임.
mod <- lm(n ~ wday, data = daily)
grid <- daily %>% 
  data_grid(wday) %>% 
  add_predictions(mod, "n")
ggplot(daily, aes(wday, n)) + 
  geom_boxplot() +
  geom_point(data = grid, colour = "red", size = 4)
daily <- daily %>% 
  add_residuals(mod)
daily %>% 
  ggplot(aes(date, resid)) + 
  geom_ref_line(h = 0) + 
  geom_line()
#위의 모형(mod)에 요일 효과를 넣었기 때문에, 이제 위의 플롯(y축이 잔차)을 볼 때 요일과 관련된 영향은 제외된 것임.
#플롯을 보면 6월부터 이상함(잔차가 크게 나타남). 
#또한 6월 이후로 규칙적인 패턴이 보임. -> 우리는 이러한 패턴을 모형(mod)에 넣지 못한것임.
ggplot(daily, aes(date, resid, colour = wday)) + 
  geom_ref_line(h = 0) + 
  geom_line()
#일요일을 보면 우리 모형이 잘 예측하지 못함을 볼 수 있음.
#여기서 주로 여름에는 예측보다 비행이 많고, 가을에는 예측보다 비행이 적게 나타남.
# Exercises
#3
daily <- daily %>%
  mutate(wday2 = 
         case_when(.$wday == "Sat" & .$term == "summer" ~ "Sat-summer",
         .$wday == "Sat" & .$ term == "fall" ~ "Sat-fall",
         .$wday == "Sat" & .$term == "spring" ~ "Sat-spring",
         TRUE ~ as.character(.$wday)))
mod4 <- lm(n ~ wday2, data = daily)
daily %>% 
  gather_residuals(sat_term = mod4, all_interact = mod2) %>% 
  ggplot(aes(date, resid, colour = model)) +
    geom_line(alpha = 0.75)
#4
daily <- daily %>%
  mutate(wday3 = 
         case_when(
           .$date %in% lubridate::ymd(c(20130101, # new years
                                        20130121, # mlk
                                        20130218, # presidents
                                        20130527, # memorial
                                        20130704, # independence
                                        20130902, # labor
                                        20131028, # columbus
                                        20131111, # veterans
                                        20131128, # thanksgiving
                                        20131225)) ~
             "holiday",
           .$wday == "Sat" & .$term == "summer" ~ "Sat-summer",
           .$wday == "Sat" & .$ term == "fall" ~ "Sat-fall",
           .$wday == "Sat" & .$term == "spring" ~ "Sat-spring",
           TRUE ~ as.character(.$wday)))

mod5 <- lm(n ~ wday3, data = daily)

daily %>% 
  spread_residuals(mod5) %>%
  arrange(desc(abs(resid))) %>%
  slice(1:20) %>% select(date, wday, resid)


#7
flights %>% 
  mutate(date = make_date(year, month, day),
         wday = wday(date, label = TRUE)) %>%
  group_by(wday) %>%
  summarise(dist_mean =  mean(distance),
            dist_median = median(distance)) %>%
  ggplot(aes(y = dist_mean, x = wday)) +
  geom_point()

#8
monday_first <- function(x) {
  forcats::fct_relevel(x, levels(x)[-1])  
}

aily <- daily %>% 
  mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(monday_first(wday), n)) + 
  geom_boxplot() +
  labs(x = "Day of Week", y = "Number of flights")

CHAPTER20_조현선

“Many Models with purr and broom”"

Hans Rosling 의 ted 강연보기 (한글 자막) https://youtu.be/NVwB_of8ZYs

powerful ideas

• 많은 수의 간단한 모델은 복잡한 데이터셋을 더 잘 이해하게 해준다.
• list-columns를 이용하여 데이터를 구조화시킨다.
• broom 패키지를 이용해서 tidy data를 만들고, 이로 모델링을 한다.

library(modelr)
library(tidyverse)
library(gapminder)

we will use the “gapminder”data.

gapminder

Q. [“How does life expectancy(lifeExp) change over time(year) for each country(country)?”]

: 3가지 변수에 초점을 맞추어 분석하고자 함.

gapminder %>% 
  ggplot(aes(year, lifeExp, group = country)) +
    geom_line(alpha = 1/3)

1700개의 관측치와 3개의 변수로 plotting해 본 결과, 많은 수가 아님에도 데이터의 흐름을 파악하기가 쉽지 않다. 전반적으로 증가하는 양상을 볼 수는 있으나, 몇 개의 국가들은 이를 따르지 않는 다는 것을 볼 수 있다. 따라서 “strong signal”을 사용해서 트랜드를 분석해보자. 선형 회귀를 적합시킨 것과 남은 residual에 대해 plotting을 해보자.

nz <- filter(gapminder, country == "New Zealand")
nz %>% 
  ggplot(aes(year, lifeExp)) + 
  geom_line() + 
  ggtitle("Full data = ")

nz_mod <- lm(lifeExp ~ year, data = nz)
nz %>% 
  add_predictions(nz_mod) %>%
  ggplot(aes(year, pred)) + 
  geom_line() + 
  ggtitle("Linear trend + ")

nz %>% 
  add_residuals(nz_mod) %>% 
  ggplot(aes(year, resid)) + 
  geom_hline(yintercept = 0, colour = "white", size = 3) + 
  geom_line() + 
  ggtitle("Remaining pattern")

new zealand 말고도, 모든 나라에 대해 쉽게 fitting 시키려면 어떻게 해야 할까? 방금과 같은 코드를 나라만 바꿔가며 복사 붙여넣기를 하는 방법도 있겠지만, 더 효율적인 방법이 있다.

Nested Data

데이터셋을 나라별로 뽑아서 반복을 하면 되는데, 이를 위해 자료의 구조를 새롭게 정의해야 한다. 이를 nested data frame이라고 한다. 이를 위한 코드는 다음과 같다.

by_country <- gapminder %>% 
  group_by(country, continent) %>% 
  nest()
by_country

변형된 자료를 파악하기 위해서, 하나의 원소를 골라서 어떤 정보를 어떻게 담고 있는지를 확인해보자. (Afghanistan을 고른 예시)

by_country$data[[1]]

위에서 처럼, nested data frame에서는 각 열 자체가 그룹임을 확인 할 수 있다. 이제, 이 nested data frame으로 모델을 fitting 시켜 보는 함수를 만들어보자.

country_model <- function(df) {
  lm(lifeExp ~ year, data = df)
}

이 함수와 데이터를 이제, purrr패키지의 map 함수를 이용해서 연결시켜보자.

library(purrr)
models <- map(by_country$data, country_model)

그리고, 이 결과를 다른 새로운 object로 만들어서 global environment에 저장하는 것이 아니라, by_country라는 데이터 프레임에 새로운 변수를 만들어서 저장하는 방식을 택하는 것이 더 낫다. 이를 위해 dplyr패키지의 mutate 함수를 이용한다.

by_country <- by_country %>% 
  mutate(model = map(data, country_model))
by_country

이렇게 작업을 하면, 모든 관련 개체가 함께 저장되므로 필터링하거나 정렬 할 때 수동으로 동기화 할 필요가 없다는 강력한 장점이 있다. 만약 이렇게 하지 않고, 별개의 객체로 저장을 해놓으면, 벡터를 재정의 하거나 벡터를 subset시킬때 하위집합을 계속 동기화해주어야 한다. 이를 잊어버린다면, 코드는 돌아는 가겠지만, 잘못된 답을 출력해주게 된다.

Unnesting

142개의 나라가 있으므로, 우리는 142개의 데이터프레임과 142개의 모델을 갖고 있다. 잔차를 계산하기 위해서, add_residuals()을 이용한다.

by_country <- by_country %>% 
  mutate(
    resids = map2(data, model, add_residuals)
  )
by_country

Q. How can you plot a list of data frames?

–> nested된 데이터 프레임을 unnest()를 이용하여 정규 데이터 프레임 형식으로 바꾼 후 plot을 그린다. 그럼 다음과 같이 나라별로 묶여있던 데이터프레임이 풀리게 됨을 확인할 수 있다.

resids <- unnest(by_country, resids)
resids

이제, 잔차도를 그릴 수 있게 되었다.

resids %>% 
  ggplot(aes(year, resid)) +
    geom_line(aes(group = country), alpha = 1 / 3)  +
    geom_smooth(se = FALSE)

대륙별로 보면, 대륙마다의 특별한 특징이 더 드러나게 된다.

resids %>% 
  ggplot(aes(year, resid, group = country)) +
    geom_line(alpha = 1 / 3) + 
    facet_wrap(~continent)

잔차그림의 모양을 보니, fitting 시킨 모형이 적절하지 않았다고 판단할 수 있다. 특히 africa에서 큰 잔차들을 볼 수 있다. 다음 섹션에서 이를 다른 방향으로 접근하여 분석해보자.

Model Quality

모델로부터 잔차를 구하여 보는 것보다, model quality에 대한 보편적인 측정치를 보는 방법에 대해 생각해보자. broom 패키지로, 모델을 tidy데이터로 만들 수 있다. 또한, broom 패키지의 glance()함수를 써서 모형 평가 지표를 데이터 프레임 형태로 뽑아낼 수 있다.

broom::glance(nz_mod)

여기서 구한 이 정보를 mutate()와 unnest()를 이용해서 각 나라별의 행으로 다시 넣은 데이터 프레임을 만들 수 있다.

by_country %>% 
  mutate(glance = map(model, broom::glance)) %>% 
  unnest(glance)

그러나 이렇게, 모든 list 열을 포함하는 것은 우리가 원하던 결과가 아니므로, unnest에서 .drop=TRUE로 다시 옵션을 설정하여 바꿔준다.

glance <- by_country %>% 
  mutate(glance = map(model, broom::glance)) %>% 
  unnest(glance, .drop = TRUE)
glance

이렇게 얻어진 데이터프레임으로 드디어, 잘 맞지 않는 모델이 무엇인지 찾을 수 있게 되었다. 모형 평가 지표중의 하나인 r square값에 대해 오름차순으로 sorting을 해보자.

glance %>% 
  arrange(r.squared)

그 결과, Africa에서 모델이 제일 안맞음을 확인 할 수 있었다. 이를 다시 재확인해보기 위해 그림을 그려보려 한다. 관측수가 적고, 이산형 변수임으로, 지터그림을 보는 것이 효과적이다.

glance %>% 
  ggplot(aes(continent, r.squared)) + 
    geom_jitter(width = 0.5)

우리는 다음과 같이 낮은 R-square값을 갖고 있는 국가들만을 가지고 다시 그림을 그려볼 수 있다.

bad_fit <- filter(glance, r.squared < 0.25)
gapminder %>% 
  semi_join(bad_fit, by = "country") %>% 
  ggplot(aes(year, lifeExp, colour = country)) +
    geom_line()

우리는 여기서 HIV / AIDS 전염병과 르완다 대학살, 이 두가지 비극의 효과를 그림으로 확인 할 수 있었다. (동영상 내용 참고)

=============================================================================================

List-Columns

지금까지 많은 모델을 다루는 것에 대한 전반적인 흐름을 보았다. 지금부터는 list-column 데이터 구조에 대해 더 자세히 살펴 볼 것이다. 데이터 프레임은 동일한 길이 벡터가 모여 명명되고, list-column은 데이터프레임의 정의안에 내포되어 있는 개념이다. list는 벡터이므로, list를 데이터 프레임의 열로 추가하는 것은 괜찮다. 그러나, data.frame()에서는 list를 열의 list로 다루므로, 기존의 R에서는 이러한 작업이 쉽지 않았다. 예를 통해 살펴보자.

data.frame(x = list(1:3, 3:5))

이렇게 하는 것은 우리가 원했던 결과가 아니므로, 이때는 I()를 사용하여 아래와 같이 나타내야 한다. 그러나 완전히 잘 출력해준다고는 볼 수 없다.

data.frame(
  x = I(list(1:3, 3:5)), 
  y = c("1, 2", "3, 4, 5")
)

tibble(
  x = list(1:3, 3:5), 
  y = c("1, 2", "3, 4, 5")
)

tibble은 입력값을 복사하는 개념이 아니므로, 더 나은 출력 방법을 제공하여 위의 예시들에서 발생한 문제들을 해결한다.

tribble(
   ~x, ~y,
  1:3, "1, 2",
  3:5, "3, 4, 5"
)

이렇게 하면 내가 필요로 하는 리스트를 더 쉽게 자동적으로 tibble을 이용하여 만들 수 있다. 이렇게, list-columns는 즉각적인 데이터 구조로 활용되기에 유용하다. 대부분의 R함수들은 atomic 벡터와 데이터 프레임을 이용해서 함수가 짜여져 있으므로, tibble을 이용해서 바로 작업을 수행하기는 어려울 수 있지만, 관련있는 항목들을 연결하여 쓸 수 있다는 장점이 있으므로 그만한 가치가 있다고 본다.

Creating List-Columns

통상적으로, tibble()을 이용해서 list-columns를 바로 만들지는 않는다. 대신에, regular columns를 이용해서 만드는데, 다음 세가지 방법들 중의 하나로 한다. - (1) tidyr::nest() 로 그룹화된 데이터 프레임을 nested 데이터 프레임으로 바꾼다. - (2) mutate() 와 벡터화 시키는 함수를 이용해 list를 반환한다. - (3) summerize() 와 요약해주는 함수를 통해 여러 결과를 출력한다. 대체적으로, tibble::enframe()을 사용하면, 명명된 리스트로 부터 만들어 낼 수도 있다. 일반적으로, list-columns를 만들때는, 그들이 homogeneous한 것을 확실하게 해야한다. 즉, 같은 요소를 담고 있어야 한다는 것이다. purrr에서 type-stable 함수를 배웠으므로, 이를 기억하고 자연스럽게 작업을 할 수 있을 것이다.

(1) With Nesting

nest()는 리스트 열을 갖는 nested 데이터 프레임을 만들어 준다. nest()를 쓰는 방법에는 두 가지가 있는데, 지금까지는 그룹화된 데이터 프레임에서 사용해왔었다. 이 때 nest()는 그룹화된 열을 그 그대로 유지시켜주면서 모든 것을 다발로 그 리스트 열 안에 담아 주었었다. 다음 예시를 확인해보면 이해가 될 것이다.

gapminder %>% 
  group_by(country, continent) %>% 
  nest()

우리는 이것을 그룹화되지 않았던 데이터 프레임에도 사용할 수 있고, 중첩 할 열을 특별히 지정할 수도 있다.

gapminder %>% 
  nest(year:gdpPercap)

(2) From Vectorized Functions

atomic 벡터를 리스트로 바꿔주는 유용한 함수들에 대해서 배웠었다(11장 참고). mutate를 사용하면, list-column을 만들 수 있다. 다음 예시를 확인해보자.

df <- tribble(
  ~x1,
  "a,b,c", 
  "d,e,f,g"
) 
df %>% 
 mutate(x2 = stringr::str_split(x1, ","))

이를 unnest()를 사용하면, x2에 할당된 리스트를 벡터로 바꿀 수 있게 된다.

df %>% 
  mutate(x2 = stringr::str_split(x1, ",")) %>% 
  unnest()

이러한 유형은 purrr 패키지의 map(), map2(), pmap()등을 사용하여 비슷한 작업을 할 수 있다.

(3) From Multivalued Summeries

summerize()의 제약조건은, 한 값에 대해서만 출력을 해주는 함수라는 것이다. 쉽게 말해, 우리는 이것을 임의의 길이를 가진 벡터 즉, quantile()같은 함수에는 사용을 할 수 없다는 것이다. 때문에 다음 코드에서 에러가 나타난다.

mtcars %>% 
  group_by(cyl) %>% 
  summarise(q = quantile(mpg))

Error in summarise_impl(.data, dots) : expecting a single value

에러가 나지 않게 하려면 아래처럼 결과를 list로 묶어서 돌리면 돌아가게 된다. 왜냐하면 summary는 길이가 1이 된 리스트에 적용되기 때문이다.

mtcars %>% 
  group_by(cyl) %>% 
  summarise(q = list(quantile(mpg)))

이 결과를 unnest()를 이용하여 다시 풀어주면, 우리가 원하던 결과와 확률값을 얻게 된다.

probs <- c(0.01, 0.25, 0.5, 0.75, 0.99)
mtcars %>% 
  group_by(cyl) %>% 
  summarise(p = list(probs), q = list(quantile(mpg, probs))) %>% 
  unnest()

(4) From a Named List

list열을 갖고 있는 데이터 프레임은 “목록의 내용과 요소를 모두 반복하려는 경우 어떻게 해야 하나?”라는 흔한 문제에 대한 해결책을 제시한다. 모든 것을 하나의 객체로 묶는 대신 데이터 프레임을 만드는 것이 더 쉽다. 한 열에는 element가 포함될 수 있고 또 한 열에는 list가 포함될 수 있다. list로부터 그러한 데이터 프레임을 만드는 쉬운 방법은 tibble :: enframe ()를 쓰면 된다. 다음 예시를 보자.

x <- list(
  a = 1:5,
  b = 3:4, 
  c = 5:6
) 
df <- enframe(x)
df

이러한 구조의 장점은 직설적인 방식으로 일반화된다는 것이다. name은 문자형 변수가 meta data로 있으면 유용하지만, 다른 유형의 데이터나 여러개의 벡터가 이름에 들어가있으면 유용하지는 않은 표현방식이다. 어쨌든, 이제 이름과 값을 함께 반복하려면, map2()를 사용하면 된다.

df %>% 
  mutate(
    smry = map2_chr(name, value, ~ stringr::str_c(.x, ": ", .y[1]))
  )

List-Columns를 간단히 하기

• (1). 리스트를 벡터화 시키기 : 단일 값을 원할 경우 map_lgl (), map_int (), map_dbl () 및 map_chr ()에 mutate ()를 사용하여 원자 벡터를 만든다.

df <- tribble(
  ~x,
  letters[1:5],
  1:3,
  runif(5)
)
  
df %>% mutate(
  type = map_chr(x, typeof),
  length = map_int(x, length)
)

• (2). Unnesting : 많은 값을 원하면 unnest ()를 사용하여 목록 열을 일반 열로 다시 변환하고 필요한만큼 행을 반복한다. 예시1

tibble(x = 1:2, y = list(1:4, 1)) %>% unnest(y)

예시2

df1 <- tribble(
  ~x, ~y,           ~z,
   1, c("a", "b"), 1:2,
   2, "c",           3
)
df1

df1 %>% unnest(y, z)

Broom을 이용해 Tidy data 만들기

• broom::glance(model)
• broom::tidy(model)
• broom::augment(model, data)

Broom은 가장 인기인는 모델링 패키지들을 이용해 다양한 모델들로 작동한다. 최근에 제공해주는 모델들에 대한 정보는 다음의 깃허브를 확인해보자. https://github.com/tidyverse/broom