히말라야 등산 탐험 모형에서 계급 불균형을 처리하다⛰

최근에 나는 screencasts 프레임워크를 어떻게 사용하는지 시범을 보였다. 최초의 모델링 절차부터 더욱 복잡한 모델링을 조정했다.오늘의 스크린 프레젠테이션은 처음부터 끝까지 상세한 모델 분석으로 중요한 특징 공사 절차와 몇 가지 모델 유형을 포함하여 이번 주 tidymodels 히말라야 등산 탐험을 사용한다.
이것은 내가 동영상에서 사용하는 코드로 동영상이나 동영상 이외의 책을 즐겨 읽는 사람들에게 적용된다.

#TidyTuesday 데이터 세트 데이터 찾아보기

우리의 모델링 목표는 히말라야 탐험대원들의 생존 또는 사망 확률 을 예측하는 것이다.이 데이터 집합은 우리에게 특징적인 공사 절차를 토론할 수 있는 기회를 제공했다. 예를 들어 계급이 불균형한 하위 표본 추출(살아있는 사람이 죽은 사람보다 많다)과 부족한 데이터 삽입(예를 들어 많은 탐험대원들이 나이가 부족하다).
데이터를 읽는 것부터 시작합시다.

library(tidyverse)
members <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-09-22/members.csv")

members


## # A tibble: 76,519 x 21
## expedition_id member_id peak_id peak_name year season sex age
## <chr> <chr> <chr> <chr> <dbl> <chr> <chr> <dbl>
## 1 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 40
## 2 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 41
## 3 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 27
## 4 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 40
## 5 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 34
## 6 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 25
## 7 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 41
## 8 AMAD78301 AMAD7830… AMAD Ama Dabl… 1978 Autumn M 29
## 9 AMAD79101 AMAD7910… AMAD Ama Dabl… 1979 Spring M 35
## 10 AMAD79101 AMAD7910… AMAD Ama Dabl… 1979 Spring M 37
## # … with 76,509 more rows, and 13 more variables: citizenship <chr>,
## # expedition_role <chr>, hired <lgl>, highpoint_metres <dbl>, success <lgl>,
## # solo <lgl>, oxygen_used <lgl>, died <lgl>, death_cause <chr>,
## # death_height_metres <dbl>, injured <lgl>, injury_type <chr>,
## # injury_height_metres <dbl>

영상에서 나는 skimr::skim()의 결과를 훑어본 결과 어떤 변수가 데이터가 부족하고 공민 신분이나 산봉우리 등 변수가 얼마나 유일한 값이 있는지 알아차렸다.
시간의 추이에 따라 탐험대의 성공률과 대원의 사망률은 어떤 변화가 생겼습니까?

members %>%
  group_by(year = 10 * (year %/% 10)) %>%
  summarise(
    died = mean(died),
    success = mean(success)
  ) %>%
  pivot_longer(died:success, names_to = "outcome", values_to = "percent") %>%
  ggplot(aes(year, percent, color = outcome)) +
  geom_line(alpha = 0.7, size = 1.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = NULL, y = "% of expedition members", color = NULL)

based on characteristics of the person and climbing expedition from this week’s #TidyTuesday dataset
탐험대원의 나이는 탐험의 성공이나 죽음과 관계가 있습니까?우리는 같은 코드를 사용할 수 있지만, year 를 age 로 전환하기만 하면 된다.

members %>%
  group_by(age = 10 * (age %/% 10)) %>%
  summarise(
    died = mean(died),
    success = mean(success)
  ) %>%
  pivot_longer(died:success, names_to = "outcome", values_to = "percent") %>%
  ggplot(aes(age, percent, color = outcome)) +
  geom_line(alpha = 0.7, size = 1.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = NULL, y = "% of expedition members", color = NULL)

사람들이 성공하지 못한 탐험으로 죽을 가능성이 더 높습니까?

members %>%
  count(success, died) %>%
  group_by(success) %>%
  mutate(percent = scales::percent(n / sum(n))) %>%
  kable(
    col.names = c("Expedition success", "Died", "Number of people", "% of people"),
    align = "llrr"
  )

탐험 성공
사망
인원수
% 명 수
틀렸어
틀렸어
46452
98%
틀렸어
사실에 부합했어
868
2%
사실에 부합했어
틀렸어
28961
99%
사실에 부합했어
사실에 부합했어
238
1%
우리는 비슷한 방법으로 히말라야산의 산봉우리마다 사망률이 얼마나 다른지 관찰할 수 있다.

members %>%
  filter(!is.na(peak_name)) %>%
  mutate(peak_name = fct_lump(peak_name, prop = 0.05)) %>%
  count(peak_name, died) %>%
  group_by(peak_name) %>%
  mutate(percent = scales::percent(n / sum(n))) %>%
  kable(
    col.names = c("Peak", "Died", "Number of people", "% of people"),
    align = "llrr"
  )

봉우리
사망
인원수
% 명 수
아마다 브라운 봉
틀렸어
8374
100%
아마다 브라운 봉
사실에 부합했어
32
0%
Cho Oyu
틀렸어
8838
99%
Cho Oyu
사실에 부합했어
52
1%
에베레스트 산
틀렸어
21507
99%
에베레스트 산
사실에 부합했어
306
1%
마나스루
틀렸어
4508
98%
마나스루
사실에 부합했어
85
2%
따로
틀렸어
32171
98%
따로
사실에 부합했어
631
2%
마지막 탐색적인 줄거리를 만들어 계절을 봅시다.사계절 간의 생존율은 얼마나 차이가 있습니까?

members %>%
  filter(season != "Unknown") %>%
  count(season, died) %>%
  group_by(season) %>%
  mutate(
    percent = n / sum(n),
    died = case_when(
      died ~ "Died",
      TRUE ~ "Did not die"
    )
  ) %>%
  ggplot(aes(season, percent, fill = season)) +
  geom_col(alpha = 0.8, position = "dodge", show.legend = FALSE) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_wrap(~died, scales = "free") +
  labs(x = NULL, y = "% of expedition members")

그리고 #Tidy Tuesday EDA에 대한 좋은 예는 탐색할 수 있습니다!현재, 우리는 일부 변수를 필터하고 일부 변수를 인자로 바꾸어 모델링하는 데이터 집합을 만듭니다.나이는 여전히 NA치가 많지만, 우리는 그것을 보충할 것이다.

members_df <- members %>%
  filter(season != "Unknown", !is.na(sex), !is.na(citizenship)) %>%
  select(peak_id, year, season, sex, age, citizenship, hired, success, died) %>%
  mutate(died = case_when(
    died ~ "died",
    TRUE ~ "survived"
  )) %>%
  mutate_if(is.character, factor) %>%
  mutate_if(is.logical, as.integer)

members_df


## # A tibble: 76,507 x 9
## peak_id year season sex age citizenship hired success died    
## <fct> <dbl> <fct> <fct> <dbl> <fct> <int> <int> <fct>   
## 1 AMAD 1978 Autumn M 40 France 0 0 survived
## 2 AMAD 1978 Autumn M 41 France 0 0 survived
## 3 AMAD 1978 Autumn M 27 France 0 0 survived
## 4 AMAD 1978 Autumn M 40 France 0 0 survived
## 5 AMAD 1978 Autumn M 34 France 0 0 survived
## 6 AMAD 1978 Autumn M 25 France 0 0 survived
## 7 AMAD 1978 Autumn M 41 France 0 0 survived
## 8 AMAD 1978 Autumn M 29 France 0 0 survived
## 9 AMAD 1979 Spring M 35 USA 0 0 survived
## 10 AMAD 1979 Spring M 37 W Germany 0 1 survived
## # … with 76,497 more rows

모델 구축

우리는tidymodels 패키지를 불러오는 것부터 시작하여 데이터를 훈련 집합과 테스트 집합으로 나눌 수 있습니다.

library(tidymodels)

set.seed(123)
members_split <- initial_split(members_df, strata = died)
members_train <- training(members_split)
members_test <- testing(members_split)

우리는 를 사용하여 모델의 성능을 평가할 것이기 때문에 이러한 샘플링 집합을 준비할 것이다.

set.seed(123)
members_folds <- vfold_cv(members_train, strata = died)
members_folds


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

다음으로, 우리는 데이터를 미리 처리하는 설계도를 구축한다.

우선, 우리는 recipe() 우리의 모델이 무엇인지(여기 공식을 사용함)와 우리의 훈련 데이터가 무엇인지 알려야 한다.

다음으로 우리는 훈련 데이터가 집중된 연령 중위수로 age의 결핍치를 추산한다.더 복잡한 resampling도 있지만, 우리는 간단한 옵션을 계속 사용할 것이다.

다음으로 우리는 step_other() 최고치와 공민의 등급을 축소하는 것을 사용한다.이 단계 이전에 변수마다 수백 개의 값이 있습니다.

그 다음에 우리는 비디지털 분류 값에 대한 지표 변수를 만들 수 있지만, 결과died를 제외하고는 이를 하나의 인자로 보존해야 한다.

마지막으로 탐험에서 살아남은 사람은 죽은 사람steps available for imputation보다 훨씬 많다.

대상step_smote()은 데이터에 따라 훈련되지 않은 설계도(예를 들어 접어야 할 분류 단계가 계산되지 않은 것)이다.

library(themis)

members_rec <- recipe(died ~ ., data = members_train) %>%
  step_medianimpute(age) %>%
  step_other(peak_id, citizenship) %>%
  step_dummy(all_nominal(), -died) %>%
  step_smote(died)

members_rec


## Data Recipe
## 
## Inputs:
## 
## role #variables
## outcome 1
## predictor 8
## 
## Operations:
## 
## Median Imputation for age
## Collapsing factor levels for peak_id, citizenship
## Dummy variables from all_nominal(), -died
## SMOTE based on died

우리는 amembers_rec에서 이 레시피를 사용할 것이기 때문에 workflow() 여부를 너무 강조할 필요가 없다.도안이 데이터에 대한 작용을 탐색하려면 먼저 prep() 도안을 작성하여 각 단계에 필요한 매개 변수를 추정한 다음에 prep() 이 절차를 응용한 훈련 데이터를 추출할 수 있습니다.
논리 회귀 모델과 무작위 삼림 모델 두 가지를 비교합시다.이것은 내가 사용하는 두 가지 모델we will use bake(new_data = NULL) to balance the classes이다.우리는 모델 규범을 만드는 것부터 시작한다.

glm_spec <- logistic_reg() %>%
  set_engine("glm")

glm_spec


## Logistic Regression Model Specification (classification)
## 
## Computational engine: glm


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

rf_spec


## Random Forest Model Specification (classification)
## 
## Main Arguments:
## trees = 1000
## 
## Computational engine: ranger

다음은 티디모델스workflow()를 조합해 봅시다. 이것은 보조 대상으로 모델링 파이프를 관리하고 각 부분을 레고 블록처럼 조합하는 데 사용됩니다.참고, 아직 모델이 없습니다.Model: None.

members_wf <- workflow() %>%
  add_recipe(members_rec)

members_wf


## ══ Workflow ═══════════════════════════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: None
## 
## ── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## ● step_medianimpute()
## ● step_other()
## ● step_dummy()
## ● step_smote()

이제 우리는 모델을 추가하고 모든 샘플링을 의합할 수 있다.우선, 우리는 논리 회귀 모델을 의합할 수 있다.민감도와 특이성을 추가하기 위해 비기본 도량 집합을 설정합니다.

members_metrics <- metric_set(roc_auc, accuracy, sensitivity, specificity)

doParallel::registerDoParallel()
glm_rs <- members_wf %>%
  add_model(glm_spec) %>%
  fit_resamples(
    resamples = members_folds,
    metrics = members_metrics,
    control = control_resamples(save_pred = TRUE)
  )

glm_rs


## # Resampling results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 5
## splits id .metrics .notes .predictions       
## <list> <chr> <list> <list> <list>             
## 1 <split [51.6K/5.7K… Fold01 <tibble [4 × 3… <tibble [0 × … <tibble [5,739 × 5…
## 2 <split [51.6K/5.7K… Fold02 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 3 <split [51.6K/5.7K… Fold03 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 4 <split [51.6K/5.7K… Fold04 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 5 <split [51.6K/5.7K… Fold05 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 6 <split [51.6K/5.7K… Fold06 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 7 <split [51.6K/5.7K… Fold07 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 8 <split [51.6K/5.7K… Fold08 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 9 <split [51.6K/5.7K… Fold09 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 10 <split [51.6K/5.7K… Fold10 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…

그 다음으로 우리는 무작위 삼림 모형을 작성할 수 있다.

rf_rs <- members_wf %>%
  add_model(rf_spec) %>%
  fit_resamples(
    resamples = members_folds,
    metrics = members_metrics,
    control = control_resamples(save_pred = TRUE)
  )

rf_rs


## # Resampling results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 x 5
## splits id .metrics .notes .predictions       
## <list> <chr> <list> <list> <list>             
## 1 <split [51.6K/5.7K… Fold01 <tibble [4 × 3… <tibble [0 × … <tibble [5,739 × 5…
## 2 <split [51.6K/5.7K… Fold02 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 3 <split [51.6K/5.7K… Fold03 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 4 <split [51.6K/5.7K… Fold04 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 5 <split [51.6K/5.7K… Fold05 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 6 <split [51.6K/5.7K… Fold06 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 7 <split [51.6K/5.7K… Fold07 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 8 <split [51.6K/5.7K… Fold08 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 9 <split [51.6K/5.7K… Fold09 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…
## 10 <split [51.6K/5.7K… Fold10 <tibble [4 × 3… <tibble [0 × … <tibble [5,738 × 5…

우리는 이미 모든 후보 모델을 다시 샘플링한 훈련집과 일치시켰다!

파머 펭귄 팀의 댓글에. 평가 모델

이제 우리가 어떻게 하는지 봅시다.

collect_metrics(glm_rs)


## # A tibble: 4 x 5
## .metric .estimator mean n std_err
## <chr> <chr> <dbl> <int> <dbl>
## 1 accuracy binary 0.619 10 0.00230
## 2 roc_auc binary 0.705 10 0.00721
## 3 sens binary 0.678 10 0.0139 
## 4 spec binary 0.619 10 0.00243

좋아, 이것은 중등이지만, 적어도 정류와 음류에 대해서는 기본적으로 일치한다.함수collect_metrics()는 샘플링 결과에서 .metrics열을 추출하여 포맷합니다. 우리가 여기서 본 것처럼.

collect_metrics(rf_rs)


## # A tibble: 4 x 5
## .metric .estimator mean n std_err
## <chr> <chr> <dbl> <int> <dbl>
## 1 accuracy binary 0.972 10 0.000514
## 2 roc_auc binary 0.746 10 0.00936 
## 3 sens binary 0.164 10 0.0125  
## 4 spec binary 0.984 10 0.000499

정확도가 높은데 그 감도가... 아이고!설령 샘플링 전략을 채택했다 하더라도 무작위 삼림 모델은 이 두 종류를 어떻게 식별하는지 잘 배우지 못했다.이 모델들이 어떻게 이 점을 하는지 더욱 깊이 있게 이해합시다.예를 들어 그들은 어떻게 이 두 종류를 예측합니까?

glm_rs %>%
  conf_mat_resampled()


## # A tibble: 4 x 3
## Prediction Truth Freq
## <fct> <fct> <dbl>
## 1 died died 55.5
## 2 died survived 2157. 
## 3 survived died 26.5
## 4 survived survived 3499.


rf_rs %>%
  conf_mat_resampled()


## # A tibble: 4 x 3
## Prediction Truth Freq
## <fct> <fct> <dbl>
## 1 died died 13.5
## 2 died survived 89.9
## 3 survived died 68.5
## 4 survived survived 5566.

무작위 삼림 모델은 어떤 탐험대원의 사망을 확정하는 데 상당히 나쁘지만, 논리 회귀 모델은 이 두 부류에 대한 작용이 대체적으로 같다.
ROC 커브를 그릴 수도 있습니다.

glm_rs %>%
  collect_predictions() %>%
  group_by(id) %>%
  roc_curve(died, .pred_died) %>%
  ggplot(aes(1 - specificity, sensitivity, color = id)) +
  geom_abline(lty = 2, color = "gray80", size = 1.5) +
  geom_path(show.legend = FALSE, alpha = 0.6, size = 1.2) +
  coord_equal()

우리는 마침내 시험집으로 돌아갈 때가 되었다.전체 분석 과정에서 우리는 아직 테스트 집합을 사용하지 않았음을 주의한다.우리는 모델을 비교하고 평가하기 위해 훈련집의 샘플링을 사용했다.훈련 데이터를 다시 한 번 작성하고 함수last_fit()를 사용하여 테스트 데이터를 평가합시다.

members_final <- members_wf %>%
  add_model(glm_spec) %>%
  last_fit(members_split)

members_final


## # Resampling results
## # Manual resampling 
## # A tibble: 1 x 6
## splits id .metrics .notes .predictions .workflow
## <list> <chr> <list> <list> <list> <list>   
## 1 <split [57.4K… train/test… <tibble [2 … <tibble [0… <tibble [19,126… <workflo…

이곳의 지표와 예측은 테스트 데이터를 바탕으로 한다.

collect_metrics(members_final)


## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.619
## 2 roc_auc binary 0.689


collect_predictions(members_final) %>%
  conf_mat(died, .pred_class)


## Truth
## Prediction died survived
## died 196 7204
## survived 90 11636

훈련 데이터로 계수를 추정했다.만약 우리가 사용한다면tidy(), 우리는 우세비례가 있다.

members_final %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy(exponentiate = TRUE) %>%
  arrange(estimate) %>%
  kable(digits = 3)

학기
어림잡다
표준 오차
통계 자료
p, 가치
(차단)
0
0.944
-57.309
0
peak\u id\u MANA
0.220
0.042
-35.769
0
peak\u id\u 기타
0.223
0.034
-43.635
0
peak\u id\u EVER
0.294
0.036
-33.641
0
고용
0.497
0.054
-12.928
0
성별
0.536
0.029
-21.230
0
공민 신분과 기타
0.612
0.032
-15.299
0
일본 시민
0.739
0.038
-7.995
0
겨울
0.747
0.041
-7.180
0
네팔 시민
0.776
0.062
-4.128
0
peak\u id\u CHOY
0.865
0.043
-3.404
0.001
봄철
0.905
0.016
-6.335
0
나이.
0.991
0.001
-12.129
0
연중
1.029
0
59.745
0
미국 시민
1.334
0.043
6.759
0
영국 시민
1.419
0.045
7.858
0
성공
2.099
0.016
46.404
0
여름
7.142
0.092
21.433
0
우리는 또 이런 결과를 가시화할 수 있다.

members_final %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy() %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(estimate, fct_reorder(term, estimate))) +
  geom_vline(xintercept = 0, color = "gray50", lty = 2, size = 1.2) +
  geom_errorbar(aes(
    xmin = estimate - std.error,
    xmax = estimate + std.error
  ),
  width = .2, color = "gray50", alpha = 0.7
  ) +
  geom_point(size = 2, color = "#85144B") +
  labs(y = NULL, x = "Coefficent from logistic regression")

계수는 정상적인 특징(예를 들어 여름 등반, 성공적인 탐험 또는 영국이나 미국에서 온 것)과 생존과 관련이 있다.

계수가 마이너스인 특징(예를 들어 에베레스트 산을 포함한 특정 산봉우리를 등반하거나 탐험대의 고용 구성원 중 하나 또는 남성으로서)은 사망과 관련이 있다.

기억해라. 우리는 반드시 우리 모델의 예측 정밀도에 따라 이런 모델 계수를 해석해야 한다. 이것은 약간 중등이다.우리가 모델에서 직접 해석한 요소에 비해 이런 탐험에서 살아남은 요소가 더 많다.또 주의해야 할 것은 우리가 이 모델에서 본 증거에 의하면 탐험대원으로 네팔에 고용된 현지 샤를바 등산객이 얼마나 위험한지, as.