library(tidyverse)library(tidymodels)options(conflicts.policy = "depends.ok")
path_data <- "Psych_752/Application Exams/Data"Tips Analysis
Setup Chunk
Reading in data
data <- read_csv(here::here(path_data, "tips.csv")) |>
glimpse()Rows: 304 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (7): customer_sex, smoker, day, time, server_sex, any_children, ordered_...
dbl (5): total_bill, group_size, tip_percentage, customer_age, number_drinks
ℹ 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: 304
Columns: 12
$ total_bill <dbl> 7.25, 9.60, 3.07, 11.61, 23.17, 14.31, 7.51, 16.32, 13…
$ customer_sex <chr> "Male", "Female", "Female", "Male", "Male", "Female", …
$ smoker <chr> "Yes", "Yes", "Yes", "No", "Yes", "Yes", "No", "Yes", …
$ day <chr> "Sun", "Sun", "Sat", "Sat", "Sun", "Sat", "Thur", "Fri…
$ time <chr> "Dinner", "Dinner", "Dinner", "Dinner", "Dinner", "Din…
$ group_size <dbl> 2, 2, 1, 2, 4, 2, 2, 2, 2, 2, 2, 5, 2, 2, 2, 3, 2, 3, …
$ tip_percentage <dbl> 71.03448, 41.66667, 32.57329, 29.19897, 28.05352, 27.9…
$ server_sex <chr> "Female", "Male", "Male", "Male", "Female", "Male", "M…
$ customer_age <dbl> 63.96420, 55.12560, 45.32490, 47.81858, 43.99298, 53.9…
$ any_children <chr> "Yes", "No", "No", "No", "No", "No", "No", "No", "Yes"…
$ ordered_dessert <chr> "No", "No", "No", "No", "No", "No", "Yes", "No", "Yes"…
$ number_drinks <dbl> 0, 0, 0, 2, 0, 0, 0, 5, 1, 0, 1, 0, 0, 4, 0, 0, 0, 0, …Simple data reclassing
data <- data |>
mutate(
customer_sex = factor(if_else(customer_sex == "Male", 1, 0),
levels = c(0, 1),
labels = c("Female", "Male")),
smoker = factor(if_else(smoker == "Yes", 1, 0),
levels = c(0, 1),
labels = c("No", "Yes")),
time = factor(if_else(time == "Dinner", 1, 0),
levels = c(0, 1),
labels = c("Lunch", "Dinner")),
server_sex = factor(if_else(server_sex == "Male", 1, 0),
levels = c(0, 1),
labels = c("Female", "Male")),
any_children = factor(if_else(any_children == "Yes", 1, 0),
levels = c(0, 1),
labels = c("No", "Yes")),
ordered_dessert = factor(if_else(ordered_dessert == "Yes", 1, 0),
levels = c(0, 1),
labels = c("No", "Yes")),
day = factor(day, ordered = FALSE))Direction
We want to zero in on the effect of age on tip percentage, not simply determine the most predictive model. Let’s first look at the simple effect of age on tip percentage and determine if there seems to be a linear relationship.
data |>
ggplot(aes(x = customer_age, y = tip_percentage)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = TRUE, color = "blue") +
geom_smooth(method = "lm", se = FALSE, color = "red", linetype = "dashed") +
labs(
x = "Customer Age",
y = "Tip Percentage",
title = "Relationship Between Age and Tip Percentage") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Definitely doesn’t appear to have a linear relationship at first glance, looks more quadratic, though we want to be careful about interpreting noise. Similarly, if we were to ignore a few outliers, it seems as if the relationship between age and tip percentage could actually be quite linear. When identifying the best model configuration we’ll want to consider each of these as a possibility.
Next, we want to consider potential covariates and interactions. Let’s first check univariate distributions and then we’ll check the relationship between age and our other variables.
data |>
pivot_longer(where(is.numeric)) |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30) +
facet_wrap(~ name, scales = "free")
Total bill is heavily positively skewed, we’ll add a log-transformed version to correct for this. We also know we’ll want to look at potential quadratic effects of age on tip percentage so we’ll make a squared version of age here as well. We’ll no longer consider number of drinks as a covariate given the extreme pooling into a singular response category.
data <- data |>
mutate(age_sq = customer_age^2,
log_total_bill = log(total_bill))Now lets look at numeric relationships between age and other variables.
data |>
select(where(is.numeric)) |>
cor(use = "pairwise.complete.obs") |>
corrplot::corrplot.mixed()
data |>
pivot_longer(cols = c(customer_sex, smoker, time, server_sex, any_children, ordered_dessert, day),
names_to = "variable",
values_to = "value") |>
ggplot(aes(x = value, y = customer_age)) +
geom_boxplot() +
facet_wrap(~ variable, scales = "free_x") +
theme_minimal()
We see the largest correlations for age are with tip percentage (.24), total bill (-.16), and the log of the total bill (-.15).
Let’s now look into the categorical variables. If there is any notable relationship with tip percentage then they will be considered as covariate candidates.
data |>
pivot_longer(cols = c(customer_sex, smoker, time, server_sex, any_children, ordered_dessert, day),
names_to = "variable",
values_to = "value") |>
ggplot(aes(x = value, y = tip_percentage)) +
geom_boxplot() +
facet_wrap(~ variable, scales = "free_x") +
theme_minimal()
There are some interesting distributions within the day variable as Tuesday and Wednesday seem to be associated with lower tips on average to a considerable degree. Any children, customer sex, dessert, and server sex seem to have a considerable amount of overlap between class distributions, but smoker and time look to have potentially meaningful shifts such that smokers tip less on average and lunch gets higher tips while dinner has many outliers on both sides of the distribution. Let’s check their relationships with our focal variable, age.
data |>
pivot_longer(cols = c(customer_sex, smoker, time, server_sex, any_children, ordered_dessert, day),
names_to = "variable",
values_to = "value") |>
ggplot(aes(x = value, y = customer_age)) +
geom_boxplot() +
facet_wrap(~ variable, scales = "free_x") +
theme_minimal()
Doesn’t appear that any groups/categories have major systematic differences in age. So these variables may be less likely to help (as controls/covariates) to isolate the age effect. The previously mentioned variables appeared to be stronger candidates for inclusion as covariates.
data |>
ggplot(aes(x = day, y = tip_percentage)) +
geom_boxplot() +
theme_minimal()
data |>
ggplot(aes(x = customer_age, y = tip_percentage, color = day)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = FALSE) +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
data |>
ggplot(aes(x = day, y = tip_percentage)) +
geom_boxplot(outlier.shape = NA) + # hide duplicate outliers
geom_jitter(aes(color = customer_age), width = 0.2, alpha = 0.7) +
scale_color_viridis_c() +
labs(
x = "Day",
y = "Tip Percentage",
color = "Customer Age",
title = "Tip Percentage by Day with Age Coloring") +
theme_minimal()
Looking into day more leads us towards including it as a covariate, doesn’t seem to be major differences in age across days, but definitely appears that certain days are reliably associated with lower tips.
Now we need to consider which statistical algorithm might best capture the relationship between customer age and tip percentage. We’ll definitely want something interpretable since we want to capture the age effect and not just optimize for prediction while obscuring the relationship. This leads us to conclude that we’ll stick with general linear models in which we can still include a quadratic term. Through our investigations we’ve determined we want to focus in on the following variables in our search for the optimal model configuration: Raw age, age as a quadratic, total bill, log-transformed total bill, day (potentially looking specifically at Tuesday and Wednesday vs others to capture small-tip days vs big-tip days), smoker, and time.
Before we move to model configurations, we’ll look at outliers.
data |>
ggplot(aes(x = tip_percentage)) +
geom_histogram(bins = 30) +
theme_minimal() +
labs(title = "Distribution of Tip Percentage")
While we do see a handful of extreme outliers, we have no reason to believe these are invalid/incorrect responses. Further, it is common to have major outliers on tipping percentage metrics as that is an inherent part of tipping behavior/culture; it is commonplace for people who have an exceptional experience to tip much, much larger than average. It is also plenty typical for customers with smaller checks to tip much larger proportionally, so removing outliers here would artificially bias our estimates away from true behavior towards average behavior. Therefore we won’t manipulate or exclude any outliers in our analysis.
Using what we’ve found so far, we’ll now compare across many different model configurations to determine which one best generalizes to new data. Since we need to do this many times, we’ll use repeated k-fold cross-validation to get stable estimates of held-out performance for each configuration.
set.seed(123)
folds <- vfold_cv(data, v = 5, repeats = 5)Model configurations need to address the following, is the age effect better captured by:
Considering quadratic form of age
Log-transformed total bill as covariate
Day as covariate or interaction variable
Smoker as covariate or interaction variable
Time as covariate
We’ll fit a compact model with raw age as our baseline. Then we’ll have a series of models where we have age as a quadratic and we’ll add in covariates. We expect to see that error will decrease from the compact to the quadratic, and then we can compare the quadratic model to those with differing covariates.
Since we want to compare across many model configurations, we’ve constructed a pipeline that simplifies this process. We’ll specify each of the model configurations we want to compare and add them into a list, and then each one will be trained on 4 folds of the data before evaluating on a 5th, held-out fold. Our function does this for each configuration and then computes the rmse of each configuration (25 times since 5-fold repeated 5 times), which we can average together for the models and then directly compare across each configuration.
models <- list(
m1 = tip_percentage ~ customer_age,
m2 = tip_percentage ~ customer_age + age_sq,
m3 = tip_percentage ~ customer_age + age_sq + total_bill,
m4 = tip_percentage ~ customer_age + age_sq + log_total_bill,
m5 = tip_percentage ~ customer_age + age_sq + day,
m6 = tip_percentage ~ customer_age + age_sq + smoker)get_rmse <- function(split, formula) {
train <- analysis(split)
test <- assessment(split)
model <- lm(formula, data = train)
preds <- predict(model, newdata = test)
rmse_vec(test$tip_percentage, preds)}
results <- map_dfr(names(models), function(name) {
formula <- models[[name]]
rmse_vals <- map_dbl(folds$splits, get_rmse, formula)
tibble(
model = name,
rmse = mean(rmse_vals))
})
results |> arrange(rmse)# A tibble: 6 × 2
model rmse
<chr> <dbl>
1 m5 6.71
2 m4 6.81
3 m3 6.88
4 m1 7.08
5 m2 7.11
6 m6 7.13Surprisingly, including age as a quadratic did not decrease model error compared to the compact model. So we’ll have to rerun some comparisons in which we use raw age with each of the covariates.
models <- list(
m1 = tip_percentage ~ customer_age,
m2 = tip_percentage ~ customer_age + log_total_bill,
m3 = tip_percentage ~ customer_age + day,
m4 = tip_percentage ~ customer_age + smoker,
m5 = tip_percentage ~ customer_age + log_total_bill + day)
results <- map_dfr(names(models), function(name) {
formula <- models[[name]]
rmse_vals <- map_dbl(folds$splits, get_rmse, formula)
tibble(
model = name,
rmse = mean(rmse_vals))})
results |> arrange(rmse)# A tibble: 5 × 2
model rmse
<chr> <dbl>
1 m5 6.31
2 m3 6.61
3 m2 6.80
4 m1 7.08
5 m4 7.09We see here that some of our hypotheses were supported, including log-transformed total bill and day as single covariates decreased error, but including both of them reduced error further, so they account for different variance. Including smoker as a covariate did not help. Now we can look at interactive models:
models <- list(
m1 = tip_percentage ~ customer_age * log_total_bill + day,
m2 = tip_percentage ~ customer_age * day + log_total_bill,
m3 = tip_percentage ~ customer_age * customer_sex + log_total_bill + day,
m4 = tip_percentage ~ customer_age * smoker + log_total_bill + day)
results <- map_dfr(names(models), function(name) {
formula <- models[[name]]
rmse_vals <- map_dbl(folds$splits, get_rmse, formula)
tibble(
model = name,
rmse = mean(rmse_vals))})
results |> arrange(rmse)# A tibble: 4 × 2
model rmse
<chr> <dbl>
1 m4 6.26
2 m3 6.37
3 m1 6.49
4 m2 6.56Interestingly, the best model still included log-transformed total bill and day as covariates, but it also included the age X smoker interaction. It seems entirely plausible that the effect of age on tipping differs for smokers and non-smokers, potentially with smoking being some sort of behavioral proxy (mood differences, appetite, impulsivity, alcohol consumption, etc.).
Before we identify this model as our best configuration, I want to verify that age as a quadratic truly doesn’t better capture the DGP.
models2 <- list(
m1 = tip_percentage ~ customer_age * smoker + log_total_bill + day,
m2 = tip_percentage ~ (customer_age + age_sq) * smoker + log_total_bill + day)results2 <- map_dfr(names(models2), function(name) {
formula <- models2[[name]]
rmse_vals <- map_dbl(folds$splits, get_rmse, formula)
tibble(
model = name,
rmse = mean(rmse_vals))
})
results2 |> arrange(rmse)# A tibble: 2 × 2
model rmse
<chr> <dbl>
1 m1 6.26
2 m2 6.44This is further evidence that quadratic age doesn’t improve generalization, so we can now move forward to fitting what we’ve identified as the best model configuration from our cross-validation analyses:
\(TipPercentage = b_0 + b_1(CustomerAge) + b_2(Smoker) + b_3(CustomerAge * Smoker) + b_4(LogTotalBill) + b_5(day) + e\)
best_model <- lm(tip_percentage ~ customer_age * smoker + log_total_bill + day,
data = data)
summary(best_model)
Call:
lm(formula = tip_percentage ~ customer_age * smoker + log_total_bill +
day, data = data)
Residuals:
Min 1Q Median 3Q Max
-11.235 -3.298 -0.569 2.332 52.940
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.32998 4.84798 5.637 4.05e-08 ***
customer_age 0.02669 0.08659 0.308 0.758088
smokerYes -27.75948 6.15682 -4.509 9.43e-06 ***
log_total_bill -4.44772 0.85898 -5.178 4.17e-07 ***
daySat 0.07995 1.61053 0.050 0.960443
daySun 1.34829 1.66801 0.808 0.419561
dayThur -0.05215 1.68383 -0.031 0.975312
dayTue -4.67319 1.81427 -2.576 0.010489 *
dayWed -7.51937 1.90793 -3.941 0.000101 ***
customer_age:smokerYes 0.62239 0.13786 4.515 9.19e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.224 on 294 degrees of freedom
Multiple R-squared: 0.3135, Adjusted R-squared: 0.2925
F-statistic: 14.92 on 9 and 294 DF, p-value: < 2.2e-16Before concluding, we want to also tune a linear regression model on all features while using L1 regularization. This final check is a wholly data-driven feature selection approach to use simply as a cross-reference with our explanatory-approach to finding an optimal model configuration. If we see extremely different feature selection and much lower RMSE from the resulting configurations, we will likely have to consider reassessing our earlier approch. We will similarly use repeated 5-fold cross validation here, dummy-coding nominal predictors and normalizing numerics since we must apply the penalty uniformly between features.
l1_rec <- recipe(
tip_percentage ~ customer_age + age_sq + total_bill + log_total_bill +
customer_sex + smoker + day + time + server_sex + any_children +
ordered_dessert + group_size + number_drinks, data = data) |>
step_dummy(all_nominal_predictors()) |>
step_normalize(all_numeric_predictors())
# lasso specification
lasso_mod <- linear_reg(
penalty = tune(),
mixture = 1) |>
set_engine("glmnet")
# workflow
lasso_wf <- workflow() |>
add_recipe(l1_rec) |>
add_model(lasso_mod)
# tuning grid
lambda_grid <- grid_regular(
penalty(range = c(-5, 1)),
levels = 50)
# tune with repeated CV
lasso_results <- tune_grid(
lasso_wf,
resamples = folds,
grid = lambda_grid,
metrics = metric_set(rmse))
# best lambda by RMSE
show_best(lasso_results, metric = "rmse")# A tibble: 5 × 7
penalty .metric .estimator mean n std_err .config
<dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
1 0.339 rmse standard 6.34 25 0.448 pre0_mod38_post0
2 0.256 rmse standard 6.35 25 0.443 pre0_mod37_post0
3 0.193 rmse standard 6.35 25 0.438 pre0_mod36_post0
4 0.146 rmse standard 6.36 25 0.432 pre0_mod35_post0
5 0.450 rmse standard 6.37 25 0.449 pre0_mod39_post0# select best model
best_lambda <- select_best(lasso_results, metric = "rmse")
# finalize workflow
final_lasso_wf <- finalize_workflow(
lasso_wf,
best_lambda)
# fit on full dataset
final_lasso_fit <- fit(final_lasso_wf, data = data)
coef(final_lasso_fit$fit$fit$fit,s = best_lambda$penalty)18 x 1 sparse Matrix of class "dgCMatrix"
s=0.3393222
(Intercept) 14.7772650
customer_age .
age_sq 1.3738993
total_bill .
log_total_bill -1.6723356
group_size .
number_drinks .
customer_sex_Male .
smoker_Yes .
day_Sat .
day_Sun 0.5933588
day_Thur .
day_Tue -1.0259372
day_Wed -1.7496789
time_Dinner -0.4028048
server_sex_Male -0.1124560
any_children_Yes .
ordered_dessert_Yes -0.4142917This provided us with some very interesting results. First, the top performing model included relatively similar features to what we found through exploratory analyses. When comparing the top performing LASSO configuration to our top performing manually constructed configuration, we saw that LASSO retained log-transformed total bill and certain days (Tuesday and Wednesday as expected low-tip days, but also retained Sunday as high-tip). This aligns quite well with our analysis. However, squared age (and not raw age) was retained, along with time, sex of the server, and dessert, all of which we had not retained in our earlier configuration. But most interestingly, the RMSE of the tuned L1 models never surpassed our minimum RMSE of 6.26, instead stagnating around 6.35, with 6.34 being the lowest. We’ll employ the same pipeline one more time while now allowing it to consider a few potentially notable interactions to see if there are further RMSE improvements or any otherwise major model configuration changes.
l1_rec_interactions <- recipe(
tip_percentage ~ customer_age + age_sq + total_bill + log_total_bill +
customer_sex + smoker + day + time + server_sex +
any_children + ordered_dessert + group_size +
number_drinks,
data = data) |>
step_dummy(all_nominal_predictors()) |>
step_interact(
terms = ~ customer_age:smoker_Yes +
customer_age:customer_sex_Male +
customer_age:log_total_bill) |>
step_normalize(all_numeric_predictors())
# lasso specification
lasso_mod <- linear_reg(
penalty = tune(),
mixture = 1) |>
set_engine("glmnet")
# workflow
lasso_wf <- workflow() |>
add_recipe(l1_rec_interactions) |>
add_model(lasso_mod)
# tuning grid
lambda_grid <- grid_regular(
penalty(range = c(-5, 1)),
levels = 50)
# tune with repeated CV
lasso_results <- tune_grid(
lasso_wf,
resamples = folds,
grid = lambda_grid,
metrics = metric_set(rmse))
# best lambda by RMSE
show_best(lasso_results, metric = "rmse")# A tibble: 5 × 7
penalty .metric .estimator mean n std_err .config
<dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
1 0.339 rmse standard 6.35 25 0.447 pre0_mod38_post0
2 0.193 rmse standard 6.36 25 0.434 pre0_mod36_post0
3 0.450 rmse standard 6.37 25 0.449 pre0_mod39_post0
4 0.256 rmse standard 6.37 25 0.441 pre0_mod37_post0
5 0.146 rmse standard 6.38 25 0.425 pre0_mod35_post0# select best model
best_lambda <- select_best(lasso_results,metric = "rmse")
# finalize workflow
final_lasso_wf <- finalize_workflow(lasso_wf, best_lambda)
# fit on full dataset
final_lasso_fit <- fit(
final_lasso_wf,
data = data)
# inspect coefficients
coef(final_lasso_fit$fit$fit$fit, s = best_lambda$penalty)21 x 1 sparse Matrix of class "dgCMatrix"
s=0.3393222
(Intercept) 14.7772650
customer_age .
age_sq 1.3738993
total_bill .
log_total_bill -1.6723356
group_size .
number_drinks .
customer_sex_Male .
smoker_Yes .
day_Sat .
day_Sun 0.5933588
day_Thur .
day_Tue -1.0259372
day_Wed -1.7496789
time_Dinner -0.4028048
server_sex_Male -0.1124560
any_children_Yes .
ordered_dessert_Yes -0.4142917
customer_age_x_smoker_Yes .
customer_age_x_customer_sex_Male .
customer_age_x_log_total_bill . We see very similar feature selection here, with no interactions being retained. All RMSE values are still higher than our manual model. Given the simplicity and interpretability of our manual model, we take this as additional evidence in favor of our earlier model configuration.
Summary
Candidate model configurations:
We considered a series of increasingly complex linear model configurations to determine how best to capture the relationship between customer age and tip percentage. We varied several analytic decisions: 1) Whether age should be modeled linearly or quadratically, 2) Whether total bill should be included in raw or log-transformed form, 3) Whether categorical variables such as day, smoker, and time should be included as covariates, and 4) Whether interactions between age and these variables improved model performance. We were pointed towards considering these specific decisions through our initial EDA.
We began with a compact baseline model containing only customer age, then compared models including quadratic age terms, total bill, log-transformed total bill, day, smoker, and combinations of these predictors. We additionally evaluated interaction models such as age × smoker, age × day, and age × log(total bill). Early exploratory analyses visually suggested that the age relationship may be nonlinear, but cross-validation ultimately suggested that a simpler linear age model generalized better. The best-performing configuration included customer age, smoking status, the age × smoker interaction, log-transformed total bill, and day of the week.
Resampling approach:
We used repeated 5-fold cross-validation, repeating five times for a total of 25 held-out performance estimates per model configuration. We chose this approach because we wanted low variance performance estimates while still retaining enough training data within each fold to fit relatively complex interaction models. A single train/test split would have produced estimates that were much higher variance, whereas repeated cross-validation reduces that variance by averaging across many held-out samples. This is particularly useful when comparing many analytic decisions that could otherwise show occasional optimization bias when they happen to fit the noise within a single dataset. The primary cost is computational complexity and runtime, though this was manageable given the relatively small dataset size.
Model selection:
We used root mean squared error (RMSE) as our primary model selection metric. RMSE measures the average magnitude of prediction error in the original units of the outcome variable (tip percentage), making it easily interpretable. More importantly, RMSE was used over a metric like sum of squared error (SSE) because SSE depends directly on sample size which may make comparing across folds less valid. RMSE therefore provides a valid, interpretable, and standardized metric for comparing model configurations.
Conclusions:
We have answered the explanatory question to a meaningful degree. We isolated the effect of age on tip percentage, and did so through effective cross-validation to select an optimal model configuration, though it is possible that we did not identify the “true” optimal model to isolate the age effect. In hopes to further explore this possibility, we also tuned a linear regression with an L1 penalty to compare feature selection and RMSE, where we obtained further evidence supporting our EDA-informed model configuration. However, it is still entirely possible that an alternative untested model configuration could better capture the underlying DGP and better isolate the age effect.
Our final model suggested that the relationship between age and tipping differs between smokers and non-smokers. Specifically, age showed little relationship with tip percentage among non-smokers (b = 0.027, p = 0.76), whereas the significant positive interaction between age and smoking status (b = 0.622, p < .001) suggested that the relationship between age and tip percentage was substantially stronger among smokers. Further, log-transformed total bill and day of the week consistently improved model performance, suggesting that contextual factors significantly influence tipping behavior as Tuesday and Wednesday were associated with notably lower tip percentages relative to other days. Regularized model tuning further indicated Sunday as a feature of interest, suggesting that it’s associated with notably higher tip percentages while controlling for certain features. Additional data, modeling, and replication would strengthen our confidence in these findings. Future analyses could also investigate more flexible nonparametric approaches, mixed-effects models, or simply seek to replicate across different restaurants or geographic regions to reach a more robust, grounded conclusion.
Reflection:
Using cross-validation to compare candidate model configurations allowed us to evaluate analytic decisions empirically rather than relying entirely on subjective judgment, arbitrary preregistration choices, or wholly automated feature selection. This reduced the likelihood of selecting overly complex, uninterpretable models that fit noise in the dataset and it also allows us to investigate which transformations, covariates, and interactions meaningfully improve generalization. However, the process did highlight how exploratory visual analysis can sometimes be misleading. Initially, the relationship between age and tip percentage appeared potentially nonlinear, which led us down the wrong path before easing up on this assumption. This demonstrated the importance of evaluating candidate explanations systematically rather than relying solely on visual inspection. Further, repeatedly comparing many candidate models still introduces researcher degrees of freedom, particularly when model configurations are influenced by exploratory analysis. Cross-validation also evaluates predictive generalization rather than causal validity, so strong predictive performance does not necessarily imply that the identified relationships are causal or otherwise “true”. Finally, selecting models partially based on predictive performance may favor variables that improve prediction without necessarily being theoretically meaningful.