1 파이썬 scikit-learn ¹

기계학습 모형 개발을 위해서 전용 라이브러리 scikit-learn이 잘 준비되어 있다. 물론 R에서 기계학습을 위해서 The caret Package 팩키지가 있지만 개발된지 거의 20년이 되어가고 현재와 같은 상황을 반영하여 설계된 것이 아니라 RStudio로 자리를 옮긴 후에 본격적인 재설계와 새로운 바탕에서 새로 개발되고 있고 그 결과물 중 하나가 parsnip이다.

predictive models

2 `tidymodels` → `parsnip`

2.1 데이터

캐글 hr-comma-sep, “HR Analytics for employee rentension” 데이터셋을 가져와서 직원 이탈에 대한 예측모형을 개발한다.

library(tidyverse)
library(tidymodels)
library(skimr)
library(janitor)

hr_dat <- read_csv("data/HR_comma_sep.csv") %>% 
  clean_names()

hr_dat %>% 
  skim()

Skim summary statistics
 n obs: 14999 
 n variables: 10 

── Variable type:character ──────────────────────────────────────────────────────────────────────
    variable missing complete     n min max empty n_unique
 departments       0    14999 14999   2  11     0       10
      salary       0    14999 14999   3   6     0        3

── Variable type:numeric ────────────────────────────────────────────────────────────────────────
              variable missing complete     n    mean    sd    p0    p25
  average_montly_hours       0    14999 14999 201.05  49.94 96    156   
       last_evaluation       0    14999 14999   0.72   0.17  0.36   0.56
                  left       0    14999 14999   0.24   0.43  0      0   
        number_project       0    14999 14999   3.8    1.23  2      3   
 promotion_last_5years       0    14999 14999   0.021  0.14  0      0   
    satisfaction_level       0    14999 14999   0.61   0.25  0.09   0.44
    time_spend_company       0    14999 14999   3.5    1.46  2      3   
         work_accident       0    14999 14999   0.14   0.35  0      0   
    p50    p75 p100     hist
 200    245     310 ▁▇▇▆▆▇▆▂
   0.72   0.87    1 ▂▇▇▆▆▇▇▇
   0      0       1 ▇▁▁▁▁▁▁▂
   4      5       7 ▅▇▁▇▅▁▂▁
   0      0       1 ▇▁▁▁▁▁▁▁
   0.64   0.82    1 ▅▂▅▆▆▇▇▇
   3      4      10 ▇▂▁▁▁▁▁▁
   0      0       1 ▇▁▁▁▁▁▁▂

2.2 데이터 전처리

범주형 변수를 명시적으로 정의한다. 대표적으로 departments, left, work_accident를 범주형 변수로 빼고 나머지는 숫자형이라 예측모형을 돌리기 위한 최소한의 기준을 만족시켰다.

hr_df <- hr_dat %>% 
  mutate(left = factor(left, levels=c(0,1), labels=c("stay", "left"))) %>%
  mutate(departments = factor(departments),
         work_accident = factor(work_accident))

hr_df %>% 
  skim()

Skim summary statistics
 n obs: 14999 
 n variables: 10 

── Variable type:character ──────────────────────────────────────────────────────────────────────
 variable missing complete     n min max empty n_unique
   salary       0    14999 14999   3   6     0        3

── Variable type:factor ─────────────────────────────────────────────────────────────────────────
      variable missing complete     n n_unique
   departments       0    14999 14999       10
          left       0    14999 14999        2
 work_accident       0    14999 14999        2
                                top_counts ordered
 sal: 4140, tec: 2720, sup: 2229, IT: 1227   FALSE
              sta: 11428, lef: 3571, NA: 0   FALSE
                  0: 12830, 1: 2169, NA: 0   FALSE

── Variable type:numeric ────────────────────────────────────────────────────────────────────────
              variable missing complete     n    mean    sd    p0    p25
  average_montly_hours       0    14999 14999 201.05  49.94 96    156   
       last_evaluation       0    14999 14999   0.72   0.17  0.36   0.56
        number_project       0    14999 14999   3.8    1.23  2      3   
 promotion_last_5years       0    14999 14999   0.021  0.14  0      0   
    satisfaction_level       0    14999 14999   0.61   0.25  0.09   0.44
    time_spend_company       0    14999 14999   3.5    1.46  2      3   
    p50    p75 p100     hist
 200    245     310 ▁▇▇▆▆▇▆▂
   0.72   0.87    1 ▂▇▇▆▆▇▇▇
   4      5       7 ▅▇▁▇▅▁▂▁
   0      0       1 ▇▁▁▁▁▁▁▁
   0.64   0.82    1 ▅▂▅▆▆▇▇▇
   3      4      10 ▇▂▁▁▁▁▁▁

3 `tidymodels` 예측모형

3.1 훈련/시험 데이터 구분

가장 먼저 훈련/시험 데이터 구분한다. 이를 위해서 rsample 팩키지를 사용하고 initial_split() 함수를 사용해서 7:3으로 훈련/시험 데이터로 쪼갠다.

train_test_split <- rsample::initial_split(
  data = hr_df,     
  prop = 0.70)

train_test_split

<10500/4499/14999>

훈련/시험 데이터프레임으로 training(), testing() 함수를 각각 사용해서 준비한다.

train_df <- train_test_split %>% training() 
test_df  <- train_test_split %>% testing()

3.2 피처 공학(Feature engineering)

예측모형을 개발할 때 피처 공학(Feature engineering)을 통해 예측모형 base table을 제작해야 한다. tidymodels에서는 이를 요리법(recipe)에 비유하여 요리 레시피를 준비해야 한다. 대표적인 요리법에는 결측값 제거, 이상점 처리, 다공선성을 갖는 변수 제거, 중심화와 척도 조정을 통한 정규화, 가변수 처리 등이 포함된다.

범주형변수는 가변수화(step_dummy())하고 연속형 변수는 정규화하는 요리법을 hr_recipe로 준비한다. 그리고 이를 prep()으로 준비하여 bake()해서 예측모형에 넣을 수 있도록 구워둔다.

hr_recipe <- function(df) {
  recipe(left ~ ., data = df) %>%
    step_dummy(all_nominal(), -all_outcomes()) %>%
    step_center(all_numeric()) %>%
    step_scale(all_numeric()) %>%
    prep(data = df)
}

recipe_prepped <- hr_recipe(df = train_df)

train_baked <- bake(recipe_prepped, new_data = train_df)
test_baked  <- bake(recipe_prepped, new_data = test_df)

train_baked

# A tibble: 10,500 x 19
   satisfaction_le… last_evaluation number_project average_montly_…
              <dbl>           <dbl>          <dbl>            <dbl>
 1           -0.937          -1.09          -1.47            -0.892
 2            0.753           0.844          0.960            1.22 
 3            0.431           0.903          0.960            0.435
 4           -0.816          -1.26          -1.47            -0.973
 5           -2.06            0.317          1.77             0.917
 6           -0.655          -1.03          -1.47            -1.33 
 7            0.914           1.20           0.148            0.656
 8           -0.937          -1.03          -1.47            -1.17 
 9           -0.655          -1.44          -1.47            -0.832
10            0.673           1.61           0.148            1.08 
# … with 10,490 more rows, and 15 more variables:
#   time_spend_company <dbl>, left <fct>, promotion_last_5years <dbl>,
#   work_accident_X1 <dbl>, departments_hr <dbl>, departments_IT <dbl>,
#   departments_management <dbl>, departments_marketing <dbl>,
#   departments_product_mng <dbl>, departments_RandD <dbl>,
#   departments_sales <dbl>, departments_support <dbl>,
#   departments_technical <dbl>, salary_low <dbl>, salary_medium <dbl>

3.3 모형 적합

드디어 parsnip의 가장 인기있는 부분을 살펴볼 수 있는데 앞서 준비한 base table 데이터프레임을 다음 3단계로 준비하여 훈련데이터를 예측모형에 적합시킨다.

예측모형 유형: 로지스틱 회귀모형(logistic_reg), 모드(classification)
엔진(set_engine): glm
적합(fit): fit

예를 들어, logistic_reg에 사용가능한 엔진은 glm, glmnet, stan, spark, keras 등이 있다.

hr_glm <-
  logistic_reg(mode = "classification") %>%
  set_engine("glm") %>%
  fit(left ~ ., data = train_baked)

3.4 모형 성능평가

yardstick 팩키지를 사용하게 되면 예측모형의 성능을 수월히 평가할 수 있다.

3.4.1 시험데이터 예측

먼저 predict() 함수로 예측값을 뽑아내서 직원 잔존과 이탈여부를 데이터프레임으로 만들 수 있다.

hr_glm_pred_df <- hr_glm %>%
  predict(new_data = test_baked) %>%
  bind_cols(test_baked %>% select(left))

hr_glm_pred_df

# A tibble: 4,499 x 2
   .pred_class left 
   <fct>       <fct>
 1 left        left 
 2 left        left 
 3 stay        left 
 4 stay        left 
 5 stay        left 
 6 left        left 
 7 stay        left 
 8 left        left 
 9 stay        left 
10 left        left 
# … with 4,489 more rows

3.4.2 모형성능

conf_mat() 함수를 사용해서 이진 분류 성능예측모형을 확인할 수 있다. 물론, yardstick 팩키지 metrics() 함수를 사용해도 된다.

hr_glm_pred_df %>%
  conf_mat(left, .pred_class) %>% 
  pluck(1) %>%
  as_tibble() %>%
  spread(Truth, n)

# A tibble: 2 x 3
  Prediction  left  stay
  <chr>      <int> <int>
1 left         361   271
2 stay         668  3199

hr_glm_pred_df %>%
  metrics(left, .pred_class)

# A tibble: 2 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.791
2 kap      binary         0.316

실제 예측모형을 현업에 적용시킬 때, 정확도(accuracy) 보다는 precision, recall, 혹은 f1_score가 더 적합한 측도가 될 수 있다.

clf_metric_df <- tibble(
  "precision" = 
     hr_glm_pred_df %>% precision(left, .pred_class) %>%
     select(.estimate),
  "recall" = 
     hr_glm_pred_df %>% recall(left, .pred_class) %>%
     select(.estimate),
  "f1_score" = 
     hr_glm_pred_df %>% f_meas(left, .pred_class) %>%
     select(.estimate)
) %>%
  unnest()

clf_metric_df

# A tibble: 1 x 3
  precision recall f1_score
      <dbl>  <dbl>    <dbl>
1     0.827  0.922    0.872

4 `tidymodels` - `ranger` 예측모형

GLM 계열 예측모형은 2015년 전까지 가장 검증된 예측모형으로 자리를 잡았으나 그 후 random forest, xgboost 계열 앙상블 모형으로 대체되고 있지만, 기준 성능을 나타내는 지표를 제공한다는 점에서 나름 굳건히 자리를 지켜나가고 있다.

random forest 예측모형을 C/C++로 성능을 대폭 향상시킨 ranger 팩키지를 많이 사용하는데 사용업은 크게 두가지를 추가로 설정해야 된다. 하나는 교차검증(CV, cross-validation), set_engine을 설정하는 부분이다.

4.1 교차검증(CV)

vfold_cv() 함수를 사용해서 훈련데이터를 훈련/검증(train/validation) 데이터로 쪼갠다. CV를 10조각 내는 경우 훈련데이터를 9:1비율로 자동으로 쪼갠다.

훈련/시험 데이터 쪼개기

train_test_split

<10500/4499/14999>

훈련/교차검증 데이터 쪼개기

cross_val_df <- vfold_cv(train_df, v = 10)
cross_val_df

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

cross_val_df$splits %>% 
  pluck(7)

<9450/1050/10500>

4.2 피처 공학(feature engineering)

교차검증으로 준비된 훈련/검증 데이터셋에 대하여 피처 공학을 적용시키는데, 특별히 더 추가하지 않고 앞서 준비한 hr_recipe를 재사용한다.

hr_recipe

function(df) {
  recipe(left ~ ., data = df) %>%
    step_dummy(all_nominal(), -all_outcomes()) %>%
    step_center(all_numeric()) %>%
    step_scale(all_numeric()) %>%
    prep(data = df)
}

4.3 `ranger` 예측모형 적합

fit_ranger() 함수로 작성하고, 교차검정 데이터와 초모수 mtry와 tree를 넣어주는 함수를 작성시켜 적합시킨다.

fit_ranger <- function(split, id, mtry, tree) {
   
  analysis_set <- split %>% analysis()
  analysis_prepped <- analysis_set %>% hr_recipe()
  analysis_baked <- analysis_prepped %>% bake(new_data = analysis_set)
  
  model_rf <- rand_forest(
      mode = "classification",
      mtry = mtry,
      trees = tree) %>%
    set_engine("ranger",
      importance = "impurity") %>%
    fit(left ~ ., data = analysis_baked)
  
  assessment_set <- split %>% assessment()
  assessment_prepped <- assessment_set %>% hr_recipe()
  assessment_baked <- assessment_prepped %>% bake(new_data = assessment_set)
  
  tibble(
    "id" = id,
    "truth" = assessment_baked$left,
    "prediction" = model_rf %>%
      predict(new_data = assessment_baked) %>%
      unlist()
  )
}

pred_rf <- map2_df(
  .x = cross_val_df$splits,
  .y = cross_val_df$id,
  ~ fit_ranger(split = .x, id = .y, mtry = 3, tree = 200)
)

pred_rf

# A tibble: 10,500 x 3
   id     truth prediction
   <chr>  <fct> <fct>     
 1 Fold01 left  left      
 2 Fold01 left  left      
 3 Fold01 left  left      
 4 Fold01 left  left      
 5 Fold01 left  left      
 6 Fold01 left  left      
 7 Fold01 left  left      
 8 Fold01 left  left      
 9 Fold01 left  left      
10 Fold01 left  stay      
# … with 10,490 more rows

4.4 `ranger` 예측모형 성능

conf_mat() 함수와 summary() 함수를 파이프로 연결시키면 예측모형 성능에 대한 모든 지표를 추출할 수 있다. 가장 많이 사용하는 “accuracy”, “precision”, “recall”, “f_meas” 측도만 추출하여 향상된 성능을 파악한다.

pred_rf %>%
  conf_mat(truth, prediction) %>%
  summary() %>%
  select(-.estimator) %>%
  filter(.metric %in% c("accuracy", "precision", "recall", "f_meas"))

# A tibble: 4 x 2
  .metric   .estimate
  <chr>         <dbl>
1 accuracy      0.966
2 precision     0.961
3 recall        0.996
4 f_meas        0.978

5 `tidymodels` - `xgboost` 예측모형

5.1 `XGBoost` 모형적합

fit_rf() 함수에 XGBoost 팩키지 초모수 learn_rate, tree_depth, sample_size을 적용시키도록 예측모형을 boost_tree로 변경시키고, set_engine도 xgboost로 바꿔 훈련시킨다.

fit_xgb <- function(split, id, learn_rate, tree_depth, sample_size) {
   
  analysis_set <- split %>% analysis()
  analysis_prepped <- analysis_set %>% hr_recipe()
  analysis_baked <- analysis_prepped %>% bake(new_data = analysis_set)
  
  model_xgb <- boost_tree(
      mode = "classification",
      learn_rate = learn_rate,
      tree_depth = tree_depth,
      sample_size = sample_size) %>%
    set_engine("xgboost") %>%
    fit(left ~ ., data = analysis_baked)
  
  assessment_set <- split %>% assessment()
  assessment_prepped <- assessment_set %>% hr_recipe()
  assessment_baked <- assessment_prepped %>% bake(new_data = assessment_set)
  
  tibble(
    "id" = id,
    "truth" = assessment_baked$left,
    "prediction" = model_xgb %>%
      predict(new_data = assessment_baked) %>%
      unlist()
  )
}

pred_xgb <- map2_df(
  .x = cross_val_df$splits,
  .y = cross_val_df$id,
  ~ fit_xgb(split = .x, id = .y, learn_rate=0.3, tree_depth=5, sample_size=1)
)

pred_xgb

# A tibble: 10,500 x 3
   id     truth prediction
   <chr>  <fct> <fct>     
 1 Fold01 left  left      
 2 Fold01 left  left      
 3 Fold01 left  left      
 4 Fold01 left  left      
 5 Fold01 left  left      
 6 Fold01 left  left      
 7 Fold01 left  left      
 8 Fold01 left  left      
 9 Fold01 left  left      
10 Fold01 left  left      
# … with 10,490 more rows

XGBoost 예측모형 성능을 동일한 방식으로 확인한다.

pred_xgb %>%
  conf_mat(truth, prediction) %>%
  summary() %>%
  select(-.estimator) %>%
  filter(.metric %in% c("accuracy", "precision", "recall", "f_meas"))

# A tibble: 4 x 2
  .metric   .estimate
  <chr>         <dbl>
1 accuracy      0.963
2 precision     0.962
3 recall        0.990
4 f_meas        0.976

5.2 격자 탐색(grid search) ²

XGBoost 초모수 튜닝을 위한 격자탐색에 초모수 임의 탐색(random search)을 반영한 model_xgb XGBoost 모형과 xgb_grid 탐색 격자를 준비하여 이를 merge 시켜 spec_df를 준비시킨다.

더불어 vfold_cv() 함수를 통해 교차검증 데이터셋도 함께 준비시킨다.

## 초모수 -----
xgb_grid <- grid_random(
  learn_rate  %>% range_set(c(0.1,  0.01)), 
  tree_depth  %>% range_set(c( 3,  10)), 
  sample_size %>% range_set(c(0.5, 1)), 
  size = 5)

xgb_grid

# A tibble: 5 x 3
  learn_rate tree_depth sample_size
       <dbl>      <int>       <dbl>
1     0.0366          7         0.5
2     0.0608          6         0.5
3     0.0387          7         0.5
4     0.0602          4         0.5
5     0.0849          4         0.5

## XGBoost 모형 -----
model_xgb <- boost_tree(
    mode = "classification",
    learn_rate = varing(),
    tree_depth = varing(),
    sample_size = varing()) 

spec_df <- tibble(spec = merge(model_xgb, xgb_grid)) %>% 
  mutate(model_id = row_number())

spec_df

# A tibble: 5 x 2
  spec      model_id
  <list>       <int>
1 <spec[+]>        1
2 <spec[+]>        2
3 <spec[+]>        3
4 <spec[+]>        4
5 <spec[+]>        5

## CV 데이터셋 -----
cross_val_df <- vfold_cv(train_df, v = 3)
cross_val_df

#  3-fold cross-validation 
# A tibble: 3 x 2
  splits            id   
  <named list>      <chr>
1 <split [7K/3.5K]> Fold1
2 <split [7K/3.5K]> Fold2
3 <split [7K/3.5K]> Fold3

spec_cv_df <- crossing(spec_df, cross_val_df)
spec_cv_df

# A tibble: 15 x 4
   spec      model_id splits            id   
   <list>       <int> <named list>      <chr>
 1 <spec[+]>        1 <split [7K/3.5K]> Fold1
 2 <spec[+]>        1 <split [7K/3.5K]> Fold2
 3 <spec[+]>        1 <split [7K/3.5K]> Fold3
 4 <spec[+]>        2 <split [7K/3.5K]> Fold1
 5 <spec[+]>        2 <split [7K/3.5K]> Fold2
 6 <spec[+]>        2 <split [7K/3.5K]> Fold3
 7 <spec[+]>        3 <split [7K/3.5K]> Fold1
 8 <spec[+]>        3 <split [7K/3.5K]> Fold2
 9 <spec[+]>        3 <split [7K/3.5K]> Fold3
10 <spec[+]>        4 <split [7K/3.5K]> Fold1
11 <spec[+]>        4 <split [7K/3.5K]> Fold2
12 <spec[+]>        4 <split [7K/3.5K]> Fold3
13 <spec[+]>        5 <split [7K/3.5K]> Fold1
14 <spec[+]>        5 <split [7K/3.5K]> Fold2
15 <spec[+]>        5 <split [7K/3.5K]> Fold3

5.3 격자 탐색(grid search) 적합

앞서 준비한 초모수 XGBoost 모형과 CV 교차검증 데이터셋에 대한 적합을 시도함.

fit_on_fold <- function(df) {
  fit(left ~ ., data=df)
}

spec_cv_fit <- spec_cv_df %>%
  mutate(
    analysis_set = map(splits, rsample::analysis),
    prepped = map(splits, hr_recipe),
    baked   = map2(prepped, analysis_set, bake)
  ) %>% 
  mutate(
    xgb_fit = map(baked, ~fit_on_fold)
  )

xgb_predict <- function(model, newdata) {
  predict(model, newdata)
}

spec_cv_pred <- spec_cv_fit %>% 
  mutate(
    assess_set = map(splits, rsample::assessment),
    assess_prepped = map(splits, hr_recipe),
    assess_baked   = map2(prepped, assess_set, bake)
  ) %>% 
  mutate(
    pred = map2(xgb_fit, assess_baked, ~xgb_predict)
  )

xwMOOC 모형 - `tidymodels`

`caret` → `parsnip`

xwMOOC

2019-08-18

1 파이썬 scikit-learn ¹

2 `tidymodels` → `parsnip`

2.1 데이터

2.2 데이터 전처리

3 `tidymodels` 예측모형

3.1 훈련/시험 데이터 구분

3.2 피처 공학(Feature engineering)

3.3 모형 적합

3.4 모형 성능평가

3.4.1 시험데이터 예측

3.4.2 모형성능

4 `tidymodels` - `ranger` 예측모형

4.1 교차검증(CV)

4.2 피처 공학(feature engineering)

4.3 `ranger` 예측모형 적합

4.4 `ranger` 예측모형 성능

5 `tidymodels` - `xgboost` 예측모형

5.1 `XGBoost` 모형적합

5.2 격자 탐색(grid search) ²

5.3 격자 탐색(grid search) 적합

xwMOOC 모형 - tidymodels

caret → parsnip

xwMOOC

2019-08-18

1 파이썬 scikit-learn 1

2 tidymodels → parsnip

2.1 데이터

2.2 데이터 전처리

3 tidymodels 예측모형

3.1 훈련/시험 데이터 구분

3.2 피처 공학(Feature engineering)

3.3 모형 적합

3.4 모형 성능평가

3.4.1 시험데이터 예측

3.4.2 모형성능

4 tidymodels - ranger 예측모형

4.1 교차검증(CV)

4.2 피처 공학(feature engineering)

4.3 ranger 예측모형 적합

4.4 ranger 예측모형 성능

5 tidymodels - xgboost 예측모형

5.1 XGBoost 모형적합

5.2 격자 탐색(grid search) 2

5.3 격자 탐색(grid search) 적합

xwMOOC 모형 - `tidymodels`

`caret` → `parsnip`

1 파이썬 scikit-learn ¹

2 `tidymodels` → `parsnip`

3 `tidymodels` 예측모형

4 `tidymodels` - `ranger` 예측모형

4.3 `ranger` 예측모형 적합

4.4 `ranger` 예측모형 성능

5 `tidymodels` - `xgboost` 예측모형

5.1 `XGBoost` 모형적합

5.2 격자 탐색(grid search) ²