ヒマラヤを登る遠征モデルの不均衡を扱う⛰

最近、私は出版されましたscreencasts 使い方のデモtidymodels フレームワークは、より複雑なモデルを調整するために最初のモデリングステップを開始から.今日のスクリーンキャストは最初から最後まで詳細なモデル分析を通して、重要な機能のエンジニアリングステップといくつかのモデルタイプで、今週の #TidyTuesday dataset ヒマラヤ登り探検
ここでは、ビデオの代わりに、またはビデオに加えて読書を好む人々のために、私はビデオで使用されるコードです.

データを調査する

我々のモデリング目標は、ヒマラヤの遠征メンバーの生存確率を予測することであるbased on characteristics of the person and climbing expedition from this week’s #TidyTuesday dataset . このデータセットは、クラスの不均衡(多くの人々が死ぬより生き残る)のためのサブサンプリングのような機能工学のステップについて話す機会を与えて、不足しているデータ(例えば多くの探検隊員が年齢がなくなっている)を与えます.
データを読み始めましょう.

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)

遠征メンバーの年齢と遠征か死の関係はありますか?我々は同じコードを使用することができますがyear for 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 %
趙鴎
偽
8838年
99 %
趙鴎
真
52
1 %
エベレスト
偽
21507年
99 %
エベレスト
真
306
1 %
マナスル
偽
4508年
98 %
マナスル
真
85
2 %
その他
偽
32171年
98 %
その他
真
631年
2 %
つの最後の探検的な計画をして、季節を見ましょう.どのような違いは、4シーズン間での生存には何ですか?

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

そこに多くのより多くの偉大な例があります.では、変数の一部をフィルタリングして変数を変数に変換することでモデリングに使用するデータセットを作成しましょう.まだたくさんある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

モデルを作る

私たちは、TidyModel Metapackageをロードし、トレーニングとテストセットにデータを分割することから始めることができます.

library(tidymodels)

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

私たちは使用するつもりですresampling モデル性能を評価するために、それらの再標本化されたセットを準備しましょう.

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 トレーニングデータセットにおける中央値年齢の使用もっと複雑steps available for imputation , しかし、我々はここで簡単なオプションに固執します.

次に、私たちはstep_other() ピークと市民権のための分類レベルを崩壊させる.このステップの前に、変数の何百もの値がありました.

この後、結果を除いて、非数値、カテゴリの値のインジケータ変数を作成できますdied 我々は、要因として維持する必要があります.

最後に、死んだ人よりも多くの人々が彼らの遠征を生き延びたwe will use step_smote() to balance the classes .

オブジェクトmembers_rec データにまだ訓練されていないレシピ(例えば、どのカテゴリー的なレベルが崩壊しなければならないかは計算されませんでした).

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

我々は、このレシピを使用するつもりですworkflow() だから私たちはどうするかについて多くのストレスをする必要はありませんprep() かどうか.あなたがレシピをあなたのデータにしているものを調査したいならば、あなたは最初にそうすることができますprep() 各ステップごとに必要なパラメータを推定するためのレシピbake(new_data = NULL) これらのステップを適用してトレーニングデータを引き出す.
つの異なるモデル、ロジスティック回帰モデルとランダム森モデルを比較しましょうこれらは私が使用したのと同じモデルですin the post on the Palmer penguins . モデル仕様を作成することから始めます.

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

次は一緒にtidymodelを入れ始める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()

今、我々はモデルを加えて、各々のResamplesに合うものを加えることができます.まず,ロジスティック回帰モデルに適合できる.感度と特異性を加えることができるように、デフォルトのメトリックセットを設定しましょう.

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…

我々は、我々のResampledトレーニングセットの各々の候補モデルに合いました!

評価モデル

さあ、どうやってチェックしましょう.

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 我々がここにあるもののようなResampling結果からのコラム.

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

精度は素晴らしいですが、その感度…!ランダムフォレストモデルは、どのように両方のクラスを認識する学習の素晴らしい仕事をしていない、我々のオーバーサンプリング戦略でも.より深くこれらのモデルがこれを見るためにしている方法に掘りましょう.たとえば、2つのクラスをどのように予測していますか?

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

最終的にテストセットに戻る時間です.この分析の間、まだテストセットを使用していないことに注意してくださいモデルを比較し評価するために,訓練セットの再サンプリングを使用した.トレーニングデータにもう1回合いましょう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() ) トレーニングデータを使用して推定されている.我々が使うならばexponentiate = TRUE , オッズ比があります.

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

用語
見積り
stdエラー
統計量
p値
(迎撃)
0.000
0.944
- 57.309
0.000
ペクキンスイディントマナ
0.220
0.042
- 35.769
0.000
他
0.223
0.034
- 43.635
0.000
これまでに
0.294
0.036
- 33.641
0.000
雇われる
0.497
0.054
- 12.928
0.000
雌雄
0.536
0.029
- 21.230
0.000
その他
0.612
0.032
- 15.299
0.000
日本
0.739
0.038
- 7.995
0.000
季節冬
0.747
0.041
- 7.180
0.000
ネパール
0.776
0.062
- 4.128
0.000
ペコッキ・アイランド・チョイ
0.865
0.043
- 3.404
0.001
季節風泉
0.905
0.016
- 6.335
0.000
年齢
0.991
0.001
- 12.129
0.000
年
1.029
0.000
59.745
0.000
米国
1.334
0.043
6.759
0.000
イギリス
1.419
0.045
7.858
0.000
成功
2.099
0.016
46.404
0.000
夏季
7.142
0.092
21.433
0.000
これらの結果を可視化することもできる.

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

ポジティブな面での係数(夏に登るように、成功した探検であるか、英国か我々からであるように)の特徴は、生き残ることと関連します.

負の側の係数(エベレストを含む登山特定のピークのように、遠征の雇われたメンバーの一人であるか、男性である)は、死に関連しています.

このようにモデル係数を解釈する必要があることを覚えておいてください.誰がモデルで直接説明したものよりも、これらの遠征を生き残る人の遊びのより多くの要因があります.また、私たちは、このモデルで証拠を見ることができます.それは、ネパールのネイティブ・シェルパクライマーであり、遠征メンバーとして雇われていることですpointed out in this week’s #TidyTuesday README .