[Tidy Modeling with R] 11. Model Comparison (모델 비교)

11 Comparing Models with Resampling | Tidy Modeling with R

 

 

 

만약 모델을 두개 이상 만들었다면, 이제는 최고의 모델을 선택하기 위해서 모델의 성능을 서로 비교해야 합니다.

 

모델 비교(Model Comparison)는 크게 모델 내 비교(within-model comparison)와 모델 간 비교(between-model comparison)의 2가지로 구분할 수 있으며 전자의 경우에는 서로 다른 변수선택 또는 전처리를 동일한 모델에 사용하여 나온 결과를 비교하는 것을 뜻하며, 후자의 경우에는 선형회귀와 랜덤포레스트와 같은 서로 다른 모델을 비교하는 경우를 뜻합니다.

 

두 가지 경우 모두 비교를 하기 위해서 RMSE, AUC와 같은 리샘플링 요약 통계량(Resampled Summary statistics)을 이용해야 합니다.

 

 

1. Creating Mulitple models with workflow sets

첫번째 코드는 리샘플링 데이터를 준비하는 코드로 이전 포스팅에서 다룬 것 처럼 8:2로 Train-Test 분할 후, 분할된 Train data에 대해 10-Fold CV를 적용하여 10개의 Analysis-Assessment set을 구성하였습니다. 

library(tidymodels)
library(tidyverse)

ames <- modeldata::ames |> mutate(Sale_Price = log10(Sale_Price))
ames_split <- ames |> rsample::initial_split(strata = Sale_Price, prop = 0.8)
ames_train <- ames_split |> training()
ames_test <- ames_split |> testing()
ames_vfolds <- ames_train |> vfold_cv(v=10)

Dataset 구성

 

 

두번째 코드는 기본이 되는 레시피 전처리기(Preprocessing)에 하나씩 다른 요소를 추가하여 3가지의 전처리기를 생성하였고, PCA를 적용한 전처리기를 합해 총 4가지의 레시피를 생성하였습니다.

### recipe
# 기본 전처리 : 반응변수, 설명변수 세팅 / Gr_Liv_Area 로그 변환 / Neighborhood 빈도 적은 변수 처리 / 명목변수 더미화
basic_recipe <- 
  recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area + 
           Year_Built + Bldg_Type + Latitude + Longitude, 
         data = ames_train) |> 
  step_log(Gr_Liv_Area, base=10) |> 
  step_other(Neighborhood, threshold = 0.01) |> 
  step_dummy(all_nominal_predictors())

# 기본 전처리에 상호작용항 추가 
interaction_recipe <- basic_recipe |> 
  step_interact(~ Gr_Liv_Area : starts_with("Bldg_Type_"))

# 상호작용 전처리에 경도와 위도의 ns spline 추가
spline_recipe <- interaction_recipe |> 
  step_ns(Latitude, Longitude, deg_free = 20)

# PCA Recipe 
pca_recipe <- 
  recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area + 
           Year_Built + Bldg_Type + Latitude + Longitude, 
         data = ames_train) |> 
  step_log(Gr_Liv_Area, base=10) |> 
  step_other(Neighborhood, threshold = 0.01) |> 
  step_pca(all_numeric_predictors(), num_comp = 5) |> 
  step_dummy(all_nominal_predictors())

# 전처리 저장
preproc <- list(basic = basic_recipe,
                interact = interaction_recipe,
                splines = spline_recipe,
                pca = pca_recipe)

 

 

이제 Recipe 객체를 사용한 전처리기를 정의하였으니 workflow_set 함수를 사용하여 워크플로의 집합인 workflow set 객체를 생성하여 봅시다.

lm_models <- workflow_set(preproc = preproc,
                          models = list(lm = linear_reg() |> set_mode("regression")))

워크플로우 셋에는 4개의 워크플로가 있습니다.

 

 

10장에서는 워크플로에 fit_resamples 함수를 사용해 리샘플링 학습과 평가를 할 수 있다는 사실 기억하시나요?

4개의 워크플로에 각각 fit_resamples를 적용하여 결과를 얻을 수 있지만, 코드 중복이 생기고 번거로울 수 있는데요.
(지금은 4개지만,  4개가 아니라 10개, 100개라면... 하나하나 치고 있기 너무 힘들죠)

따라서, purrr::map 함수처럼 워크플로 객체마다 동일한 함수를 실행하는 workflow_map 함수를 사용해서 일괄적으로 fit_resamples 함수를 적용해 보겠습니다.

# workflow_map(fn="fit_resamples")
# >> 각각의 workflow에 fit_resamples를 적용
keep_pred <- control_resamples(save_pred = T, save_workflow = T, verbose = T)
lm_models <- lm_models |> 
  workflow_map(fn = "fit_resamples", verbose = T, seed=312,
               resamples = ames_folds, control = keep_pred)

• [syntax] workflowsets::workflow_map(object, fn="tune_grid", verbose=F, seed, ...)
• option 열은 fn함수에 주어지는 리샘플링과 컨트롤 옵션을 포함하고 있습니다.
• result 열은 각 object에 fn함수를 적용한 결과를 포함하고 있습니다.

 

 

또한, workflow set에 대해 collect_metrics 함수나 collect_predictions 함수를 사용해 편리하게 성과지표와 예측값을 출력할 수 있습니다.

lm_models %>% collect_metrics(summarize = F)
lm_models %>% collect_metrics(summarize = T)
lm_models %>% collect_predictions()

collect_metrics(summarise=F)인 경우 개별 폴드에 대한 성과 지표 추정치를 제공한다.

 

collect_metrics(summarise=T)인 경우 리샘플링 요약 통계량을 제공한다.

 

collect_predictions인 경우 개별 워크플로의 예측치를 제공한다.

 

 

이제 랜덤포레스트 모델을 함께 비교하기 위해 랜덤포레스트 워크플로를 생성하고 리샘플링 표본에 대해 학습 및 평가를 진행해 봅시다.

### Random Forest Model
rf_model <- rand_forest(trees=1000, engine = "ranger", mode = "regression")

rf_workflow <- workflow() %>% 
  add_variables(outcomes = "Sale_Price",
                predictors = c("Neighborhood", "Gr_Liv_Area", "Year_Built",
                               "Bldg_Type", "Latitude", "Longitude")) |> 
  add_model(rf_model)
  
rf_res <- rf_workflow |> fit_resamples(resamples = ames_folds, control = keep_pred)

 

 

fit_resamples를 사용하면 성과지표와 예측값이(control_resamples 옵션 사용한 경우) 생성되는데 as_workflow_set 함수를 사용해 앞서 미리 적합한 선형회귀 결과들과 묶어서 사용할 수 있습니다.

또한 $R^2$ 성능으로 모델을 비교하기 위해서 autoplot 함수를 사용한 시각화를 진행하였으며 랜덤포레스트의 결과가 제일 우수한 것으로 판단할 수 있겠습니다. 반면에 PCA를 적용한 경우에는 기본 레시피보다 순위가 낮을 것을 확인할 수 있겠습니다.

### 모델 성능 시각화
total_model <- as_workflow_set(random_forest = rf_res) |> 
  bind_rows(lm_models)

autoplot(total_model, id = "workflow_set")

workflowsets::autoplot(total_model, metric="rsq") +
  ggrepel::geom_label_repel(aes(label=wflow_id),
                            nudge_x = 0.5, nudge_y = 0.005) +
  theme(legend.position = "none")

Random Forest의 성능이 제일 좋습니다.

 

참고로 시각화를 보면 basic_recipe 보다는 splines_recipe를 적용한 결과가 조금 더 나아보입니다. 비록 차이는 작아보이지만 통계적으로 유의할 수 있는데요. basic_recipe에서 추가되는 항이 $R^2$를 통계적으로 향상시키는지 검정할 수 있습니다.

 

 

 

 

2. Comparing Resampled performance statstics

리샘플링을 사용하여 모델 간 비교를 수행할 때, 리샘플링의 폴드별로 통계량(=성과지표)이 유사한 경향이 있습니다.

(즉, 특정 폴드에서 랜덤포레스트의 결과가 좋다면 회귀모델의 결과도 좋게 나올 가능성이 크며 반대로 또다른 특정 폴드에서 랜덤포레스트의 결과가 나쁘게 나온다면 회귀모델의 결과도 나쁘게 나올 가능성이 큽니다. 이를 resample-to-resample component of variation 라고 합니다)

이를 살펴보기 위해서 전체 모델의 리샘플링 결과가 있는 total_model를 시각화에 적절한 형태로 만들기 위해 아래와 같이 코드를 작성하였습니다.

rsq_result <- total_model |> 
  collect_metrics(summarize = F) |> 
  filter(.metric == "rsq")
  
  rsq_result |> group_by(wflow_id) |> 
  slice_head(n=2)

 

 

ggridges 패키지를 이용한 시각화에서 워크플로가 다른 경우에도 등락이 폴드별로 유사한 경향을 보이는 것을 알 수 있습니다.

따라서, 시각화에서 리샘플링 상관(Resample Correlation)이 존재한다고 볼 수 있겠습니다.

library(ggplot2)
library(ggridges)
theme_set(theme_gray(base_family = "AppleGothic"))
rsq_result |> 
  ggplot(aes(x=id, y=wflow_id, height=.estimate, group=wflow_id, fill=as.factor(wflow_id))) + 
  geom_ridgeline(stat = "identity", scale=5) + 
  theme(legend.position = "none") + 
  labs(x=NULL, y=NULL, 
       title = latex2exp::TeX("Workflow별 Fold에 따른 $R^2$ 비교"))

 

 

통계적 검정 방법인 cor.test를 통해 상관계수의 신뢰구간을 확인할 수 있는데 basic_recipe를 사용한 회귀모델과 랜덤포레스트 모델을 비교해 봅시다.

rsq_metric <- rsq_result %>% 
  select(wflow_id, .estimate, id) %>% 
  pivot_wider(names_from = wflow_id, values_from = .estimate)

cor.test(rsq_metric$random_forest, rsq_metric$basic_lm)

• 상관계수 추정값은 0.9440이며 95% 신뢰구간은 [0.7751, 0.9869]입니다.
• 신뢰구간이 0을 포함하지 않으므로 (or p-value가 0.05보다 작으므로) basic_recipe를 사용한 선형회귀모델과 랜덤포레스트 모델의 폴드별 리샘플링 통계량($R^2$)에는 상관관계가 존재한다고 할 수 있습니다.

 

리샘플링 상관성이 존재한다면 두 모델간의 리샘플링 통계량 차이에 대한 검정에 문제를 발생시킬 수 있으며 이는 곧 리샘플링 통계량을 통해 모델을 비교하는데 있어서 resample-to-resample effect를 무시한다면 잘못된 결과를 얻을 수 있다는 뜻입니다. 본문에서 언급하는 해결 방법으로는 모델을 비교하기 전 Relevant practical effect size를 설정하고 통계적으로 유의하더라도 통계량의 차이가 Relevant practical effect size보다 작다면 실용적으로 중요하지 않다고 판단하는 방법을 소개하고 있습니다. 

 

 

 

 

3. Simple Hypothesis Testing Methods

모델 비교를 위해서 유명한 ANOVA를 적용할 수 있는데 $(p+1)$개의 모델 간 비교를 위해서 아래와 같은 선형 통계 모델을 생각해봅시다.

$$y_{ij} = \beta_0 + \beta_1x_{i1} + \ldots + \beta_px_{ip} + \epsilon_{ij}$$

$y_{ij}$를 우리가 관심있는 리샘플링 통계량(=성과지표)으로, 각 $x_{ij}$를 다른 모델 또는 항을 나타내는 더미변수라고 생각해 봅시다.

이런 경우에는 회귀계수 $\beta$의 검정을 통해 어떤 모델이 다른 모델과 다른지 테스트할 수 있을 것 입니다.

 

만약 $x_1, x_2, \ldots, x_n$이 전부 0인 경우 basic 모델, $x_1=1$이고, $x_2, x_3, \ldots, x_n$이 0인 경우를 상호작용항을 나타내는 모델이라고 생각해봅시다. 상호작용항이 실제로 유의해서 리샘플링 통계량에 영향을 준다면 상호작용항을 나타내는 더미변수의 회귀계수($\beta_1$)은 통계적으로 유의하다는 결론을 얻을 것 입니다. 반대로 상호작용항이 실제로 유의하지 않아 통계량에 영향을 주는 것이라고 생각하기 힘들다고한다면 $\beta_1 = 0$이라고 볼 수 있습니다.

 

다만 간단하게 P-value를 통해 모델이나 특정 항의 유의성을 판단할 수 있지만 resample-to-resample effect를 고려해야합니다. 이런 효과를 추가하기 위해 리샘플링 그룹을 나타내는 항을 추가하거나 리샘플(Resample)을 랜덤효과(Random Effect)로 처리하는 방식을 고려해볼 수 있습니다. 특히 랜덤효과를 포함하는 ANOVA 모델을 적합하기 위해서는 선형혼합모델(Linear Mixed Model)이 있습니다.

특히 두 모델을 비교하는 방법은 리샘플링 간 효과를 Matching(짝짓기)을 통해 제거해 검정할 수 있으므로 동일 리샘플링으로 생성된 $R^2$의 차이가 0과 같은지 검정해서 알아볼 수 있습니다. (대응표본 T검정과 비슷한 맥락)

 

 

 

 

4. Bayesian Methods

3절에서 소개한 Anova 검정을 위한 모델 $y_{ij} = \beta_0 + \beta_1x_{i1} + \cdots + \beta_px_{ip} + \epsilon_{ij}$에서 잔차항 $\epsilon_{ij}$는 서로 독립이고 정규분포를 따른다고 가정하는데요. 이번에 소개할 베이지안 방법에서는 잔차항의 분산($\sigma^2$) 또는 회귀계수의 분포($\beta_j$)에 대한 사전 분포(Prior distribution)를 가정합니다. (참고 : https://moogie.tistory.com/131)

 

리샘플링을 적절하게 모델에 반영하는 베이지안 모델을 적용하기 위해 임의절편모델(random intercept model)에 대해 생각해볼 수 있는데요. 리샘플링이 회귀계수를 변화시키는 것이 아닌 단지 절편에만 영향을 준다는 가정을 하고 있습니다.

$$y_{ij} = (\beta_0 + b_i) + \beta_1*x_{i1} + \beta_2*x_{i2} + \beta_3*x_{i3} + e_{ij}$$

해당 모델에서는 random effect인 $b_i$ 대한 사전 분포를 지정해야 합니다. 본문에는 $b_i \sim t(1)$이라고 가정했네요.

 

 

R에서 리샘플링 모델을 비교할 목적으로 tidyposterior 패키지를 사용할 수 있습니다.

주요 함수는 perf_mod 함수를 사용하며 workflow set에 대해서는 각 워크플로우에 대한 ANOVA 모델을 생성하며 모델 간 비교를 수행합니다. 또한 리샘플링 표본을 사용해 튜닝된 단일 모델 객체에 대해서는 모델 내 비교를 수행합니다. 어떤 객체가 주어지든 perf_mod 함수는 적절한 베이지안 모델을 결정하고 리샘플링 통계량을 사용해 적합합니다. 

library(tidyposterior)
library(rstanarm)
rsq_anova <- perf_mod(total_model, metric = "rsq", 
                      prior_intercept = student_t(df=1),
                      chains = 5,
                      iter = 10000,
                      seed = 1)

• 임의절편에 대한 분포를 제외하고는 default prior distribution을 사용했습니다
• 임의절편에 대한 사전분포는 T분포를 사용했습니다.

 

 

베이지안 통계에서 보통 관심이 있는 분포는 사후확률(Posterior Probability)이므로 아래와 같이 각 모델에 대한 $R^2$의 사후분포를 시각화합니다.

model_post <- rsq_anova %>% tidy()

model_post %>% 
  mutate(model = fct_inorder(model)) %>% 
  ggplot(mapping=aes(x=posterior)) + 
  geom_histogram(bins=50, color="white", fill="blue", alpha=0.4) + 
  facet_wrap(~model, ncol=1) + 
  labs(title = "Posterior Distribution of R-Square", y=NULL) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) + 
  scale_x_continuous(labels = scales::percent_format())