랜덤 포레스트로 시에라리온의 물 가용성 예측 🚰

이것은 screencasts 패키지를 사용하는 방법을 보여 주는 최신 시리즈tidymodels로서 시작부터 많은 하이퍼파라미터로 더 복잡한 모델을 조정하는 것까지 포함됩니다. 오늘의 스크린캐스트는 수원에 대한 이번 주 #TidyTuesday dataset와 함께 무작위 숲 모델을 훈련하고 평가하는 방법을 안내합니다. 🚰



다음은 비디오 대신 또는 비디오에 추가하여 읽기를 선호하는 사람들을 위해 비디오에서 사용한 코드입니다.

데이터 탐색

우리의 모델링 목표는 방문 중에 관찰된 수원의 특성을 기반으로 water source에 실제로 사용 가능한 물이 있는지 예측하는 것입니다. 데이터를 읽는 것부터 시작하겠습니다.

library(tidyverse)
water_raw <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-05-04/water.csv")

이 모델 분석을 시에라리온의 수원과 물 가용성에 대해 "y"또는 "n"으로 분류된 수원으로만 제한하겠습니다. 이러한 수원은 시에라리온 전역에 어떻게 분포되어 있습니까?

water_raw %>%
  filter(
    country_name == "Sierra Leone",
    lat_deg > 0, lat_deg < 15, lon_deg < 0,
    status_id %in% c("y", "n")
  ) %>%
  ggplot(aes(lon_deg, lat_deg, color = status_id)) +
  geom_point(alpha = 0.1) +
  coord_fixed() +
  guides(color = guide_legend(override.aes = list(alpha = 1)))

앞으로 사용할 새 water 데이터 세트를 만들고 pay 변수를 처리해 보겠습니다.

water <- water_raw %>%
  filter(
    country_name == "Sierra Leone",
    lat_deg > 0, lat_deg < 15, lon_deg < 0,
    status_id %in% c("y", "n")
  ) %>%
  mutate(pay = case_when(
    str_detect(pay, "^No") ~ "no",
    str_detect(pay, "^Yes") ~ "yes",
    is.na(pay) ~ pay,
    TRUE ~ "it's complicated"
  )) %>%
  select(-country_name, -status, -report_date) %>%
  mutate_if(is.character, as.factor)

수원이 설치된 시기에 따라 물 가용성의 차이가 보이나요?

water %>%
  ggplot(aes(install_year, y = ..density.., fill = status_id)) +
  geom_histogram(position = "identity", alpha = 0.5) +
  labs(fill = "Water available?")

수원지의 납부현황별은 어떻게 되나요?

water %>%
  ggplot(aes(y = pay, fill = status_id)) +
  geom_bar(position = "fill") +
  labs(fill = "Water available?")

이는 데이터 자체의 문제일 수 있습니다. 방문 당시 물이 없었던 수원에 대한 지불에 대한 정보가 많지 않을 수 있습니다. 이것은 확실히 더 배울 가치가 있습니다!

모델 구축

"데이터 예산"을 설정하여 모델링을 시작하겠습니다.

library(tidymodels)

set.seed(123)
water_split <- initial_split(water, strata = status_id)
water_train <- training(water_split)
water_test <- testing(water_split)

set.seed(234)
water_folds <- vfold_cv(water_train, strata = status_id)
water_folds


## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 2
## splits id    
## <list> <chr> 
## 1 <split [36985/4110]> Fold01
## 2 <split [36985/4110]> Fold02
## 3 <split [36985/4110]> Fold03
## 4 <split [36985/4110]> Fold04
## 5 <split [36985/4110]> Fold05
## 6 <split [36986/4109]> Fold06
## 7 <split [36986/4109]> Fold07
## 8 <split [36986/4109]> Fold08
## 9 <split [36986/4109]> Fold09
## 10 <split [36986/4109]> Fold10

이 분석을 위해 usemodels 패키지를 사용하여 모델링 코드를 빠르게 설정했습니다.

usemodels::use_ranger(status_id ~ ., data = water_train)

이 코드는 내가 들어가서 편집할 수 있는 코드를 생성했습니다. 특히 모델링 코드의 기능 엔지니어링 부분에 추가해야 했습니다.

library(themis)
ranger_recipe <-
  recipe(formula = status_id ~ ., data = water_train) %>%
  update_role(row_id, new_role = "id") %>%
  step_unknown(all_nominal_predictors()) %>%
  step_other(all_nominal_predictors(), threshold = 0.03) %>%
  step_impute_linear(install_year) %>%
  step_downsample(status_id)

ranger_spec <-
  rand_forest(trees = 1000) %>%
  set_mode("classification") %>%
  set_engine("ranger")

ranger_workflow <-
  workflow() %>%
  add_recipe(ranger_recipe) %>%
  add_model(ranger_spec)

doParallel::registerDoParallel()
set.seed(74403)
ranger_rs <-
  fit_resamples(ranger_workflow,
    resamples = water_folds,
    control = control_resamples(save_pred = TRUE)
  )

모델링 완료! 랜덤 포레스트는 일반적으로 충분한 트리를 제공하면 꽤 잘 작동하기 때문에 랜덤 포레스트를 조정하지 않기로 선택했습니다.

결과 살펴보기

랜덤 포레스트 모델은 어떻게 수행되었습니까?

collect_metrics(ranger_rs)


## # A tibble: 2 x 6
## .metric .estimator mean n std_err .config             
## <chr> <chr> <dbl> <int> <dbl> <chr>               
## 1 accuracy binary 0.893 10 0.00179 Preprocessor1_Model1
## 2 roc_auc binary 0.951 10 0.00145 Preprocessor1_Model1

10개의 교차 검증 폴드 세트에 대한 ROC 곡선을 볼 수 있습니다.

collect_predictions(ranger_rs) %>%
  group_by(id) %>%
  roc_curve(status_id, .pred_n) %>%
  autoplot()

conf_mat_resampled()를 사용하여 재표본에서 혼동 행렬을 만들 수도 있습니다.

conf_mat_resampled(ranger_rs, tidy = FALSE) %>%
  autoplot()

이 결과는 (pay 기능에 대한 주의 사항과 함께) 꽤 괜찮아 보입니다. 마지막으로 한 번에 전체 교육 세트에 맞추고 테스트 세트에서 한 번 평가해 보겠습니다. 테스트 세트를 사용한 것은 이번이 처음입니다.

final_fitted <- last_fit(ranger_workflow, water_split)
collect_metrics(final_fitted) ## metrics on the *testing* set


## # A tibble: 2 x 4
## .metric .estimator .estimate .config             
## <chr> <chr> <dbl> <chr>               
## 1 accuracy binary 0.892 Preprocessor1_Model1
## 2 roc_auc binary 0.951 Preprocessor1_Model1

테스트 세트에 대한 예측을 수집하고 ROC 곡선 또는 여기에 표시된 대로 혼동 행렬을 만들 수 있습니다.

collect_predictions(final_fitted) %>%
  conf_mat(status_id, .pred_class) %>%
  autoplot()

변수 중요도는 어떻습니까? vip 패키지를 사용하여 이 랜덤 포레스트에 대한 컴퓨터 변수 중요성을 알아봅시다. importance = "permutation"와 수동으로 사전 처리한 데이터 세트imp_data를 사용하여 모델을 다시 맞춰야 합니다.

library(vip)

imp_data <- ranger_recipe %>%
  prep() %>%
  bake(new_data = NULL) %>%
  select(-row_id)

ranger_spec %>%
  set_engine("ranger", importance = "permutation") %>%
  fit(status_id ~ ., data = imp_data) %>%
  vip(geom = "point")

pay가 가장 중요하다는 것은 놀라운 일이 아닙니다! 다음으로 가장 중요한 변수는 수원에서 사용되는 기술과 이를 설치한 사람입니다. 전처리된 데이터를 사용하여 플롯을 하나 더 만들어 데이터가 어떻게 분포되어 있는지 살펴보겠습니다.

imp_data %>%
  select(status_id, pay, water_tech, installer) %>%
  pivot_longer(pay:installer, names_to = "feature", values_to = "value") %>%
  ggplot(aes(y = value, fill = status_id)) +
  geom_bar(position = "fill") +
  facet_grid(rows = vars(feature), scales = "free_y", space = "free_y") +
  theme(legend.position = "top") +
  scale_fill_brewer(type = "qual", palette = 7) +
  scale_x_continuous(expand = expansion(mult = c(0, .01)), labels = scales::percent) +
  labs(
    x = "% of water sources", y = NULL, fill = "Water available?",
    title = "Water availability by source characteristic in Sierra Leone",
    subtitle = "Water sources with no payment information are likely to have no water available"
  )