In this vignette we demonstrate how the various packages in the
TidyConsultant universe work together in a data science workflow by
analyzing insurance data set, which consists of a population of people
described by their characteristics and their insurance charge.
Analyze data properties with validata
insurance %>%
diagnose()
#> # A tibble: 7 × 6
#> variables types missing_count missing_percent unique_count unique_rate
#> <chr> <chr> <int> <dbl> <int> <dbl>
#> 1 age integer 0 0 47 0.0351
#> 2 sex character 0 0 2 0.00149
#> 3 bmi numeric 0 0 548 0.410
#> 4 children integer 0 0 6 0.00448
#> 5 smoker character 0 0 2 0.00149
#> 6 region character 0 0 4 0.00299
#> 7 charges numeric 0 0 1337 0.999
insurance %>%
diagnose_numeric()
#> # A tibble: 4 × 11
#> variables zeros minus infs min mean max `|x|<1 (ratio)` integer_ratio
#> <chr> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 age 0 0 0 18 3.92e1 6.4 e1 0 1
#> 2 bmi 0 0 0 16.0 3.07e1 5.31e1 0 0.0224
#> 3 children 574 0 0 0 1.09e0 5 e0 0.429 1
#> 4 charges 0 0 0 1122. 1.33e4 6.38e4 0 0.000747
#> # ℹ 2 more variables: mode <dbl>, mode_ratio <dbl>
insurance %>%
diagnose_category(everything(), max_distinct = 7) %>%
print(width = Inf)
#> # A tibble: 14 × 4
#> column level n ratio
#> <chr> <chr> <int> <dbl>
#> 1 sex male 676 0.505
#> 2 sex female 662 0.495
#> 3 children 0 574 0.429
#> 4 children 1 324 0.242
#> 5 children 2 240 0.179
#> 6 children 3 157 0.117
#> 7 children 4 25 0.0187
#> 8 children 5 18 0.0135
#> 9 smoker no 1064 0.795
#> 10 smoker yes 274 0.205
#> 11 region southeast 364 0.272
#> 12 region northwest 325 0.243
#> 13 region southwest 325 0.243
#> 14 region northeast 324 0.242
We can infer that each row represents a person, uniquely identified
by 5 characteristics. The charges column is also unique, but simply
because it is a high-precision double.
insurance %>%
determine_distinct(everything())
#> database has 1 duplicate rows, and will eliminate them
Analyze data relationships with autostats
insurance %>%
auto_cor(sparse = TRUE) -> cors
cors
#> # A tibble: 20 × 6
#> x y cor p.value significance method
#> <chr> <chr> <dbl> <dbl> <chr> <chr>
#> 1 smoker_no charges -0.787 8.27e-283 *** pearson
#> 2 smoker_yes charges 0.787 8.27e-283 *** pearson
#> 3 charges age 0.299 4.89e- 29 *** pearson
#> 4 region_southeast bmi 0.270 8.71e- 24 *** pearson
#> 5 charges bmi 0.198 2.46e- 13 *** pearson
#> 6 region_northeast bmi -0.138 3.91e- 7 *** pearson
#> 7 region_northwest bmi -0.136 5.94e- 7 *** pearson
#> 8 bmi age 0.109 6.19e- 5 *** pearson
#> 9 smoker_no sex_female 0.0762 5.30e- 3 ** pearson
#> 10 smoker_yes sex_female -0.0762 5.30e- 3 ** pearson
#> 11 smoker_yes sex_male 0.0762 5.30e- 3 ** pearson
#> 12 smoker_no sex_male -0.0762 5.30e- 3 ** pearson
#> 13 region_southeast charges 0.0740 6.78e- 3 ** pearson
#> 14 region_southeast smoker_no -0.0685 1.22e- 2 * pearson
#> 15 region_southeast smoker_yes 0.0685 1.22e- 2 * pearson
#> 16 charges children 0.0680 1.29e- 2 * pearson
#> 17 sex_male charges 0.0573 3.61e- 2 * pearson
#> 18 sex_female charges -0.0573 3.61e- 2 * pearson
#> 19 sex_male bmi 0.0464 9.00e- 2 . pearson
#> 20 sex_female bmi -0.0464 9.00e- 2 . pearson
Analyze the relationship between categorical and continuous
variables
insurance %>%
auto_anova(everything(), baseline = "first_level") -> anovas1
anovas1 %>%
print(n = 50)
#> # A tibble: 32 × 12
#> target predictor level estimate target_mean n std.error level_p.value
#> <chr> <chr> <chr> <dbl> <dbl> <int> <dbl> <dbl>
#> 1 age region (Inter… 3.93e+1 39.3 324 0.781 4.75e-310
#> 2 age region northw… -7.16e-2 39.2 325 1.10 9.48e- 1
#> 3 age region southe… -3.29e-1 38.9 364 1.07 7.59e- 1
#> 4 age region southw… 1.87e-1 39.5 325 1.10 8.66e- 1
#> 5 age sex (Inter… 3.95e+1 39.5 662 0.546 0
#> 6 age sex male -5.86e-1 38.9 676 0.768 4.46e- 1
#> 7 age smoker (Inter… 3.94e+1 39.4 1064 0.431 0
#> 8 age smoker yes -8.71e-1 38.5 274 0.952 3.60e- 1
#> 9 bmi region (Inter… 2.92e+1 29.2 324 0.325 0
#> 10 bmi region northw… 2.63e-2 29.2 325 0.459 9.54e- 1
#> 11 bmi region southe… 4.18e+0 33.4 364 0.447 3.28e- 20
#> 12 bmi region southw… 1.42e+0 30.6 325 0.459 1.99e- 3
#> 13 bmi sex (Inter… 3.04e+1 30.4 662 0.237 0
#> 14 bmi sex male 5.65e-1 30.9 676 0.333 9.00e- 2
#> 15 bmi smoker (Inter… 3.07e+1 30.7 1064 0.187 0
#> 16 bmi smoker yes 5.67e-2 30.7 274 0.413 8.91e- 1
#> 17 charges region (Inter… 1.34e+4 13406. 324 671. 7.67e- 78
#> 18 charges region northw… -9.89e+2 12418. 325 949. 2.97e- 1
#> 19 charges region southe… 1.33e+3 14735. 364 923. 1.50e- 1
#> 20 charges region southw… -1.06e+3 12347. 325 949. 2.64e- 1
#> 21 charges sex (Inter… 1.26e+4 12570. 662 470. 1.63e-126
#> 22 charges sex male 1.39e+3 13957. 676 661. 3.61e- 2
#> 23 charges smoker (Inter… 8.43e+3 8434. 1064 229. 1.58e-205
#> 24 charges smoker yes 2.36e+4 32050. 274 506. 8.27e-283
#> 25 children region (Inter… 1.05e+0 1.05 324 0.0670 1.26e- 50
#> 26 children region northw… 1.01e-1 1.15 325 0.0947 2.84e- 1
#> 27 children region southe… 3.15e-3 1.05 364 0.0921 9.73e- 1
#> 28 children region southw… 9.52e-2 1.14 325 0.0947 3.15e- 1
#> 29 children sex (Inter… 1.07e+0 1.07 662 0.0469 2.66e- 98
#> 30 children sex male 4.14e-2 1.12 676 0.0659 5.30e- 1
#> 31 children smoker (Inter… 1.09e+0 1.09 1064 0.0370 1.25e-147
#> 32 children smoker yes 2.29e-2 1.11 274 0.0817 7.79e- 1
#> # ℹ 4 more variables: level_significance <chr>, predictor_p.value <dbl>,
#> # predictor_significance <chr>, conclusion <chr>
Use the sparse option to show only the most significant effects
insurance %>%
auto_anova(everything(), baseline = "first_level", sparse = T, pval_thresh = .1) -> anovas2
anovas2 %>%
print(n = 50)
#> # A tibble: 5 × 8
#> target predictor level estimate target_mean n level_p.value
#> <chr> <chr> <chr> <dbl> <dbl> <int> <dbl>
#> 1 bmi region southeast 4.18 33.4 364 3.28e- 20
#> 2 bmi region southwest 1.42 30.6 325 1.99e- 3
#> 3 bmi sex male 0.565 30.9 676 9.00e- 2
#> 4 charges sex male 1387. 13957. 676 3.61e- 2
#> 5 charges smoker yes 23616. 32050. 274 8.27e-283
#> # ℹ 1 more variable: level_significance <chr>
From this we can conclude that smokers and males incur higher
insurance charges. It may be beneficial to explore some interaction
effects, considering that BMI’s effect on charges will be
sex-dependent.
create dummies out of categorical variables
insurance %>%
create_dummies(remove_most_frequent_dummy = T) -> insurance1
explore model-based variable contributions
insurance1 %>%
tidy_formula(target = charges) -> charges_form
charges_form
#> charges ~ age + bmi + children + sex_female + smoker_yes + region_northeast +
#> region_northwest + region_southwest
#> <environment: 0x11dfd80e0>
insurance1 %>%
auto_variable_contributions(formula = charges_form)
insurance1 %>%
auto_model_accuracy(formula = charges_form, include_linear = T)
create bins with tidybins
insurance1 %>%
bin_cols(charges) -> insurance_bins
insurance_bins
#> # A tibble: 1,338 × 10
#> charges_fr10 age bmi children charges sex_female smoker_yes
#> <int> <int> <dbl> <int> <dbl> <int> <int>
#> 1 8 19 27.9 0 16885. 1 1
#> 2 1 18 33.8 1 1726. 0 0
#> 3 3 28 33 3 4449. 0 0
#> 4 9 33 22.7 0 21984. 0 0
#> 5 2 32 28.9 0 3867. 0 0
#> 6 2 31 25.7 0 3757. 1 0
#> 7 5 46 33.4 1 8241. 1 0
#> 8 4 37 27.7 3 7282. 1 0
#> 9 4 37 29.8 2 6406. 0 0
#> 10 9 60 25.8 0 28923. 1 0
#> # ℹ 1,328 more rows
#> # ℹ 3 more variables: region_northeast <int>, region_northwest <int>,
#> # region_southwest <int>
Here we can see a summary of the binned cols. The with quintile
“value” bins we can see that the top 20% of charges comes from the top x
people.
insurance_bins %>%
bin_summary()
#> # A tibble: 10 × 14
#> column method n_bins .rank .min .mean .max .count .uniques
#> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <int>
#> 1 charges equal freq 10 10 34806. 42267. 63770. 136 136
#> 2 charges equal freq 10 9 20421. 26063. 34780. 130 130
#> 3 charges equal freq 10 8 13823. 16757. 20297. 135 135
#> 4 charges equal freq 10 7 11412. 12464. 13770. 134 134
#> 5 charges equal freq 10 6 9386. 10416. 11397. 134 134
#> 6 charges equal freq 10 5 7372. 8386. 9378. 134 134
#> 7 charges equal freq 10 4 5484. 6521. 7358. 134 134
#> 8 charges equal freq 10 3 3994. 4752. 5478. 133 133
#> 9 charges equal freq 10 2 2353. 3140. 3990. 134 134
#> 10 charges equal freq 10 1 1122. 1797. 2332. 134 133
#> # ℹ 5 more variables: relative_value <dbl>, .sum <dbl>, .med <dbl>, .sd <dbl>,
#> # width <dbl>
Create and analyze xgboost model for binary classification
Let’s see if we can predict whether a customer is a smoker.
Xgboost considers the first level of a factor to be the “positive
event” so let’s ensure that “yes” is the first level using
framecleaner::set_fct
insurance %>%
set_fct(smoker, first_level = "yes") -> insurance
insurance %>%
create_dummies(where(is.character), remove_first_dummy = T) -> insurance_dummies
insurance_dummies %>%
diagnose
#> # A tibble: 9 × 6
#> variables types missing_count missing_percent unique_count unique_rate
#> <chr> <chr> <int> <dbl> <int> <dbl>
#> 1 age integ… 0 0 47 0.0351
#> 2 bmi numer… 0 0 548 0.410
#> 3 children integ… 0 0 6 0.00448
#> 4 smoker factor 0 0 2 0.00149
#> 5 charges numer… 0 0 1337 0.999
#> 6 sex_male integ… 0 0 2 0.00149
#> 7 region_northwest integ… 0 0 2 0.00149
#> 8 region_southeast integ… 0 0 2 0.00149
#> 9 region_southwest integ… 0 0 2 0.00149
And create a new formula for the binary classification.
insurance_dummies %>%
tidy_formula(target = smoker) -> smoker_form
smoker_form
#> smoker ~ age + bmi + children + charges + sex_male + region_northwest +
#> region_southeast + region_southwest
#> <environment: 0x139c34188>
Now we can create a basic model using tidy_xgboost.
Built in heuristic will automatically recognize this as a binary
classification task. We can tweak some parameters to add some
regularization and increase the number of trees. Xgboost will
automatically output feature importance on the training set, and a
measure of accuracy tested on a validation set.
insurance_dummies %>%
tidy_xgboost(formula = smoker_form,
mtry = .5,
trees = 100L,
loss_reduction = 1,
alpha = .1,
sample_size = .7) -> smoker_xgb_classif
#> # A tibble: 15 × 3
#> .metric .estimate .formula
#> <chr> <dbl> <chr>
#> 1 accuracy 0.957 TP + TN / total
#> 2 kap 0.867 <NA>
#> 3 sens 0.91 TP / actually P
#> 4 spec 0.969 TN / actually N
#> 5 ppv 0.879 TP / predicted P
#> 6 npv 0.977 TN / predicted N
#> 7 mcc 0.868 <NA>
#> 8 j_index 0.879 <NA>
#> 9 bal_accuracy 0.939 sens + spec / 2
#> 10 detection_prevalence 0.206 predicted P / total
#> 11 precision 0.879 PPV, 1-FDR
#> 12 recall 0.91 sens, TPR
#> 13 f_meas 0.894 HM(ppv, sens)
#> 14 baseline_accuracy 0.801 majority class / total
#> 15 roc_auc 0.987 <NA>
Obtain predictions
smoker_xgb_classif %>%
tidy_predict(newdata = insurance_dummies, form = smoker_form) -> insurance_fit
Analyze the probabilities using tidybins. We can find the top 20% of
customers most likely to be smokers.
names(insurance_fit)[length(names(insurance_fit)) - 1] -> prob_preds
insurance_fit %>%
bin_cols(prob_preds, n_bins = 5) -> insurance_fit1
insurance_fit1 %>%
bin_summary()
#> # A tibble: 5 × 14
#> column method n_bins .rank .min .mean .max .count .uniques
#> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <int>
#> 1 smoker_preds_prob… equal… 5 5 0.643 0.921 0.987 268 263
#> 2 smoker_preds_prob… equal… 5 4 0.00620 0.0988 0.641 267 250
#> 3 smoker_preds_prob… equal… 5 3 0.00481 0.00527 0.00620 267 240
#> 4 smoker_preds_prob… equal… 5 2 0.00449 0.00462 0.00481 269 232
#> 5 smoker_preds_prob… equal… 5 1 0.00420 0.00436 0.00449 267 190
#> # ℹ 5 more variables: relative_value <dbl>, .sum <dbl>, .med <dbl>, .sd <dbl>,
#> # width <dbl>
Evaluate the training error. eval_preds uses both the
probability estimates and binary estimates to calculate a variety of
metrics.
insurance_fit1 %>%
eval_preds()
#> # A tibble: 4 × 5
#> .metric .estimator .estimate model target
#> <chr> <chr> <dbl> <chr> <chr>
#> 1 accuracy binary 0.987 smoker_xgb_classif smoker
#> 2 f_meas binary 0.970 smoker_xgb_classif smoker
#> 3 precision binary 0.945 smoker_xgb_classif smoker
#> 4 roc_auc binary 1.00 smoker_xgb_classif smoker
Traditional yardstick confusion matrix can be created
manually.
names(insurance_fit)[length(names(insurance_fit))] -> class_preds
insurance_fit1 %>%
yardstick::conf_mat(truth = smoker, estimate = class_preds) -> conf_mat_sm
conf_mat_sm
#> Truth
#> Prediction yes no
#> yes 273 16
#> no 1 1048