18장 modelr을 이용한 모델의 기초

18.1 들어가기

모델의 목표는 데이터셋에 대한 낮은 차원의 간단한 요약을 제공하는 것이며 모델에는 2가지 부분이 있다.

  1. 먼저 수집하고자 하는 정확하면서 일반적인 패턴을 표현하는 모델 모음을 정의한다. 예를 들어 패턴이 직선 또는 2차 곡선이라고 하면 y = a_1 * x + a_2 또는 y = a_1 * x ^ a_2와 같은 방적식으로 모델 모음을 표현할 것이다. 여기서 x와y는 데이터에서 명시된 변수이며 a_1과 a_2 포착된 패턴들에 따른 다양한 파라미터 값이다.
  2. 다음으로 모델 모음에서 데이터와 가장 가까운 모델을 찾아 적합한 모델(fitted model)을 생성한다. 일반적인 모델 모음에서 y = 3 * x +7 또는 y = 9 * x ^ 2과 같이 구체적으로 만든다.

그리고 또 다른 모델의 목표는 진실을 알아내는 것이 아니라 유용하면서 간단한 근사치를 알아내는 것이다.

18.1.1준비하기

이 장에서는 파이프에서 자연스럽게 작동하기 위해 베이스 R의 모델링 함수를 둘러싸는 modelr패키지를 사용할 것이다.

library(tidyverse)
library(modelr)
options(na.action = na.warn)

18.2 간단한 모델

시뮬레이션된 데이터셋 siml을 살표보자. 이 데이터셋에는 두 개의 연속 변수 x와 y가 포함되어 있다. 변수들의 관계를 살펴보기 위해 플롯을 나타내보자.

ggplot(sim1, aes(x, y)) +
  geom_point()

  • 위 예제는 선형관계(즉, y = a_0 + a_1 * x)인 것으로 보인다.

  • 이처럼 간단한 경우에는 기울기와 y절편을 파라미터로 사용하는 geom_abline()을 사용할 수 있고 뒷부분에서는 모든 모델에 동작하는 좀 더 일반적인 기법에 대해 배울 것이다.
    (geom_abline() : 이 함수는 기준선을 생성합니다.)

models <- tibble(
  a1 = runif(250, -20, 40),
  a2 = runif(250, -5, 5)
)

(runif : 이 함수는 최소에서 최대까지의 간격에 대한 균등 분포에 대한 정보를 제공합니다. dunif는 밀도를 제공하고 punif는 분포 함수를 제공합니다. qunif는 분위수 함수를 제공하고 runif는 무작위 편차를 생성합니다.)

ggplot(sim1, aes(x, y)) +
  geom_abline(
    aes(intercept = a1, slope = a2),
    data = models, alpha = 1/4
  ) +
  geom_point()

  • 이 플롯에는 250개의 모델이 존재하지만 실제로 대부분은 잘 맞지 않는다. 좋은 모델은 데이터와 ‘유사’하기에, 이러한 직관을 명확하게 하기 위해서라도 좋은 모델을 찾아야만 한다. 그런 다음 데이터와 가장 작은 차이를 가진 모델을 생성하는 a_0과 a_1 값을 찾아 모델을 적합할 수 있다.

  • 이를 위한 한 가지 쉬운 방법은 아래 도표에서처럼 각 점과 모델 사이의 수직거리를 계산하는 것이다. (x 값을 조금씩 이동하여 개별 거리를 확인할 수 있도록 하였다.)

  • 이 차이는 모델에 의해 주어진 y 값(예측값)과 데이터의 실제 y 값(반응값) 간의 차이일 뿐이다.
  • 이 차이를 계산하기 위해 먼저 모델 모음을 R 함수로 변환해야 한다. 이것은 모델 파라미터와 데이터를 입력값으로 사용하고 모델에 의해 예측된 값을 출력값으로 제공한다.
model1 <- function(a, data) {
  a[1] + data$x * a[2]
}
model1(c(7, 1.5), sim1)
##  [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5
## [16] 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0 22.0 22.0

  • 다음으로 예측된 값과 실제 값 사이의 전체 차이를 계산하는 방법이 필요하다. 다시 말하자면 플롯에서 나타내는 30개의 차이 값을 어떻게 하면 하나의 값으로 합칠 수 있을까?

  • 통계학에서 이를 수행하는 일반적인 방법은 ‘평균제곱근 편차’를 사용하는 것이다. 실제값과 예측값의 차이를 제곱한 후 이 값들의 평균에 제곱근을 취한다.

measure_distance <- function(mod, data) {
  diff <- data$y - model1(mod, data)
  sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
## [1] 2.665212
  • 이제 이전에 정의된 모든 모델의 차이를 계산하기 위해 purrr를 사용할 수 있다. 거리 차이를 계산하는 함수에서는 모델이, 길이가 2인 수치형 벡터 형식으로 입력되어야 하므로 도우미 함수가 필요하다.
    (purrr : 함수형 프로그래밍 도구)
sim1_dist <- function(a1, a2) {
  measure_distance(c(a1, a2), sim1)
}
models <- models %>%
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models
## # A tibble: 250 x 3
##        a1      a2  dist
##     <dbl>   <dbl> <dbl>
##  1  29.9   3.16   32.0 
##  2  26.9  -2.78   14.5 
##  3  -5.94 -2.55   37.9 
##  4  34.2   1.07   24.9 
##  5  18.7   2.58   17.6 
##  6  13.6  -0.279   7.81
##  7  37.1   4.12   44.7 
##  8  -2.98  0.880  14.2 
##  9 -11.0   0.0244 27.1 
## 10  29.2   0.0862 15.4 
## # ... with 240 more rows
  • 다음으로 최적의 10가지 모델을 데이터에 겹쳐서 나타내보자, -dist를 사용해서 모델별로 색상을 지정했다. 이는 최적의 모델 (즉, 가장 차이가 적은 모델)을 가장 밝은 색상으로 나타낼 수 있는 쉬운 방법이다.
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로 동일하게 색상을 나타낸 a1과 a2의 산점도로 시각화할 수 있다. 데이터와 직접 비교한 모델은 확인할 수 없지만 한번에 여러 모델을 확인할 수 있다. 다시 말해 이전에는 가장 좋은 10개 모델만을 나타냈지만, 지금은 빨간색 원으로 10개의 모델을 강조했다.
ggplot(models, aes(a1, a2)) +
  geom_point(
    data = filter(models, rank(dist) <= 10),
    size = 4, color = "red"
  ) + 
  geom_point(aes(colour = -dist))

  • 랜덤한 모델을 많이 시도해 보는 체계적이고 균등하게 배치된 격자무늬의점(그리드 서치(grid search)라고 함)을 생성할 수 있다. 앞의 플롯에서 가장 좋은 모델이 어디에 있는지 대략 살펴보며 그리드 파라미터를 선택한다.
grid <- expand.grid(
  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(color = -dist))

  • 다음 과 같이 원데이터에 최적의 10개 모델을 겹쳐 그리면 모두 잘 맞는 것처럼 보이다.
ggplot(sim1, aes(x, y)) +
  geom_point(size = 2, color = "grey30") +
  geom_abline(
    aes(intercept = a1, slope = a2, color = -dist),
    data = filter(grid, rank(dist) <= 10)
  )

최적의 모델을 선택할 떄까지 그리드를 반복적으로 더욱 세밀하게 만드는 작업을 생각해볼 수있다. 하지만 이 문제를 해결할 수 있는 더 좋은 방법인 뉴턴 랩슨 기법이라 불리는 수치 최소화 도구가 있다. 이는 시작적을 선택하고 가장 가파른 기울기를 찾기 위해 탐색한다. 그런 다음 가장 작은 값으로 갈 수 없을 때까지 기울기를 약간씩 기울이는 작업을 반복한다. R에서는 optim()을 사용하여 이 작업을 수행 할 수 있다.

best <- optim(c(0, 0), measure_distance, data = sim1)
best$par
## [1] 4.222248 2.051204
ggplot(sim1, aes(x, y)) +
  geom_point(size = 2, color = "grey30") +
  geom_abline(intercept = best$par[1], slope = best$par[2])

optim()의 세부적인 작동 방식에 대해서는 고려하지 않아도 된다. 여기서 중요한 것은 바로 직관이다. 모델과 데이터셋 간의 차이를 정의하는 함수와 모델의 파라미터를 변경함에 따라 거리를 최소화할 수 있는 알고리즘이 있다면 최적의 모델을 찾을 수 있다. 이 접근법의 휼률한 점은 방정식을 작성할 수 있는 모든 모델 모음에 작동한다는 것이다.
이 모델은 일반적인 모델 모음(즉, 신형 모델)의 특별한 경우이므로 이 모델에 적용할 수 있는 접근법이 한 가지 더 있다. 선형 모델은 y = a_1 + a_2 * x_1 + a_3 * x_2 + ... + a_n * x_(n-1)과 같은 일반적인 형식을 가진다. 따라서 이 간단한 모델은 n이 2이고 x_1이 x인 일반 선형 모델과 같다. R에는 lm()이라는 선형 모델을 적합하기 위해 특별히 고안된 툴이 있다. lm()에는 모델 집합을 지정하는 특별한 방법인 ‘수식’이 들어있다. 수식은 y ~ x와 같은 형태이며 lm()은 y =a_1 + a_2 * x와 같은 함수로 변환해준다. 다음과 같이 모델을 적합하여 결과를 확인할 수 있다.

sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)
## (Intercept)           x 
##    4.220822    2.051533

(coef: 모델링 함수에 의해 반환 된 객체에서 모델 계수를 추출하는 일반 함수입니다.)

위 값은 optim()을 사용해서 얻은 값과 정확히 일치하는 값이다. lm()의 이면에서는 optim()을 사용하지 않는 대신, 선형 모델의 수학적 구조의 이점을 사용한다. 기하학, 미적분학, 선형대수학 간의 연결을 기반으로 한 lm()은 정교한 알고리즘을 사용하여 개별 단계에서 가장 가까운 모델을 찾는다. 이 접근법은 더 빠르며 전체의 최솟값을 보장한다.

18.3 모델 시각화하기

이전 절과 같은 간단한 모델의 경우, 모델 집합과 적합된 계수를 주의 깊게 살펴보면 모델이 포착하는 패턴을 파악할 수 있다. 모델링에 관한 통계 과정을 수강하게 된다면 이 작업을 수행하는 데 많은 시간을 할애하게 될 것이다. 그러나 여기서는 다른 방침을 취할 것이다. 즉, 모델의 예측값을 탐색함으로써 모델을 이해하는데 초점을 맞출 것이다. 이 방법은 큰 장점을 가지고 있다. 모든 유형의 예측 모델은 예측값을 생성하므로 모든 유형의 예측 모델을 이해하기 위해 같은 기법을 사용할 수 있다. 또한, 모델에서 포착하지 못한 것(데이터에서 예측값을 뺀 나머지 부분으로 흔히 잔차라고 함)을 확인할 때 유용하다. 잔차는 눈에 띄는 패턴을 제거하여 감지하기 어려운 나머지 추세를 분석할 수 있도록 해주므로 중요하다고 할 수 있다.

18.3.1 예측값

모델의 예측값을 시각화하기 위해 데이터가 존재하는 영역을 포함하는 균일한 간격의 그리드를 생성한다. 가장 쉬운 방법은 modelr::data_grid()를 사용하는 것이다. 첫 번째 인수는 데이터프레임이며, 그다음 인수에 대해서는 고유한 값에 대한 모든 조합을 생성한다.

grid <- sim1 %>%
  data_grid(x)
grid
## # A tibble: 10 x 1
##        x
##    <int>
##  1     1
##  2     2
##  3     3
##  4     4
##  5     5
##  6     6
##  7     7
##  8     8
##  9     9
## 10    10
  • 다음은 예측값을 추가한다. 데이터프레임과 모델을 인수로 갖는 modelr::add_predictions()을 사용한다. 이 함수는 모델의 예측값을 데이터 프레임의 새로운 열로 추가한다.
grid <- grid %>%
  add_predictions(sim1_mod)
grid
## # A tibble: 10 x 2
##        x  pred
##    <int> <dbl>
##  1     1  6.27
##  2     2  8.32
##  3     3 10.4 
##  4     4 12.4 
##  5     5 14.5 
##  6     6 16.5 
##  7     7 18.6 
##  8     8 20.6 
##  9     9 22.7 
## 10    10 24.7

(또한, 이 함수를 사용하여 원 데이터셋에 예측값을 추가할 수 있다.) 다음으로 예측값을 플롯에 나타낸다. geom_abline()을 사용하는 것과 비교하였을 때 다음의 추가 진행 작업에 어ᄄᅠᆫ 차이점이 있는지 궁금할 것이다. 이 접근법의 장점은 가장 간단한 모델부터 복잡한 모델까지 R의 모든 모델에 동작한다는 것이다. 오직 시각화 기술에 의해서만 제한된다. 더 복잡한 모델 유형을 시각화 하는 방법에 대해 알고 싶다면 http://vita.had.co.nz/papers/model-vis.html 을 참고하자.

ggplot(sim1, aes(x)) +
  geom_point(aes(y = y)) +
  geom_line(
    aes(y = pred),
    data = grid,
    colour = "red",
    size = 1
  )

18.3.2 잔차

예측값의 다른 측면에는 잔차가 있다. 예측값은 모델이 포착한 패턴을 나타내고, 잔차는 모델이 놓친 것을 알려준다. 잔차는 이전에 계산한 예측값과 관측값 사이의 차이값이다.
add_predictions()와 유사한 기능을 하는 add_residuals()을 사용해서 데이터에 잔차를 추가한다. 주의: 만들어진 그리드가 아니라 원래의 데이터셋을 사용한다. 이는 잔차를 계산하는데 실제 y값이 필요하기 때문이다.

sim1 <- sim1 %>%
  add_residuals(sim1_mod)
sim1
## # A tibble: 30 x 3
##        x     y    resid
##    <int> <dbl>    <dbl>
##  1     1  4.20 -2.07   
##  2     1  7.51  1.24   
##  3     1  2.13 -4.15   
##  4     2  8.99  0.665  
##  5     2 10.2   1.92   
##  6     2 11.3   2.97   
##  7     3  7.36 -3.02   
##  8     3 10.5   0.130  
##  9     3 10.5   0.136  
## 10     4 12.4   0.00763
## # ... with 20 more rows

잔차가 모델에 대해 전달해주는 것을 이해하기 위한 몇 가지 방법이 있다. 한 가지는 단순히 잔차의 분포를 이해할 수 있도록 하는 빈도 다각형을 그리는 것이다.

ggplot(sim1, aes(resid)) +
  geom_freqpoly(binwidth = 0.5)

  • 이는 모델의 품질을 보정하는 데 도움이 되며 다음과 같이 ‘관측값과 예측값이 얼마나 멀리 떨어져 있는가?’ 라는 질문에 답할 수 있다. 또한 잔차의 평균값은 항상 0이다.
    본래의 예측 변수 대신 잔차를 사용하여 플롯을 다시 그리고 싶을 수도 있다. 다음 장에서 많은 부분을 확인하게 될 것이다.
ggplot(sim1, aes(x, resid)) +
  geom_ref_line(h = 0) +
  geom_point()

  • 위 플롯은 랜덤한 잔차를 나타내며 이는 모델이 데이터셋에서 패턴을 잘 포착하였다는 것을 의미한다.

18.4 수식과 모델 모음

이전에 facet_wrap()과 facet_grid()를 사용할 때 수식을 본 적이 있다. R에서 수식은 ‘특수한 동작’을 하도록 하는 일반적인 방법을 제공한다. 변수의 값을 즉시 평가하기보다는 함수에 의해 해석될 수 있도록 수식을 보유한다. R의 모델링 함수 대부분은 수식에서 함수로의 표준 전환을 사용한다. 이미 간단한 전환 한 가지는 보았다. y ~ x는 y = a_1 * a_2 * x로 전환된다. R에서 실제로 전환하는 것을 보고 싶다면 model_matrix() 함수를 이용할 수 있다. 데이터프레임과 수식을 이용하여 모델 방정식을 정의하는 티블을 반환한다. 출력값의 각열은 모델 하나의 계수와 연결되며 함수는 항상 y = a_1 * out1 + a_2 * out_2 형식이다. 가장 간단한 케이스인 y ~x1을 통해 흥미로운 사실을 알아본다.

df <- tribble(
  ~y, ~x1, ~x2,
  4,2,5,
  5,1,6
)
model_matrix(df, y ~ x1)
## # A tibble: 2 x 2
##   `(Intercept)`    x1
##           <dbl> <dbl>
## 1             1     2
## 2             1     1

R이 모델에 y절편을 추가하는 방법은 값이 1인 열을 추가하는 것이다. 기본적으로 R은 항상 1인 열을 추가한다. 만약 이를 원하지 않는다면 –1을 사용하여 제거해야 한다.

model_matrix(df, y ~ x1 - 1)
## # A tibble: 2 x 1
##      x1
##   <dbl>
## 1     2
## 2     1

모델에 변수를 추가할 때마다 모델 매트릭스는 증가한다.

model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 x 3
##   `(Intercept)`    x1    x2
##           <dbl> <dbl> <dbl>
## 1             1     2     5
## 2             1     1     6

떄때로 이 수식의 표기법은 ‘월킨슨-로저스 표기법’이라고 불리며, 원킨슨과 로저스의 Symbolic Description of Factorial Models for Analysis of Variance에 처음 설명되었다. 모델링 대수학의 모든 세부사항을 이해하고 싶다면 원 논문을 자세히 읽어볼 가치가 있다. 다음 절에서는 이러한 수식 표기법이 범주형 변수, 상호작용 및 변형에 적용되는 방식을 설명한다.
(논문링크 : https://www.jstor.org/stable/2346786?seq=1#page_scam_tab_contents)

18.4.1 범주형 변수

수식에서 함수를 생성하는 것은 예측 변수가 연속형일 때는 간단하지만, 예측 변수가 범주형일 때는 조금 복잡해진다. y ~ sex(sex는 남성 혹은 여성을 나타내는 변수)와 같은 함수가 있다고 가정해보자. sex는 숫자가 아니어서 곱할 수 없으므로 y = x_0 + x_1 * sex와 같은 수식으로 변환하는 것은 맞지 않다. 대신 R은 수식을 y = x_0 + x_1 * sex_male로 변환하며 여기서 sex_male은 sex가 남성이면 1, 그렇지 않으면 0을 의미한다.

df <- tribble(
  ~ sex, ~ response,
  "male", 1,
  "female", 2,
  "male", 1
)
model_matrix(df, response ~ sex)
## # A tibble: 3 x 2
##   `(Intercept)` sexmale
##           <dbl>   <dbl>
## 1             1       1
## 2             1       0
## 3             1       1

왜 R이 sexfemale 열은 만들지 않는지 궁금할 것이다. 그렇게 하면 다른 열(즉, sexfemale = 1 – sexfemale)을 기반으로 완전히 예측할 수 있는 열을 만들게 되기 때문이다. 안타깝게도 이것이 왜 문제인지에 대한 정확한 세부 내용은 이 책의 범위를 벗어난다. 간략하게 말하자면 너무 유연한 모델 모음이 생성되어 데이터와 유사한 무한히 많은 모델을 생성하게 된다. 그렇지만 다행히 예측값을 시각화하는 데 초점을 맞추면 정확한 매개 변수화에 대해서 걱정할 필요가 없다. 구체적인 데이터와 모델에 대해 살펴보자. 다음과 같이 modelr의 sim2 데이터셋이 존재한다.

ggplot(sim2) +
  geom_point(aes(x, y))

이 데이터에 대해 모델을 적합하고 예측값을 생성할 수 있다.

mod2 <- lm(y ~ x, data = sim2)
grid <- sim2 %>%
  data_grid(x) %>%
  add_predictions(mod2)
grid
## # A tibble: 4 x 2
##   x      pred
##   <chr> <dbl>
## 1 a      1.15
## 2 b      8.12
## 3 c      6.13
## 4 d      1.91

사실상 범주형 변수 x를 포함한 모델은 각 범주의 평균값을 예측한다. (그 이유는 각 범주의 평균값이 평균제곱근 편차를 최소화하기 떄문이다.) 이는 원 데이터 위에 예측값을 겹쳐서 그려보면 쉽게 확인할 수 있다.

ggplot(sim2, aes(x)) +
  geom_point(aes(y = y)) +
  geom_point(
    data = grid,
    aes(y = pred),
    color = "red",
    size=4 
  )

관측하지 못한 값에 대해서는 예측할 수 없다. 실수로 이 작업을 수행하게 되면 다음과 같은 에러 메시지가 나타날 것이다.

tibble(x = "e") %>%
  add_predictions(mod2)
  
Error in eval(predvars, data, env) : 객체 'x1'를 찾을 수 없습니다

18.4.2 연속형과 범주형 변수의 상호작용

연속형 변수와 범주형 변수를 결합하면 어ᄄᅠᇂ게 되는가? sim3는 범주형 예측 변수와 연속형 예측 변수를 포함한 데이터셋이다. 이 변수들은 간단한 플롯으로 시각화할 수 있다.

ggplot(sim3, aes(x1, y)) +
  geom_point(aes(color = x2))

  • 이 데이터에 적용할 수 있는 두 가지 모델이 있다.
mod1 <- lm(y ~ x1 + x2, data = sim3)
mod2 <- lm(y ~ x1 * x2, data = sim3)

+를 사용하여 변수를 추가하면 모델은 다른 모든 변수와 독립적인 각 효과를 추정한다. *을 상용하면 상호작용이라 불리는 항을 적합할 수 있다. 예를 들어 y ~ x1 * x2는 y = a_0 + a_1 * a1 + a_2 * a2 + a_12 * a1 * a2로 변환된다. *를 사용하면 상호작용과 개별 구성요소 모두 모델에 포함된다.

이 모델을 시각화하려면 새로운 두가지 수법이 필요하다.

  • 두 개의 예측 변수에 대해 data_grid()를 적용해야 한다. 이 함수는 x1과 x2의 고유한 값에 대한 모든 조합을 생성한다.

  • 두 모델로 동시에 예측값을 생성하기 위해서 행으로 각 예측값을 추가하는 gather_predictions()를 사용할 수 있다. gather_predictions()를 보완하는 함수는 각각의 예측값을 새로운 열로 추가하는 spread_predictions()이다.

위 함수를 함꼐 사용하면 다음과 같은 결과를 얻을 수 있다.

grid <- sim3 %>%
  data_grid(x1, x2) %>%
  gather_predictions(mod1, mod2)
grid
## # A tibble: 80 x 4
##    model    x1 x2     pred
##    <chr> <int> <fct> <dbl>
##  1 mod1      1 a      1.67
##  2 mod1      1 b      4.56
##  3 mod1      1 c      6.48
##  4 mod1      1 d      4.03
##  5 mod1      2 a      1.48
##  6 mod1      2 b      4.37
##  7 mod1      2 c      6.28
##  8 mod1      2 d      3.84
##  9 mod1      3 a      1.28
## 10 mod1      3 b      4.17
## # ... with 70 more rows

면분할을 사용해서 두 모델의 결과를 하나의 플롯에 시각화할 수 있다.

ggplot(sim3, aes(x1, y, color = x2)) +
  geom_point() +
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model)

+를 사용한 모델은 각 라인의 기울기는 같지만 y 절편값은 서로 다르다. *를 사용한 모델은 기울기와 y 절편값이 모두 다르다. 어떤 모델이 이 데이터에 잘 맞을까? 이는 잔차를 통해 확인할 수 있다. 여기서는 각 그룹 내의 패턴을 쉽게 확인하기 위해 모델과 x2로 면분할하였다.

sim3 <- sim3 %>%
  gather_residuals(mod1, mod2)

ggplot(sim3, aes(x1, resid, color = x2)) +
  geom_point() +
  facet_grid(model ~ x2)

mod2의 잔차에서는 확실한 패턴이 거의 보이지 않는다. mod1의 잔차는 분명하게 모델 b에서 어떤 패턴을 놓쳤다는 것을 나타내주며 b보다는 덜하지만 c와 d에서도 여전히 패턴이 존재하는 것을 알 수 있다. mod1 또는 mod2 둘 중 어떤 모델이 더 좋은지 알 수 있는 정확한 방법이 있는지 궁금할 것이다. 다방면의 수학적 배경지식이 필요한 방법이 있긴 하지만 그 방법에 대해서는 고려하지 않을 것이다. 여기서는 모델이 관심 있는 패턴을 포착하였는지와 관련된 질적 평가에 집중할 것이다.

18.4.3 연속형 변수의 상호작용

두 개의 연속 변수로 이루어진 같은 모델을 살표보자. 초반에는 이전 예제와 거의 유사하게 진행된다.

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
## # A tibble: 50 x 4
##    model    x1    x2   pred
##    <chr> <dbl> <dbl>  <dbl>
##  1 mod1   -1    -1    0.996
##  2 mod1   -1    -0.5 -0.395
##  3 mod1   -1     0   -1.79 
##  4 mod1   -1     0.5 -3.18 
##  5 mod1   -1     1   -4.57 
##  6 mod1   -0.5  -1    1.91 
##  7 mod1   -0.5  -0.5  0.516
##  8 mod1   -0.5   0   -0.875
##  9 mod1   -0.5   0.5 -2.27 
## 10 mod1   -0.5   1   -3.66 
## # ... with 40 more rows

data_grid() 안에 seq_range()를 사용한다. x의 고유한 모든 값을 사용하는 대신 최솟값과 최대값 사이를 균일하게 간격을 나눈 5개의 값을 사용할 것이다. 여기서는 이 부분이 매우 중요하지는 않지만, 일반적으로는 유용한 기법이다. seq_range() 함수에는 3개의 유용한 인수가 있다.

  • pretty = TRUE는 ‘보기 좋은’ 시퀀스, 즉 사람의 눈에 보기 좋은 시퀀스를 생성한다. 이는 테이블로 결과물을 생성하려는 경우에 유용하다.
seq_range(c(0.0123, 0.923423), n = 5)
## [1] 0.0123000 0.2400808 0.4678615 0.6956423 0.9234230
seq_range(c(0.0123, 0.923423), n = 5, pretty = TRUE)
## [1] 0.0 0.2 0.4 0.6 0.8 1.0
  • trim = 0.1은 꼬리 값의 10%를 제거한다. 변수가 꼬리가 긴 분포를 가지고 있으며, 중심 근처의 값을 생성하고자 하는 경우에 유용하다.
x1 <- rcauchy(100)
(rcauchy = 위치 모수 위치 및 척도 모수 척도를 사용하여 코시 분포에 대한 밀도, 분포 함수, 분위 함수 및 임의 생성.) 
seq_range(x1, n = 5)
## [1] -35.048371  -5.946321  23.155729  52.257779  81.359829
seq_range(x1, n = 5, trim = 0.10)
## [1] -8.7953166 -5.5983150 -2.4013134  0.7956882  3.9926899
seq_range(x1, n = 5, trim = 0.25)
## [1] -3.0086877 -1.7109689 -0.4132500  0.8844688  2.1821876
seq_range(x1, n = 5, trim = 0.50)
## [1] -1.2103737 -0.6833697 -0.1563657  0.3706383  0.8976423
  • expand = 0.1은 어떤 의미에서는 trim()의 반대이다. 즉, 범위를 10% 확장한다.
x2 <- c(0, 1)
seq_range(x2, n = 5)
## [1] 0.00 0.25 0.50 0.75 1.00
seq_range(x2, n = 5, expand = 0.10)
## [1] -0.050  0.225  0.500  0.775  1.050
seq_range(x2, n = 5, expand = 0.25)
## [1] -0.1250  0.1875  0.5000  0.8125  1.1250
seq_range(x2, n = 5, expand = 0.50)
## [1] -0.250  0.125  0.500  0.875  1.250

다음으로 모델을 시각화해보자. 2개의 연속형 예측 변수를 가지고 있으므로 3D표면과 같은 모델을 상상해볼 수 있다. geom_tile()을 사용하여 이를 나타낼 수 있다.

ggplot(grid, aes(x1, x2)) +
  geom_tile(aes(fill = pred)) +
  facet_wrap(~ model)

위 플롯은 두 모델이 매우 다르지 않다는 것을 나타낸다. 하지만 어는 정도는 우리의 착각일 수 있다. 우리의 눈과 뇌는 색채의 농도를 정확하게 비교하는 데 그리 좋지 않다. 위에서부터 표면을 살펴보는 대신에 한 변수를 여러 조각으로 나누어 살펴볼 수 있다.

ggplot(grid, aes(x1, pred, color = x2, group = x2)) +
  geom_line() +
  facet_wrap(~ model)

ggplot(grid, aes(x2, pred, color = x1, group = x1)) +
  geom_line() +
  facet_wrap(~ model)

위 플롯은 두 개의 연속 변수 사이의 상호작용이 기본적으로 범주형과 연속형 변수의 상호작용과 같은 방식으로 동작함을 보여준다. 또한 이 플롯에서의 상호작용은 고정된 오프셋이 없다는 것을 나타낸다. 즉 y를 예측하기 위해서는 x1과 x2를 동시에 고려해야 한다. 단 두 개의 연속 변수만으로 좋은 시각화는 어렵다는 것을 알 수 있다. 하지만 그것은 합리적이다. 세 개 또는 그 이상의 변수들이 동시에 상호작용하는 방식에 대해서 이해하기가 쉬울 거라 예상하지는 않을 것이다. 다시 말하자면 모델을 탐색에 사용하면 시간을 조금 절약할 수 있고, 모델의 성능은 점점 개선할 수 있다. 모델은 완벽해야 할 필요가 없으며 단지 데이터를 더 잘 나타내는 데 도움을 주면 된다. ‘mod2가 mod1보다 더 좋은지’에 대해 확인해보기 위해 잔차를 탐색해보았다. 그 결과, 차이가 발생하기는 하지만 미묘하다. 연습문제에서 이를 확인해볼 수 있을 것이다.

18.4.4 변환

모델의 수식 내에서 변환할 수도 있다. 예를 들어 log(y) ~ sqrt(x1) + x2는 y = a_1 + a_2 * x1 * sqrt(x) + a_3 * x2로 변환된다. 변환에 +, , ^ 또는 -가 포함되어 있으면 R에서는 모델의 열거한 부분으로 처리하지 않게 하도록 I()로 묶어야 한다. 예를 들어 y ~ x + I(x ^ 2)은 y = a_1 + a_2 x +a_3 * x^2으로 변환된다. 만약 I()를 쓰지 않고 y ~ x ^ 2 + x로 명시한다면 R은 y ~ x * x + x로 계산할 것이다. x * x는 x자체의 상호작용을 의미한다. R은 자동으로 중복 변수를 제외하므로 x + x는 x가 되며 y ~ x ^ 2 + x는 y = a_1 + a_2 * x 함수가 된다. 이는 분석가가 의도한 것이 아니다. 다시 말하자면 모델이 변환되는 방식이 햇갈린다면 model_matrix()를 사용하여 lm()이 모델을 적합하는 수식이 무엇인지 정확하게 확인할 수 있다.

df <- tribble(
  ~y, ~x,
  1, 1,
  2, 2,
  3, 3
)

model_matrix(df, y ~ x^2 + x)
## # A tibble: 3 x 2
##   `(Intercept)`     x
##           <dbl> <dbl>
## 1             1     1
## 2             1     2
## 3             1     3
model_matrix(df, y ~ I(x^2) + x)
## # A tibble: 3 x 3
##   `(Intercept)` `I(x^2)`     x
##           <dbl>    <dbl> <dbl>
## 1             1        1     1
## 2             1        4     2
## 3             1        9     3

변환은 비선형 함수를 근사하는 데 유용하게 사용할 수 있다. 만약 미적분 수업을 들었다면 ‘어떤 평활 함수도 다향식의 무한한 합으로 근사시킬 수 있다’라는 테일러 정리에 대해 들어봤을 것이다. 즉, y = a_1 + a_2 * x + a_3 * x^2 + a_4 * x ^ 3과 같은 방정식을 적합하여 선형 함수를 임의의 평활함수에 가갑게 만들 수 있다는 것을 의미한다. 손으로 시퀀스를 타이핑하는 것은 번거롭기 떄문에 R에서는 도우미 함수인 poly()를 제공한다.

model_matrix(df, y ~ poly(x, 2))
## # A tibble: 3 x 3
##   `(Intercept)` `poly(x, 2)1` `poly(x, 2)2`
##           <dbl>         <dbl>         <dbl>
## 1             1     -7.07e- 1         0.408
## 2             1     -7.85e-17        -0.816
## 3             1      7.07e- 1         0.408

그러나 poly()를 사용하는 데는 한 가지 주요한 문제가 있다. 데이터의 범위를 벗어나면 다향식은 급격하게 양의 무한대 또는 음의 무한대로 발산하게 된다. 한 가지 안전한 대안은 본연의 스플라인(매끄러운 곡선)인 splins::ns()를 사용하는 것이다.

library(splines)
model_matrix(df, y ~ ns(x, 2))
## # A tibble: 3 x 3
##   `(Intercept)` `ns(x, 2)1` `ns(x, 2)2`
##           <dbl>       <dbl>       <dbl>
## 1             1       0           0    
## 2             1       0.566      -0.211
## 3             1       0.344       0.771

비선형 함수를 근사할 때 어떻게 표현되는지 살표보자.

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")

ggplot(sim5, aes(x, y)) +
  geom_point() +
  geom_line(data = grid, color = "red") +
  facet_wrap(~ model)

데이터 범위를 벗어나 추정하는 것은 명백히 좋지 않다. 이는 다항식으로 함수를 근사하는 방법에서 좋지 않은 부분이지만, 모든 모델에서 발생할 수 있는 실질적인 문제이다. 데이터의 범위 밖에서 추정을 시작하면 변화하는 추정값이 맞는지에 대한 여부를 모델에서 알려줄 수 없다. 이론과 과학에 의존해야 한다.

18.5 결측값

결측값은 변수 간의 관계에 대해 어떤 정보도 전달할 수 없으므로 모델링 함수는 결측값이 포함된 행을 삭제한다. R의 기본 동작은 자동으로 제거하는 것이지만 (전제조건에서 실행한) options(na.action = na.warn)은 경고문을 표시해준다.

df <- tribble(
  ~x, ~y,
  1, 2.2,
  2, NA,
  3, 3.5,
  4, 8.3,
  NA, 10
)

mod <- lm(y ~ x, data = df)
경고메시지(들): 
Dropping 2 rows with missing values

경고문을 표시하지 않으려면 na.action = na.exclude로 설정한다.

mod <- lm(y ~ x, data = df, na.action = na.exclude)

nobs()를 사용하면 얼마나 많은 관측값이 사용되었는지 정확하게 알 수 있다.

nobs(mod)
## [1] 3

18.6 다른 모델 모음

이 장에서는 y = a_1 * x1 +a_2 * x2 + … + a_n * xn 형식으로 된 관계를 가정하는 선형 모델 종류에 대해서만 다룬다. 또한 이 책에서는 다루지 않았지만 선형 모형은 잔차가 정규분포를 따른다고 가정한다. 여러 흥미로운 방법으로 선형 모델을 확장한 모델의 종류는 다양하며 일부 모델을 다음에서 소개한다.

  • 일반화 선형 모형: 예, stats::glm()
    선형 모형에서는 반응 변수가 연속형이며 오차가 정규분포를 따른다고 가정한다. 일반화 선형 모형은 연속형이 아닌 반응 변수(예, 이항 데이터 또는 빈도수 데이터)를 포함할 수 있도록 선형 모형을 확장한다, 일반화 선형 모형은 우도의 통계적 아이디어를 기반으로 거리 메트릭을 정의하여 동작한다.

  • 일반화 가법 모형: 예, mgcv::gam()
    일반화 선형 모형을 확장하여 임의의 평활 함수를 포함하도록 한다. 즉, y = f(x) 형식의 방정식으로 변환되는 y ~ s(x) 형식의 수식을 작성한 후, gam()을 이용하여 그 함수가 무엇인지 추정한다. 이때 문제를 다루기 쉽게 만들기 위해 평활도를 제한조건으로 설정해야 한다.

  • 벌점 선형 모형: 예, glmnet::glmnet()
    복잡한 모델에 벌점을 부과하는 패널티 항을 (모수 백터와 원점 사이의 거리로 정의되는) 차이도에 추가한다. 이 방법은 같은 모집단의 새로운 데이터셋을 좀더 일반화할 수 있는 모델을 만드는 경향이 있다.

  • 로버스트 선형 모형: 예, MASS:rlm()
    매우 멀리 떨어져 있는 가중치가 적은 점과의 차이를 조정한다. 이 모형은 이상치가 없을 EO는 크게 효과가 없지만 이상값이 존재한다면 덜 민감하도록 만든다.

  • 트리 모형: rpart::rpart()
    선형 모형과 완전히 다른 방식으로 문제를 다룬다. 트리 모형은 데이터를 계속해서 작은 부분으로 분리하면 조각별 상수 모델을 적합한다. 트리 모형 그 자체로는 매우 효과적이지는 않지만 랜덤 포레스트(예, randomForest::randomForest())나 그레이디언트 부스팅(예, xgboost::xgboost)과 같은 모델로 결합하여 사용하면 매우 강력해 진다.

이 모델들은 프로그래밍 관점에서 모두 유사하게 동작한다. 일단 선형 모형을 완전히 익히고 나면 다른 모델 종류의 메커니즘을 쉽게 이해할 수 있다. 숙련된 모델러가 된다는 것은 몇 가지 괜찮은 원칙의 결합과 기법들이 담긴 큰 도구 상자를 가진다는 것이다. 이제 일반적인 도구와 한 가지의 유용한 모델 종류를 배웠으므로 다른 자료를 통해 더 많은 내용을 배울 수 있을 것이다.