풍력 터빈 용량에 대한 결정 트리 조정 및 해석 🌬

이것은 screencasts 패키지를 사용하는 방법을 보여 주는 제 tidymodels 시리즈 중 첫 번째 모델링 단계부터 시작하여 더 복잡한 모델 튜닝에 이르는 최신 시리즈입니다. 오늘의 스크린캐스트에서는 이번 주 캐나다 풍력 터빈 #TidyTuesday dataset을 사용하여 의사 결정 트리 모델에서 조정, 적합 및 예측하는 방법을 안내합니다. 🇨🇦



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

데이터 탐색

우리의 모델링 목표는 캐나다based on other characteristics of the turbines from this week’s #TidyTuesday dataset의 풍력 터빈 용량을 예측하는 것입니다. Simon Couch가 이번 주how to use stacks for ensembling with this dataset,에 대해 설명했지만 여기서는 좀 더 간단한 접근 방식을 사용하겠습니다.

데이터를 읽는 것부터 시작하겠습니다.

library(tidyverse)

turbines <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-10-27/wind-turbine.csv")
turbines


## # A tibble: 6,698 x 15
## objectid province_territ… project_name total_project_c… turbine_identif…
## <dbl> <chr> <chr> <dbl> <chr>           
## 1 1 Alberta Optimist Wi… 0.9 OWE1            
## 2 2 Alberta Castle Rive… 44 CRW1            
## 3 3 Alberta Waterton Wi… 3.78 WWT1            
## 4 4 Alberta Waterton Wi… 3.78 WWT2            
## 5 5 Alberta Waterton Wi… 3.78 WWT3            
## 6 6 Alberta Waterton Wi… 3.78 WWT4            
## 7 7 Alberta Cowley North 19.5 CON1            
## 8 8 Alberta Cowley North 19.5 CON2            
## 9 9 Alberta Cowley North 19.5 CON3            
## 10 10 Alberta Cowley North 19.5 CON4            
## # … with 6,688 more rows, and 10 more variables:
## # turbine_number_in_project <chr>, turbine_rated_capacity_k_w <dbl>,
## # rotor_diameter_m <dbl>, hub_height_m <dbl>, manufacturer <chr>,
## # model <chr>, commissioning_date <chr>, latitude <dbl>, longitude <dbl>,
## # notes <chr>

약간의 데이터 정리 및 준비 작업을 수행해 보겠습니다.

turbines_df <- turbines %>%
  transmute(
    turbine_capacity = turbine_rated_capacity_k_w,
    rotor_diameter_m,
    hub_height_m,
    commissioning_date = parse_number(commissioning_date),
    province_territory = fct_lump_n(province_territory, 10),
    model = fct_lump_n(model, 10)
  ) %>%
  filter(!is.na(turbine_capacity)) %>%
  mutate_if(is.character, factor)

용량은 시운전 연도 또는 터빈 크기와 같은 다른 특성과 어떤 관련이 있습니까?

turbines_df %>%
  select(turbine_capacity:commissioning_date) %>%
  pivot_longer(rotor_diameter_m:commissioning_date) %>%
  ggplot(aes(turbine_capacity, value)) +
  geom_hex(bins = 15, alpha = 0.8) +
  geom_smooth(method = "lm") +
  facet_wrap(~name, scales = "free_y") +
  labs(y = NULL) +
  scale_fill_gradient(high = "cyan3")

이러한 관계는 modeling stacking Simon demonstrated이든 여기서 사용할 단일 모델이든 관계없이 모델링에 사용하려는 종류입니다.

모델 구축

tidymodels 메타패키지를 로드하고, 데이터를 교육 및 테스트 세트로 분할하고, 교차 검증 샘플을 생성하는 것으로 시작할 수 있습니다.

library(tidymodels)

set.seed(123)
wind_split <- initial_split(turbines_df, strata = turbine_capacity)
wind_train <- training(wind_split)
wind_test <- testing(wind_split)

set.seed(234)
wind_folds <- vfold_cv(wind_train, strata = turbine_capacity)
wind_folds


## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 2
## splits id    
## <list> <chr> 
## 1 <split [4.4K/488]> Fold01
## 2 <split [4.4K/487]> Fold02
## 3 <split [4.4K/486]> Fold03
## 4 <split [4.4K/486]> Fold04
## 5 <split [4.4K/486]> Fold05
## 6 <split [4.4K/486]> Fold06
## 7 <split [4.4K/486]> Fold07
## 8 <split [4.4K/486]> Fold08
## 9 <split [4.4K/485]> Fold09
## 10 <split [4.4K/484]> Fold10

다음으로 조정 가능한 결정 트리 모델 사양을 생성해 보겠습니다.

tree_spec <- decision_tree(
  cost_complexity = tune(),
  tree_depth = tune(),
  min_n = tune()
) %>%
  set_engine("rpart") %>%
  set_mode("regression")

tree_spec


## Decision Tree Model Specification (regression)
## 
## Main Arguments:
## cost_complexity = tune()
## tree_depth = tune()
## min_n = tune()
## 
## Computational engine: rpart

의사결정 트리에 대해 시도하려면 가능한 매개변수 값 세트가 필요합니다.

tree_grid <- grid_regular(cost_complexity(), tree_depth(), min_n(), levels = 4)

tree_grid


## # A tibble: 64 x 3
## cost_complexity tree_depth min_n
## <dbl> <int> <int>
## 1 0.0000000001 1 2
## 2 0.0000001 1 2
## 3 0.0001 1 2
## 4 0.1 1 2
## 5 0.0000000001 5 2
## 6 0.0000001 5 2
## 7 0.0001 5 2
## 8 0.1 5 2
## 9 0.0000000001 10 2
## 10 0.0000001 10 2
## # … with 54 more rows

이제 리샘플링된 모든 데이터세트에서 가능한 모든 매개변수 값을 시험해 보겠습니다. 기본이 아닌 몇 가지 측정항목을 사용하겠습니다.

doParallel::registerDoParallel()

set.seed(345)
tree_rs <- tune_grid(
  tree_spec,
  turbine_capacity ~ .,
  resamples = wind_folds,
  grid = tree_grid,
  metrics = metric_set(rmse, rsq, mae, mape)
)

tree_rs


## # Tuning results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 4
## splits id .metrics .notes          
## <list> <chr> <list> <list>          
## 1 <split [4.4K/488]> Fold01 <tibble [256 × 7]> <tibble [0 × 1]>
## 2 <split [4.4K/487]> Fold02 <tibble [256 × 7]> <tibble [0 × 1]>
## 3 <split [4.4K/486]> Fold03 <tibble [256 × 7]> <tibble [0 × 1]>
## 4 <split [4.4K/486]> Fold04 <tibble [256 × 7]> <tibble [0 × 1]>
## 5 <split [4.4K/486]> Fold05 <tibble [256 × 7]> <tibble [0 × 1]>
## 6 <split [4.4K/486]> Fold06 <tibble [256 × 7]> <tibble [0 × 1]>
## 7 <split [4.4K/486]> Fold07 <tibble [256 × 7]> <tibble [0 × 1]>
## 8 <split [4.4K/486]> Fold08 <tibble [256 × 7]> <tibble [0 × 1]>
## 9 <split [4.4K/485]> Fold09 <tibble [256 × 7]> <tibble [0 × 1]>
## 10 <split [4.4K/484]> Fold10 <tibble [256 × 7]> <tibble [0 × 1]>

제가 자주 보여드린 것처럼 여기서는 workflow()를 조정하지 않는다는 점에 유의하십시오. 대신 모델 사양을 조정하고 있습니다(수식 전처리기와 함께). 이는 일부 모델 평가 활동에서 베어 모델을 사용할 수 있도록 하기 위한 것입니다.

모델 평가

이제 우리가 어떻게 했는지 확인해 봅시다. 메트릭을 수집하거나 시각화할 수 있습니다.

collect_metrics(tree_rs)


## # A tibble: 256 x 9
## cost_complexity tree_depth min_n .metric .estimator mean n std_err
## <dbl> <int> <int> <chr> <chr> <dbl> <int> <dbl>
## 1 0.0000000001 1 2 mae standard 386. 10 1.50  
## 2 0.0000000001 1 2 mape standard 27.7 10 1.30  
## 3 0.0000000001 1 2 rmse standard 508. 10 1.44  
## 4 0.0000000001 1 2 rsq standard 0.303 10 0.0134
## 5 0.0000001 1 2 mae standard 386. 10 1.50  
## 6 0.0000001 1 2 mape standard 27.7 10 1.30  
## 7 0.0000001 1 2 rmse standard 508. 10 1.44  
## 8 0.0000001 1 2 rsq standard 0.303 10 0.0134
## 9 0.0001 1 2 mae standard 386. 10 1.50  
## 10 0.0001 1 2 mape standard 27.7 10 1.30  
## # … with 246 more rows, and 1 more variable: .config <chr>


autoplot(tree_rs) + theme_light(base_family = "IBMPlexSans")

이 데이터에는 상당히 복잡한 트리가 필요한 것 같습니다!

우리는 원하는 지표에 따라 선택된 최상의 매개 변수 옵션 세트를 검사하거나 선택할 수 있습니다.

show_best(tree_rs, "mape")


## # A tibble: 5 x 9
## cost_complexity tree_depth min_n .metric .estimator mean n std_err
## <dbl> <int> <int> <chr> <chr> <dbl> <int> <dbl>
## 1 0.0000000001 15 2 mape standard 0.564 10 0.0592
## 2 0.0000001 15 2 mape standard 0.564 10 0.0591
## 3 0.0000000001 15 14 mape standard 0.823 10 0.0547
## 4 0.0000001 15 14 mape standard 0.823 10 0.0547
## 5 0.0001 15 2 mape standard 0.885 10 0.0705
## # … with 1 more variable: .config <chr>


select_best(tree_rs, "rmse")


## # A tibble: 1 x 4
## cost_complexity tree_depth min_n .config              
## <dbl> <int> <int> <chr>                
## 1 0.0000000001 15 2 Preprocessor1_Model13

다음으로 이러한 "최고"의 매개변수 집합 중 하나를 사용하여 모델을 업데이트하고 마무리하겠습니다.

final_tree <- finalize_model(tree_spec, select_best(tree_rs, "rmse"))

final_tree


## Decision Tree Model Specification (regression)
## 
## Main Arguments:
## cost_complexity = 1e-10
## tree_depth = 15
## min_n = 2
## 
## Computational engine: rpart

이 모델final_tree은 업데이트되고 확정되었지만(더 이상 조정할 수 없음) 적합하지 않습니다. 모든 하이퍼파라미터가 설정되어 있지만 어떤 데이터에도 적합하지 않습니다. 이 모델에 맞추는 방법에 대한 몇 가지 옵션이 있습니다. final_tree를 사용하여 교육 데이터에 fit()를 맞추거나 last_fit()를 사용하여 테스트/훈련 분할에 맞출 수 있습니다. 그러면 맞춤 출력과 함께 다른 결과가 제공됩니다.

final_fit <- fit(final_tree, turbine_capacity ~ ., wind_train)
final_rs <- last_fit(final_tree, turbine_capacity ~ ., wind_split)

이러한 개체 중 하나에서 예측할 수 있습니다.

predict(final_fit, wind_train[144,])


## # A tibble: 1 x 1
## .pred
## <dbl>
## 1 1800


predict(final_rs$.workflow[[1]], wind_train[144,])


## # A tibble: 1 x 1
## .pred
## <dbl>
## 1 1800

터빈 용량 예측을 위한 이 결정 트리에서 가장 중요한 변수는 무엇입니까?

library(vip)

final_fit %>%
  vip(geom = "col", aesthetics = list(fill = "midnightblue", alpha = 0.8)) +
  scale_y_continuous(expand = c(0, 0))

시각화 의사 결정 트리 결과를 위한 parttree 패키지가 정말 마음에 듭니다. 하나 또는 두 개의 예측 변수가 있는 모델에서만 작동하므로 전체 모델과 완전히 동일하지 않은 예제 모델을 적합해야 합니다. 여전히 이 의사 결정 트리가 작동하는 방식을 이해하는 데 도움이 될 수 있지만 더 많은 예측 변수가 있는 전체 모델과 동일하지 않다는 점을 명심하십시오.

library(parttree)

ex_fit <- fit(
  final_tree,
  turbine_capacity ~ rotor_diameter_m + commissioning_date,
  wind_train
)

wind_train %>%
  ggplot(aes(rotor_diameter_m, commissioning_date)) +
  geom_parttree(data = ex_fit, aes(fill = turbine_capacity), alpha = 0.3) +
  geom_jitter(alpha = 0.7, width = 1, height = 0.5, aes(color = turbine_capacity)) +
  scale_colour_viridis_c(aesthetics = c("color", "fill"))

마지막으로 테스트 데이터를 살펴보겠습니다! 이러한 결과는 final_rs 에 적합 출력과 함께 저장됩니다. 테스트 데이터와 예측에 대한 메트릭을 모두 볼 수 있습니다.

collect_metrics(final_rs)


## # A tibble: 2 x 4
## .metric .estimator .estimate .config             
## <chr> <chr> <dbl> <chr>               
## 1 rmse standard 73.6 Preprocessor1_Model1
## 2 rsq standard 0.985 Preprocessor1_Model1


final_rs %>%
  collect_predictions() %>%
  ggplot(aes(turbine_capacity, .pred)) +
  geom_abline(slope = 1, lty = 2, color = "gray50", alpha = 0.5) +
  geom_point(alpha = 0.6, color = "midnightblue") +
  coord_fixed()