library(tidyverse)
library(tidymodels)
# Missingness packages
library(mice)
library(naniar, exclude = "n_complete")
source("https://github.com/jjcurtin/lab_support/blob/main/fun_eda.R?raw=true")
source("https://github.com/jjcurtin/lab_support/blob/main/fun_plots.R?raw=true")
options(conflicts.policy = "depends.ok")
path_data <- "Psych_752/Application Exams/Data"Airline Satisfaction Analysis
Setup Chunk
Loading necessary packages, setting conflict policy, and defining a path.
EDA
Data Loading & Cleaning
data <- read_csv(here::here(path_data, "airline_passenger_satisfaction.csv")) |>
janitor::clean_names() |>
glimpse()Rows: 600 Columns: 24
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): Gender, customer_type, type_of_travel, customer_class, satisfaction
dbl (19): id, age, flight_distance, inflight_wifi_service, departure_arrival...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 600
Columns: 24
$ id <dbl> 94031, 81130, 1431, 42411, 33339, 67…
$ gender <chr> "Male", "Male", "Female", NA, "Femal…
$ customer_type <chr> "Loyal Customer", "Loyal Customer", …
$ age <dbl> 32, 18, 43, 53, 25, 58, 51, 42, 44, …
$ type_of_travel <chr> "Business travel", "Business travel"…
$ customer_class <chr> "Business", "Business", "Business", …
$ flight_distance <dbl> 1389, 3456, 2899, 2694, 2419, 421, 3…
$ inflight_wifi_service <dbl> 2, 2, 0, 1, 5, 0, 5, 2, 1, 1, 2, 5, …
$ departure_arrival_time_convenient <dbl> 2, 2, 0, 1, 5, 3, 5, 2, 1, 1, 2, 5, …
$ ease_of_online_booking <dbl> 2, 2, 0, 1, 5, 0, 5, 3, 1, 1, 2, 5, …
$ gate_location <dbl> 2, 2, 5, 1, 5, 3, 5, 2, 1, 1, 2, 5, …
$ food_and_drink <dbl> 5, 4, 4, 3, 5, 3, 4, 2, 2, 5, 5, 3, …
$ online_boarding <dbl> 5, 4, 5, 5, 5, 3, 5, 4, 5, 5, 5, 5, …
$ seat_comfort <dbl> 5, 4, 5, 4, 4, 3, 4, 4, 4, 5, 4, 5, …
$ inflight_entertainment <dbl> 5, 4, 2, 4, 5, 4, 4, 4, 5, 5, 4, 5, …
$ onboard_service <dbl> 4, 4, 2, 4, 3, 4, 4, 4, 5, 3, 4, 5, …
$ leg_room_service <dbl> 5, 2, 2, 4, 3, 0, 4, 4, 5, 1, 4, 5, …
$ baggage_handling <dbl> 5, 1, 2, 4, 5, 4, 4, 4, 5, 3, 4, 5, …
$ checkin_service <dbl> 5, 3, 3, 3, 3, 3, 3, 4, 3, 2, 5, 5, …
$ inflight_service <dbl> 5, 3, 2, 4, 4, 4, 4, 4, 5, 4, 4, 5, …
$ cleanliness <dbl> 5, 4, 4, 3, 5, 3, 5, 3, 4, 5, 5, 5, …
$ departure_delay_in_minutes <dbl> 67, 0, 6, 0, 3, 0, 0, 50, 8, 21, 21,…
$ arrival_delay_in_minutes <dbl> 79, 0, 0, 0, 2, 0, 0, 42, 5, 4, 7, 1…
$ satisfaction <chr> "satisfied", "satisfied", "satisfied…data_clean <- data |>
mutate(
id = as.character(id)) |>
mutate(
gender = factor(gender,
levels = c("Male", "Female"),
labels = c("0", "1")),
customer_type = factor(customer_type,
levels = c("disloyal Customer", "Loyal Customer"),
labels = c("0", "1")),
type_of_travel = factor(type_of_travel,
levels = c("Business travel", "Personal Travel"),
labels = c("0", "1"))) |>
mutate(
customer_class = factor(customer_class,
levels = c("Eco", "Eco Plus", "Business"),
ordered = TRUE),
satisfaction = factor(satisfaction,
levels = c("neutral or dissatisfied", "satisfied"),
labels = c("0", "1"))) |>
glimpse()Rows: 600
Columns: 24
$ id <chr> "94031", "81130", "1431", "42411", "…
$ gender <fct> 0, 0, 1, NA, 1, 1, 0, NA, 0, 0, 1, 1…
$ customer_type <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ age <dbl> 32, 18, 43, 53, 25, 58, 51, 42, 44, …
$ type_of_travel <fct> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
$ customer_class <ord> Business, Business, Business, Busine…
$ flight_distance <dbl> 1389, 3456, 2899, 2694, 2419, 421, 3…
$ inflight_wifi_service <dbl> 2, 2, 0, 1, 5, 0, 5, 2, 1, 1, 2, 5, …
$ departure_arrival_time_convenient <dbl> 2, 2, 0, 1, 5, 3, 5, 2, 1, 1, 2, 5, …
$ ease_of_online_booking <dbl> 2, 2, 0, 1, 5, 0, 5, 3, 1, 1, 2, 5, …
$ gate_location <dbl> 2, 2, 5, 1, 5, 3, 5, 2, 1, 1, 2, 5, …
$ food_and_drink <dbl> 5, 4, 4, 3, 5, 3, 4, 2, 2, 5, 5, 3, …
$ online_boarding <dbl> 5, 4, 5, 5, 5, 3, 5, 4, 5, 5, 5, 5, …
$ seat_comfort <dbl> 5, 4, 5, 4, 4, 3, 4, 4, 4, 5, 4, 5, …
$ inflight_entertainment <dbl> 5, 4, 2, 4, 5, 4, 4, 4, 5, 5, 4, 5, …
$ onboard_service <dbl> 4, 4, 2, 4, 3, 4, 4, 4, 5, 3, 4, 5, …
$ leg_room_service <dbl> 5, 2, 2, 4, 3, 0, 4, 4, 5, 1, 4, 5, …
$ baggage_handling <dbl> 5, 1, 2, 4, 5, 4, 4, 4, 5, 3, 4, 5, …
$ checkin_service <dbl> 5, 3, 3, 3, 3, 3, 3, 4, 3, 2, 5, 5, …
$ inflight_service <dbl> 5, 3, 2, 4, 4, 4, 4, 4, 5, 4, 4, 5, …
$ cleanliness <dbl> 5, 4, 4, 3, 5, 3, 5, 3, 4, 5, 5, 5, …
$ departure_delay_in_minutes <dbl> 67, 0, 6, 0, 3, 0, 0, 50, 8, 21, 21,…
$ arrival_delay_in_minutes <dbl> 79, 0, 0, 0, 2, 0, 0, 42, 5, 4, 7, 1…
$ satisfaction <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …data_nona <- data_clean |>
drop_na()After loading in the data and cleaning the names, we needed to clean the dataset itself.
1) All binary predictors were coded as labeled factors (0 and 1) to keep modeling smooth. 2) Customer class was treated as an ordered factor to preserve meaningful ordering but without assuming any specific spacing between levels. 3) All Likert-type scale variables (0-5) were left as numerics to allow for flexible modeling. 4) Outcome variable (satisfaction) was coded 0 (not satisfied/neutral) and 1 (satisfied).
Then I made a new dataset from this version that undergoes casewise deletion of missing values/NAs, this left us with just over 2/3rds of the cases since many had at least 1 missing value. We’ll retain both versions for now.
Additional EDA
data_clean |>
skim_some()| Name | data_clean |
| Number of rows | 600 |
| Number of columns | 24 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 5 |
| numeric | 18 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| id | 0 | 1 | 3 | 6 | 0 | 600 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| gender | 117 | 0.80 | FALSE | 2 | 1: 262, 0: 221 |
| customer_type | 0 | 1.00 | FALSE | 2 | 1: 496, 0: 104 |
| type_of_travel | 0 | 1.00 | FALSE | 2 | 0: 480, 1: 120 |
| customer_class | 5 | 0.99 | TRUE | 3 | Bus: 359, Eco: 197, Eco: 39 |
| satisfaction | 0 | 1.00 | FALSE | 2 | 1: 400, 0: 200 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | p0 | p100 |
|---|---|---|---|---|
| age | 0 | 1.00 | 8 | 85 |
| flight_distance | 0 | 1.00 | 56 | 3989 |
| inflight_wifi_service | 0 | 1.00 | 0 | 5 |
| departure_arrival_time_convenient | 0 | 1.00 | 0 | 5 |
| ease_of_online_booking | 0 | 1.00 | 0 | 5 |
| gate_location | 0 | 1.00 | 1 | 5 |
| food_and_drink | 0 | 1.00 | 1 | 5 |
| online_boarding | 41 | 0.93 | 0 | 5 |
| seat_comfort | 0 | 1.00 | 1 | 5 |
| inflight_entertainment | 0 | 1.00 | 1 | 5 |
| onboard_service | 22 | 0.96 | 1 | 5 |
| leg_room_service | 0 | 1.00 | 0 | 5 |
| baggage_handling | 0 | 1.00 | 1 | 5 |
| checkin_service | 0 | 1.00 | 1 | 5 |
| inflight_service | 0 | 1.00 | 1 | 5 |
| cleanliness | 0 | 1.00 | 1 | 5 |
| departure_delay_in_minutes | 7 | 0.99 | 0 | 360 |
| arrival_delay_in_minutes | 3 | 1.00 | 0 | 409 |
# Quite a bit of missingness (mostly from gender), let's see if we can determine type of missingness and potentially move forward without simply excluding all cases with any missing.
d <- data_clean |> select(gender, customer_class, customer_type, age, flight_distance, online_boarding, onboard_service, departure_delay_in_minutes, arrival_delay_in_minutes)
skimr::skim(d)| Name | d |
| Number of rows | 600 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| factor | 3 |
| numeric | 6 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| gender | 117 | 0.80 | FALSE | 2 | 1: 262, 0: 221 |
| customer_class | 5 | 0.99 | TRUE | 3 | Bus: 359, Eco: 197, Eco: 39 |
| customer_type | 0 | 1.00 | FALSE | 2 | 1: 496, 0: 104 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1.00 | 39.97 | 14.08 | 8 | 28 | 40.5 | 51.0 | 85 | ▃▆▇▃▁ |
| flight_distance | 0 | 1.00 | 1318.55 | 1061.79 | 56 | 431 | 974.5 | 1998.5 | 3989 | ▇▅▂▂▂ |
| online_boarding | 41 | 0.93 | 3.56 | 1.33 | 0 | 3 | 4.0 | 5.0 | 5 | ▂▃▃▇▆ |
| onboard_service | 22 | 0.96 | 3.51 | 1.27 | 1 | 3 | 4.0 | 5.0 | 5 | ▂▃▆▇▇ |
| departure_delay_in_minutes | 7 | 0.99 | 13.82 | 37.28 | 0 | 0 | 0.0 | 9.0 | 360 | ▇▁▁▁▁ |
| arrival_delay_in_minutes | 3 | 1.00 | 13.51 | 38.21 | 0 | 0 | 0.0 | 8.0 | 409 | ▇▁▁▁▁ |
gg_miss_upset(d)Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
ℹ The deprecated feature was likely used in the UpSetR package.
Please report the issue to the authors.Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the UpSetR package.
Please report the issue to the authors.Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
ℹ The deprecated feature was likely used in the UpSetR package.
Please report the issue to the authors.
Most missingness comes from gender, 105 of the 181 variables that are missing anything are only missing gender. 32 are missing only online boarding, 18 are missing only onboard_service, and very few (13) participants had missingness on more than 1 variable, with none having more than 2. Would be great to be able to retain gender by imputing realistic “values” using correlated variables, then we could compare performance of our eventual models trained on each version and see if imputation helps our model performance.
Missingness EDA continued
d |>
mutate(across(where(is.factor), ~ as.numeric(.))) |>
cor(use = "pairwise.complete.obs") |>
corrplot::corrplot.mixed()
Seems gender is not associated with any variables, so we’ll probably just have to drop all cases missing gender. Other variables don’t have much missingness so won’t try to impute for those either. We’ll use data_nona from here on out as we also have no reason to assume any coding errors for the NAs.
Univariate EDA & Class Balance
data_nona |>
select(where(is.numeric)) |>
summary() age flight_distance inflight_wifi_service
Min. : 8.00 Min. : 67.0 Min. :0.000
1st Qu.:29.00 1st Qu.: 417.5 1st Qu.:2.000
Median :41.00 Median : 946.0 Median :3.000
Mean :40.05 Mean :1342.1 Mean :2.907
3rd Qu.:51.00 3rd Qu.:2111.0 3rd Qu.:4.000
Max. :76.00 Max. :3989.0 Max. :5.000
departure_arrival_time_convenient ease_of_online_booking gate_location
Min. :0.000 Min. :0.000 Min. :1.00
1st Qu.:2.000 1st Qu.:2.000 1st Qu.:2.00
Median :3.000 Median :3.000 Median :3.00
Mean :3.002 Mean :2.866 Mean :3.01
3rd Qu.:4.000 3rd Qu.:4.000 3rd Qu.:4.00
Max. :5.000 Max. :5.000 Max. :5.00
food_and_drink online_boarding seat_comfort inflight_entertainment
Min. :1.000 Min. :0.00 Min. :1.000 Min. :1.000
1st Qu.:2.000 1st Qu.:3.00 1st Qu.:3.000 1st Qu.:3.000
Median :3.000 Median :4.00 Median :4.000 Median :4.000
Mean :3.301 Mean :3.58 Mean :3.644 Mean :3.578
3rd Qu.:4.000 3rd Qu.:5.00 3rd Qu.:5.000 3rd Qu.:5.000
Max. :5.000 Max. :5.00 Max. :5.000 Max. :5.000
onboard_service leg_room_service baggage_handling checkin_service
Min. :1.000 Min. :0.000 Min. :1.000 Min. :1.000
1st Qu.:3.000 1st Qu.:3.000 1st Qu.:3.000 1st Qu.:3.000
Median :4.000 Median :4.000 Median :4.000 Median :3.000
Mean :3.547 Mean :3.585 Mean :3.768 Mean :3.332
3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:4.000
Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
inflight_service cleanliness departure_delay_in_minutes
Min. :1.000 Min. :1.000 Min. : 0.00
1st Qu.:3.000 1st Qu.:3.000 1st Qu.: 0.00
Median :4.000 Median :3.000 Median : 0.00
Mean :3.831 Mean :3.322 Mean : 13.97
3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.: 8.00
Max. :5.000 Max. :5.000 Max. :360.00
arrival_delay_in_minutes
Min. : 0.00
1st Qu.: 0.00
Median : 0.00
Mean : 13.45
3rd Qu.: 6.50
Max. :409.00 data_nona |>
pivot_longer(where(is.numeric)) |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30) +
facet_wrap(~ name, scales = "free")
data_nona |> count(satisfaction)# A tibble: 2 × 2
satisfaction n
<fct> <int>
1 0 138
2 1 281Ranges all look reasonable, we are seeing that we do have imbalanced classes as 67% of respondents were satisfied, so we’ll want to use something besides raw accuracy as our performance metric (F1, balanced accuracy, kappa, or auROC all considered at this point). In our specific context here, we are looking to predict passenger satisfaction. As such, we may want to minimize false positives (we predict satisfied for customers that aren’t satisfied) as we might want to reach out to potentially unsatisfied customers to improve customer relations. We will later consider increasing our classification threshold so that we miss fewer unsatisfied customers, focusing on a higher specificity. We can simply use auROC to tune, but then evaluate based on auROC and specificity. Ideally we’d speak with the rest of the team to determine beforehand what our optimizing and satisficing metrics are, along with what value to use for our satisficing metric, but unfortunately everyone else is out of office we are unable to predetermine these.
Bivariate EDA
data_nona |>
select(where(is.numeric)) |>
cor(use = "pairwise.complete.obs") |>
corrplot::corrplot.mixed()
data_nona <- data_nona |>
mutate(delay_max = pmax(departure_delay_in_minutes, arrival_delay_in_minutes))Departure delay and arrival delay have correlation of .98, we’re going to make a new variable called delay_max that simply takes the maximum delay value of either departure or arrival for each participant. This seems preferable to taking the mean of the two as we expect that satisfaction is likely to be largely driven by the worst delay-related experience. Keeping both provides very little additional information and would mean we have high multicollinearity.
Fitting & Evaluating models
Splitting data
set.seed(123)
split <- initial_split(data_nona, prop = 0.8, strata = satisfaction)
train_data <- training(split)
test_data <- testing(split)We’ll want a performance estimate of our model on held-out data, so we split the data here, will only touch the training data. Since we have 419 observations it seems setting aside 20% as held-out will be just fine here. Making our training set any smaller may inflate bias and making the held-out set smaller might give is a high variance estimate, to combat this we will additionally perform repeated kfold cross-validation on the training set, using that final held-out 20% as a pure test set/final check. We considered a single validation or nested cross-validation approach as well, but due to high variance and unnecessary complexity we decided against them. We’ll do 5 folds as we don’t want to do too many since having the sample sizes within each split be any smaller may cause variance problems, and we’ll do 10 repeats for additional variance reduction, but not any more for computational cost/time considerations.
Cross-validation setup/folds
folds <- vfold_cv(train_data, v = 5, repeats = 10, strata = satisfaction)Recipe
rec <- recipe(satisfaction ~ ., data = train_data) |>
step_rm(c("id", "departure_delay_in_minutes", "arrival_delay_in_minutes")) |>
step_dummy(all_nominal_predictors()) |>
step_zv(all_predictors()) |>
step_normalize(all_numeric_predictors())Here we make our recipe for our model. We’ll be training on the held-in training data, but first we’ve removed id (no predictive information) and also both highly redundant delay metrics as we have combined them into a single feature. We then dummy coded numerics so they can be used and removed all zero variance predictors as they also contain no useful information. Finally we normalized the numerics because we’ll be tuning a glmnet model in which large coefficients are penalized, so they all must be on the same scale to avoid unfairly penalizing coefficients for flight distance, age, etc.
Model Specification & Workflow
m1_log_reg <- logistic_reg(
penalty = tune(),
mixture = tune()) |>
set_engine("glmnet")
wf <- workflow() |>
add_recipe(rec) |>
add_model(m1_log_reg)We’re choosing to tune a glmnet as it gives us a simple, interpretable, and semi-flexible classification model (flexible in terms of regularization and potentially feature selection, though still assumes linear DGP). If this linear classsification baseline performs poorly, we’ll likely transition to a nonparametric approach. We then combine the preprocessing recipe and model specification into a workflow so that the preprocessing steps are concisely and consistently applied within our cross-validation folds.
Hyperparameter Grid
grid <- grid_regular(
penalty(range = c(-3, -1)),
mixture(),
levels = 5)We define a tuning grid for the penalty and mixture hyperparameters. We started with a more coarse grid before refining our search for optimal penalty values, we find that small to moderate penalties are ideal. This left us with -3 to -1 (log scale) as our penalty consideration range. The mixture simply tries different blends of L2 and L1 regularization, along with pure L2 and pure L1.
Tuning
# About 2 minutes
tuned_results <- tune_grid(
wf,
resamples = folds,
grid = grid,
metrics = metric_set(roc_auc, specificity, accuracy))Results
autoplot(tuned_results)
tuned_results |>
collect_metrics() |>
select(penalty, mixture, .metric, mean) |>
pivot_wider(names_from = .metric, values_from = mean) |>
arrange(desc(roc_auc)) |>
slice_head(n = 20)# A tibble: 20 × 5
penalty mixture accuracy roc_auc specificity
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0316 1 0.759 0.867 0.880
2 0.0316 0.75 0.764 0.866 0.877
3 0.1 0.5 0.750 0.865 0.902
4 0.0316 0.5 0.770 0.864 0.875
5 0.1 0.25 0.751 0.863 0.883
6 0.0316 0.25 0.780 0.861 0.871
7 0.1 0.75 0.754 0.861 0.922
8 0.01 1 0.785 0.861 0.872
9 0.1 1 0.756 0.861 0.940
10 0.01 0.75 0.786 0.860 0.871
11 0.01 0.5 0.788 0.858 0.869
12 0.0316 0 0.785 0.858 0.867
13 0.1 0 0.777 0.857 0.874
14 0.01 0.25 0.791 0.857 0.867
15 0.001 0 0.785 0.857 0.867
16 0.00316 0 0.785 0.857 0.867
17 0.01 0 0.785 0.857 0.867
18 0.00316 1 0.790 0.855 0.866
19 0.00316 0.75 0.788 0.854 0.864
20 0.00316 0.5 0.788 0.854 0.863Here we have our results, we’ve printed the top 20 in regard to auroc because we are also very interesting in specificity. Looking through these we see relatively stable auroc in that there aren’t big gaps in performance, but there are quite big gaps in specificity. Model 7 specifically stuck out to us as it is performing only .005 units lower on auroc, but has an extremely good specificity of .942, which is .06 higher than the top performing model in respect to auroc. In other words, we got a substantial gain in specificity with only a minor reduction in auroc. This model also happens to be pure L1 regularized, so we’ll have an even simpler model, and this will allow us to require less data from customers and still be able to use the model. If we aren’t as interested in specificity as I have outlined, we can easily come back and select a different model irrespective of specificity, but for now this will be our model to continue with.
Further, though the model’s overall accuracy (≈0.755) is slightly lower than the maximum achievable, this reflects a deliberate tradeoff aligned with our objective here. We prioritized correctly identifying dissatisfied customers, which leads us to favor models with higher specificity even at the cost of some accuracy. This results in a more conservative classification rule where some satisfied customers may be flagged for follow-up. However, the cost of unnecessarily contacting a satisfied customer is relatively low and still may be beneficial regarding costumer relations, whereas failing to identify a dissatisfied customer may lead to lost future business. As such, the chosen model sacrifices a small amount of overall accuracy in order to better capture customers at risk of dissatisfaction, which is more valuable in this context.
Fit our best model to the full training data and evaluate on held-out data
best_params <- tibble(
penalty = 0.1,
mixture = 1)
final_wf <- finalize_workflow(wf, best_params)
final_model <- fit(final_wf, data = train_data)
preds <- predict(final_model, test_data, type = "prob") |>
bind_cols(predict(final_model, test_data, type = "class")) |>
bind_cols(test_data)Held-out performance
roc_auc(preds, truth = satisfaction, .pred_0)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.846spec(preds, truth = satisfaction, estimate = .pred_class)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 spec binary 0.930Evaluating our model on test we see that there was some slight optimization bias in our earlier performance estimate, but only marginally so. We still had an extremely good specificity of .93 and an acceptable auROC of .846. Our model seems to be quite low variance if we compare these to the earlier values of .94 and .86. Our bias is still quite low, though I’d make the argument that we don’t mind it being where it is as there’s such a low cost to reaching out to additional customers. However, if we are providing compensation or travel vouchers based upon the model, this would change quite a bit. Depending on what we want we could make the necessary changes, but this appears to serve as a fantastic baseline model.
fullglmnet_fit <- final_model$fit$fit$fit
coefficients <- coef(fullglmnet_fit, s = 0.1)
coefficients <- as.matrix(coefficients)
coefficients <- data.frame(
term = rownames(coefficients),
estimate = coefficients[,1])
coefficients term
(Intercept) (Intercept)
age age
flight_distance flight_distance
inflight_wifi_service inflight_wifi_service
departure_arrival_time_convenient departure_arrival_time_convenient
ease_of_online_booking ease_of_online_booking
gate_location gate_location
food_and_drink food_and_drink
online_boarding online_boarding
seat_comfort seat_comfort
inflight_entertainment inflight_entertainment
onboard_service onboard_service
leg_room_service leg_room_service
baggage_handling baggage_handling
checkin_service checkin_service
inflight_service inflight_service
cleanliness cleanliness
delay_max delay_max
gender_X1 gender_X1
customer_type_X1 customer_type_X1
type_of_travel_X1 type_of_travel_X1
customer_class_1 customer_class_1
customer_class_2 customer_class_2
estimate
(Intercept) 0.759336701
age 0.000000000
flight_distance 0.006232991
inflight_wifi_service 0.000000000
departure_arrival_time_convenient 0.000000000
ease_of_online_booking 0.000000000
gate_location 0.000000000
food_and_drink 0.000000000
online_boarding 0.260454027
seat_comfort 0.000000000
inflight_entertainment 0.095860015
onboard_service 0.000000000
leg_room_service 0.000000000
baggage_handling 0.000000000
checkin_service 0.000000000
inflight_service 0.000000000
cleanliness 0.000000000
delay_max 0.000000000
gender_X1 0.000000000
customer_type_X1 0.000000000
type_of_travel_X1 -0.426792723
customer_class_1 0.065387423
customer_class_2 0.000000000coefficients_used <- coefficients[coefficients$estimate != 0, ]
coefficients_used term estimate
(Intercept) (Intercept) 0.759336701
flight_distance flight_distance 0.006232991
online_boarding online_boarding 0.260454027
inflight_entertainment inflight_entertainment 0.095860015
type_of_travel_X1 type_of_travel_X1 -0.426792723
customer_class_1 customer_class_1 0.065387423Importantly, we looked into our selected model to see what information it’s using and what we found may be extremely interesting/valuable to us. The model only retained flight distance, online boarding, inflight entertainment, type of travel, and the dummy variable coding for Economy vs Economy plus. What this means is that if we move forward with this model, these are the only datapoints we’ll need to collect and all but the inflight entertainment rating are things we are likely collecting for all passengers already. Therefore, implementing this model means there is no need for collecting the other data if doing so is expensive or time-consuming, and simply asking for a rating of inflight entertainment would suffice. Further, we may look into additionally prioritizing inflight entertainment, though we should be careful in taking away resources from any other aspects of the experience, especially without looking deeper into specific feature importance metrics like shapley values.
Prediction summary:
Spending your data:
We split the dataset into an 80/20 train–test split using stratification on the outcome to preserve class balance. The training set was used for all model fitting and model selection, while the test set was held out as a final, unbiased estimate of model performance. Within the training data, we used repeated 5-fold cross-validation (10 repeats) to evaluate model configurations. This approach balances bias and variance trade-off considerations. Holding out a test set prevents optimization bias in performance estimates, while repeated cross-validation reduces the variance of model evaluation compared to a single split. Given the moderate sample size (n = 419 after casewise deletion), this strategy allows us to efficiently use the data while still maintaining a reliable estimate of generalization performance. The main tradeoff is that reserving 20% of the data reduces the amount available for model fitting, potentially increasing bias slightly. However, this is mitigated by cross-validation, and the benefit of having a clean, unbiased test evaluation outweighs this cost.
Model configurations:
We focused on a regularized logistic regression model (glmnet), tuning two key hyperparameters: the penalty, which controls the strength of regularization, and the mixture parameter, which determines the balance between ridge (L2) and lasso (L1) penalties. We explored a range of penalty values (on a log scale) and multiple mixture values spanning ridge to lasso, beginning with a broader grid and then refining to smaller penalty values where performance was strongest. This ensured that we considered models ranging from L2 regularization with little to no penalty all the way to large penalties and L1 regularization/feature selection. We determined that we had explored a sufficient range when performance was less variable (in terms of auROC, specificity, and accuracy) across a region of hyperparameters, indicating that further expansion of the grid was unlikely to yield meaningful improvements. Earlier grid values showed sporatic and massive decreases in these metrics. Additionally, comparing ridge, lasso, and elastic net configurations ensured that both shrinkage-based and feature-selection-based models were considered.
Model selection:
Model performance during tuning was evaluated primarily using auROC, specificity, and accuracy. We didn’t overtly set values for any of these values beforehand which would’ve been ideal. If we had, we could have set either auROC or specificity as our satisficing metric and the other as our optimization metric. In practice, we basically did this post-hoc by focusing only on models with an acceptable auROC , before selecting the model out of those that maximized specificity. auROC was chosen because it is threshold-independent and provides a robust measure of the model’s ability to discriminate between satisfied and unsatisfied customers, particularly in the presence of class imbalance (67/33). Specificity was a big part of the model decision-making process as we supposed that minimizing false positives (incorrectly classifying dissatisfied customers as satisfied) was worthwhile given the potential uses for the model and the ramifications, as we explained in more depth previously. While the model with the highest auROC achieved slightly better overall discrimination, an alternative model was selected that sacrificed a small amount of auROC (≈0.005) in exchange for a substantial increase in specificity (≈0.06). This reflected a deliberate, cost-sensitive decision where maximizing overall discrimination was less important than correctly identifying dissatisfied customers.
Model characterization
The expected model performance for new data was approximated by evaluating on the held-out test set. The model achieved an auROC of 0.846 and a specificity of 0.93, indicating strong predictive discrimination and a high rate of correctly identifying dissatisfied customers. Comprehensively characterizing performance depends on our ultimate goals with the model, but we give a comprehensive performance characterization through auROC which measures overall discrimination, specificity which measures the models ability to correctly identify dissatisfied customers (potentially highly relevant to us), and accuracy which tells us the general correctness but is less informative under imbalance, though highly interpretable. In real-world terms, this model is intentionally conservative, it prioritizes identifying dissatisfied customers, even if that means sometimes flagging satisfied customers for follow-up. The validity of this choice depends entirely on what our goals with the model are, if we are simply contacting potentially dissatisfied customers for feedback/questions versus if we are offering compensation or vouchers based on predictions is paramount. The cost of contacting a satisfied customer is very low (and may even improve customer relations), whereas failing to identify a dissatisfied customer may result in lost future business. For this baseline model, we’ve assumed simple-contacting as the goal, which is reflected by the models asymmetric error structure (false positives as low cost but false negatives as potentially high cost), but this can be easily changed in future models if our objectives change.