Chapter 3: Feature & Target Engineering

Note: Some results may differ from the hard copy book due to the changing of sampling procedures introduced in R 3.6.0. See http://bit.ly/35D1SW7 for more details. Access and run the source code for this notebook here.

Hidden chapter requirements used in the book to set the plotting theme and load packages used in hidden code chunks.

# additional package required in hidden code chunks
library(purrr)
library(tidyr)
library(readr)
library(kableExtra)

# Set the graphical theme
ggplot2::theme_set(ggplot2::theme_light())

Prerequisites

This chapter leverages the following packages:

# Helper packages
library(dplyr)    # for data manipulation
library(ggplot2)  # for awesome graphics
library(visdat)   # for additional visualizations
# Feature engineering packages
library(caret)    # for various ML tasks
library(recipes)  # for feature engineering tasks

We’ll also continue working with the ames_train data set:

ames <- AmesHousing::make_ames()

# Load and split the Ames housing data using stratified sampling
set.seed(123)  # for reproducibility
split  <- rsample::initial_split(ames, prop = 0.7, strata = "Sale_Price")
ames_train  <- rsample::training(split)
ames_test   <- rsample::testing(split)

Target engineering

Figure 3.1:

models <- c("Non-log transformed model residuals", 
            "Log transformed model residuals")
list(
  m1 = lm(Sale_Price ~ Year_Built, data = ames_train),
  m2 = lm(log(Sale_Price) ~ Year_Built, data = ames_train)
) %>%
  map2_dfr(models, ~ broom::augment(.x) %>% mutate(model = .y)) %>%
  ggplot(aes(.resid)) +
    geom_histogram(bins = 75) +
    facet_wrap(~ model, scales = "free_x") +
    ylab(NULL) +
    xlab("Residuals")

transformed_response <- log(ames_train$Sale_Price)

# log transformation
ames_recipe <- recipe(Sale_Price ~ ., data = ames_train) %>%
  step_log(all_outcomes())
ames_recipe

Data Recipe

Inputs:

Operations:

Log transformation on all_outcomes

log(-0.5)

NaNs produced

[1] NaN

log1p(-0.5)

[1] -0.6931472

\[ \begin{equation} y(\lambda) = \begin{cases} \frac{Y^\lambda-1}{\lambda}, & \text{if}\ \lambda \neq 0 \\ \log\left(Y\right), & \text{if}\ \lambda = 0. \end{cases} \end{equation} \]

Figure 3.2

# Log transformation
train_log_y <- log(ames_train$Sale_Price)
test_log_y  <- log(ames_train$Sale_Price)

# Box Cox transformation
lambda  <- forecast::BoxCox.lambda(ames_train$Sale_Price)
train_bc_y <- forecast::BoxCox(ames_train$Sale_Price, lambda)
test_bc_y  <- forecast::BoxCox(ames_test$Sale_Price, lambda)

# Plot differences
levs <- c("Normal", "Log_Transform", "BoxCox_Transform")
data.frame(
  Normal = ames_train$Sale_Price,
  Log_Transform = train_log_y,
  BoxCox_Transform = train_bc_y
) %>%
  gather(Transform, Value) %>%
  mutate(Transform = factor(Transform, levels = levs)) %>% 
  ggplot(aes(Value, fill = Transform)) +
    geom_histogram(show.legend = FALSE, bins = 40) +
    facet_wrap(~ Transform, scales = "free_x")

# Log transform a value
y <- log(10)

# Undo log-transformation
exp(y)

[1] 10

# Box Cox transform a value
y <- forecast::BoxCox(10, lambda)

# Inverse Box Cox function
inv_box_cox <- function(x, lambda) {
  # for Box-Cox, lambda = 0 --> log transform
  if (lambda == 0) exp(x) else (lambda*x + 1)^(1/lambda) 
}

# Undo Box Cox-transformation
inv_box_cox(y, lambda)

[1] 10
attr(,"lambda")
[1] -0.05425403

Dealing with missingness

sum(is.na(AmesHousing::ames_raw))

[1] 13997

Figure 3.3:

AmesHousing::ames_raw %>%
  is.na() %>%
  reshape2::melt() %>%
  ggplot(aes(Var2, Var1, fill=value)) + 
    geom_raster() + 
    coord_flip() +
    scale_y_continuous(NULL, expand = c(0, 0)) +
    scale_fill_grey(name = "", 
                    labels = c("Present", 
                               "Missing")) +
    xlab("Observation") +
    theme(axis.text.y  = element_text(size = 4))

AmesHousing::ames_raw %>% 
  filter(is.na(`Garage Type`)) %>% 
  select(`Garage Type`, `Garage Cars`, `Garage Area`)

Figure 3.4:

vis_miss(AmesHousing::ames_raw, cluster = TRUE)

Imputation

ames_recipe %>%
  step_medianimpute(Gr_Liv_Area)

Data Recipe

Inputs:

Operations:

Log transformation on all_outcomes
Median Imputation for Gr_Liv_Area

ames_recipe %>%
  step_knnimpute(all_predictors(), neighbors = 6)

Data Recipe

Inputs:

Operations:

Log transformation on all_outcomes
K-nearest neighbor imputation for all_predictors

ames_recipe %>%
  step_bagimpute(all_predictors())

Data Recipe

Inputs:

Operations:

Log transformation on all_outcomes
Bagged tree imputation for all_predictors

Figure 3.5:

impute_ames <- ames_train
set.seed(123)
index <- sample(seq_along(impute_ames$Gr_Liv_Area), 50)
actuals <- ames_train[index, ]
impute_ames$Gr_Liv_Area[index] <- NA

p1 <- ggplot() +
  geom_point(data = impute_ames, aes(Gr_Liv_Area, Sale_Price), alpha = .2) +
  geom_point(data = actuals, aes(Gr_Liv_Area, Sale_Price), color = "red") +
  scale_x_log10(limits = c(300, 5000)) +
  scale_y_log10(limits = c(10000, 500000)) +
  ggtitle("Actual values")

# Mean imputation
mean_juiced <- recipe(Sale_Price ~ ., data = impute_ames) %>%
  step_meanimpute(Gr_Liv_Area) %>%
  prep(training = impute_ames, retain = TRUE) %>%
  juice()
mean_impute <- mean_juiced[index, ]
  
p2 <- ggplot() +
  geom_point(data = actuals, aes(Gr_Liv_Area, Sale_Price), color = "red") +
  geom_point(data = mean_impute, aes(Gr_Liv_Area, Sale_Price), color = "blue") +
  scale_x_log10(limits = c(300, 5000)) +
  scale_y_log10(limits = c(10000, 500000)) +
  ggtitle("Mean Imputation")

# KNN imputation
knn_juiced <- recipe(Sale_Price ~ ., data = impute_ames) %>%
  step_knnimpute(Gr_Liv_Area) %>%
  prep(training = impute_ames, retain = TRUE) %>%
  juice()
knn_impute <- knn_juiced[index, ]
  
p3 <- ggplot() +
  geom_point(data = actuals, aes(Gr_Liv_Area, Sale_Price), color = "red") +
  geom_point(data = knn_impute, aes(Gr_Liv_Area, Sale_Price), color = "blue") +
  scale_x_log10(limits = c(300, 5000)) +
  scale_y_log10(limits = c(10000, 500000)) +
  ggtitle("KNN Imputation")

# Bagged imputation
bagged_juiced <- recipe(Sale_Price ~ ., data = impute_ames) %>%
  step_bagimpute(Gr_Liv_Area) %>%
  prep(training = impute_ames, retain = TRUE) %>%
  juice()
bagged_impute <- bagged_juiced[index, ]
  
p4 <- ggplot() +
  geom_point(data = actuals, aes(Gr_Liv_Area, Sale_Price), color = "red") +
  geom_point(data = bagged_impute, aes(Gr_Liv_Area, Sale_Price), color = "blue") +
  scale_x_log10(limits = c(300, 5000)) +
  scale_y_log10(limits = c(10000, 500000)) +
  ggtitle("Bagged Trees Imputation")

gridExtra::grid.arrange(p1, p2, p3, p4, nrow = 2)

Feature filtering

Figure 3.6 (generated from benchmark data performed in the past):

model_results <- read_csv("data/feature-selection-impacts-results.csv") %>%
  mutate(type = case_when(
    model %in% c("lm", "pcr", "pls", "glmnet", "lasso") ~ "Linear models",
    model %in% c("earth", "svmLinear", "nn") ~ "Non-linear Models",
    TRUE ~ "Tree-based Models"
  )) %>%
  mutate(model = case_when(
    model == "lm" ~ "Linear regression",
    model == "earth" ~ "Multivariate adaptive regression splines",
    model == "gbm" ~ "Gradient boosting machines",
    model == "glmnet" ~ "Elastic net",
    model == "lasso" ~ "Lasso",
    model == "nn" ~ "Neural net",
    model == "pcr" ~ "Principal component regression",
    model == "pls" ~ "Partial least squares",
    model == "ranger" ~ "Random forest",
    TRUE ~ "Support vector machine"
  ))

Parsed with column specification:
cols(
  model = col_character(),
  NIP = col_double(),
  RMSE = col_double(),
  time = col_double()
)

ggplot(model_results, aes(NIP, RMSE, color = model, lty = model)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ type, nrow = 1) +
  xlab("Number of additional non-informative predictors")

Figure 3.7 (generated from benchmark data performed in the past):

model_results %>%
  group_by(model) %>%
  mutate(
    time_impact = time / first(time),
    time_impact = time_impact - 1
  ) %>%
  ggplot(aes(NIP, time_impact, color = model, lty = model)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ type, nrow = 1) +
    scale_y_continuous("Percent increase in training time", 
                       labels = scales::percent) +
    xlab("Number of additional non-informative predictors")

caret::nearZeroVar(ames_train, saveMetrics = TRUE) %>% 
  tibble::rownames_to_column() %>% 
  filter(nzv)

Numeric feature engineering

# Normalize all numeric columns
recipe(Sale_Price ~ ., data = ames_train) %>%
  step_YeoJohnson(all_numeric())                 

Data Recipe

Inputs:

Operations:

Yeo-Johnson transformation on all_numeric

Figure 3.8:

set.seed(123)
x1 <- tibble(
  variable = "x1",
  `Real value` = runif(25, min = -30, max = 5),
  `Standardized value` = scale(`Real value`) %>% as.numeric()
)

set.seed(456)
x2 <- tibble(
  variable = "x2",
  `Real value` = rlnorm(25, log(25)),
  `Standardized value` = scale(`Real value`) %>% as.numeric()
)

set.seed(789)
x3 <- tibble(
  variable = "x3",
  `Real value` = rnorm(25, 150, 15),
  `Standardized value` = scale(`Real value`) %>% as.numeric()
)

x1 %>%
  bind_rows(x2) %>%
  bind_rows(x3) %>%
  gather(key, value, -variable) %>%
  mutate(variable = factor(variable, levels = c("x3", "x2", "x1"))) %>%
  ggplot(aes(value, variable)) +
    geom_point(alpha = .6) +
    facet_wrap(~ key, scales = "free_x") +
    ylab("Feature") +
    xlab("Value")

ames_recipe %>%
  step_center(all_numeric(), -all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes())

Data Recipe

Inputs:

Operations:

Log transformation on all_outcomes
Centering for all_numeric, -, all_outcomes()
Scaling for all_numeric, -, all_outcomes()

Categorical feature engineering

count(ames_train, Neighborhood) %>% arrange(n)

count(ames_train, Screen_Porch) %>% arrange(n)

# Lump levels for two features
lumping <- recipe(Sale_Price ~ ., data = ames_train) %>%
  step_other(Neighborhood, threshold = 0.01, 
             other = "other") %>%
  step_other(Screen_Porch, threshold = 0.1, 
             other = ">0")

# Apply this blue print --> you will learn about this at 
# the end of the chapter
apply_2_training <- prep(lumping, training = ames_train) %>%
  bake(ames_train)

# New distribution of Neighborhood
count(apply_2_training, Neighborhood) %>% arrange(n)

# New distribution of Screen_Porch
count(apply_2_training, Screen_Porch) %>% arrange(n)

knitr::include_graphics("images/ohe-vs-dummy.png")

# Lump levels for two features
recipe(Sale_Price ~ ., data = ames_train) %>%
  step_dummy(all_nominal(), one_hot = TRUE)

Data Recipe

Inputs:

Operations:

Dummy variables from all_nominal

# Original categories
count(ames_train, MS_SubClass)

# Label encoded
recipe(Sale_Price ~ ., data = ames_train) %>%
  step_integer(MS_SubClass) %>%
  prep(ames_train) %>%
  bake(ames_train) %>%
  count(MS_SubClass)

ames_train %>% select(contains("Qual"))

# Original categories
count(ames_train, Overall_Qual)

# Label encoded
recipe(Sale_Price ~ ., data = ames_train) %>%
  step_integer(Overall_Qual) %>%
  prep(ames_train) %>%
  bake(ames_train) %>%
  count(Overall_Qual)

ames_train %>%
  group_by(Neighborhood) %>%
  summarize(`Avg Sale_Price` = mean(Sale_Price, na.rm = TRUE)) %>%
  head(10) %>%
  kable(caption = "Example of target encoding the Neighborhood feature of the Ames housing data set.") %>%
  kable_styling(bootstrap_options = "striped", full_width = TRUE)

Example of target encoding the Neighborhood feature of the Ames housing data set.

Neighborhood	Avg Sale_Price
North_Ames	144562.7
College_Creek	199831.7
Old_Town	122736.7
Edwards	130652.2
Somerset	227379.6
Northridge_Heights	323289.5
Gilbert	192162.9
Sawyer	136460.6
Northwest_Ames	187328.2
Sawyer_West	188644.6

ames_train %>%
  count(Neighborhood) %>%
  mutate(Proportion = n / sum(n)) %>%
  select(-n) %>%
  head(10) %>%
  kable(caption = 'Example of categorical proportion encoding the Neighborhood feature of the Ames housing data set.') %>%
  kable_styling(bootstrap_options = "striped", full_width = TRUE)

Example of categorical proportion encoding the Neighborhood feature of the Ames housing data set.

Neighborhood	Proportion
North_Ames	0.1451534
College_Creek	0.0910862
Old_Town	0.0832927
Edwards	0.0711154
Somerset	0.0623478
Northridge_Heights	0.0560156
Gilbert	0.0565027
Sawyer	0.0496834
Northwest_Ames	0.0467608
Sawyer_West	0.0414028

Dimension reduction

recipe(Sale_Price ~ ., data = ames_train) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  step_pca(all_numeric(), threshold = .95)

Data Recipe

Inputs:

Operations:

Centering for all_numeric
Scaling for all_numeric
No PCA components were extracted.

Proper implementation

Figure 3.10:

knitr::include_graphics("images/minimize-leakage.png")

blueprint <- recipe(Sale_Price ~ ., data = ames_train) %>%
  step_nzv(all_nominal())  %>%
  step_integer(matches("Qual|Cond|QC|Qu")) %>%
  step_center(all_numeric(), -all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes()) %>%
  step_pca(all_numeric(), -all_outcomes())
  
blueprint

Data Recipe

Inputs:

Operations:

Sparse, unbalanced variable filter on all_nominal
Integer encoding for matches, Qual|Cond|QC|Qu
Centering for all_numeric, -, all_outcomes()
Scaling for all_numeric, -, all_outcomes()
No PCA components were extracted.

prepare <- prep(blueprint, training = ames_train)
prepare

Data Recipe

Inputs:

Training data contained 2053 data points and no missing data.

Operations:

Sparse, unbalanced variable filter removed Street, Alley, Land_Contour, Utilities, ... [trained]
Integer encoding for Condition_1, Overall_Qual, Overall_Cond, Exter_Qual, Exter_Cond, ... [trained]
Centering for Lot_Frontage, Lot_Area, Condition_1, Overall_Qual, ... [trained]
Scaling for Lot_Frontage, Lot_Area, Condition_1, Overall_Qual, ... [trained]
PCA extraction with Lot_Frontage, Lot_Area, Condition_1, Overall_Qual, ... [trained]

baked_train <- bake(prepare, new_data = ames_train)
baked_test <- bake(prepare, new_data = ames_test)
baked_train

blueprint <- recipe(Sale_Price ~ ., data = ames_train) %>%
  step_nzv(all_nominal()) %>%
  step_integer(matches("Qual|Cond|QC|Qu")) %>%
  step_center(all_numeric(), -all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes()) %>%
  step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE)

# Specify resampling plan
cv <- trainControl(
  method = "repeatedcv", 
  number = 10, 
  repeats = 5
)

# Construct grid of hyperparameter values
hyper_grid <- expand.grid(k = seq(2, 25, by = 1))

# Tune a knn model using grid search
knn_fit2 <- train(
  blueprint, 
  data = ames_train, 
  method = "knn", 
  trControl = cv, 
  tuneGrid = hyper_grid,
  metric = "RMSE"
)

# print model results
knn_fit2

k-Nearest Neighbors 

2053 samples
  80 predictor

Recipe steps: nzv, integer, center, scale, dummy 
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 1848, 1849, 1848, 1847, 1848, 1848, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   2  35939.89  0.8048642  22437.87
   3  35010.19  0.8158644  21727.75
   4  34397.34  0.8237016  21241.75
   5  33930.03  0.8316443  20941.83
   6  33546.26  0.8383965  20777.11
   7  33336.67  0.8433273  20688.80
   8  33132.85  0.8467156  20578.56
   9  33054.25  0.8486616  20536.96
  10  32948.08  0.8505378  20529.58
  11  32902.66  0.8522363  20537.99
  12  32812.52  0.8543514  20571.13
  13  32835.43  0.8554652  20625.39
  14  32859.76  0.8564118  20685.16
  15  32890.51  0.8571297  20728.81
  16  32939.40  0.8577373  20774.40
  17  32982.90  0.8578918  20817.29
  18  33040.08  0.8579879  20910.09
  19  33113.72  0.8579053  20977.65
  20  33210.68  0.8574767  21057.82
  21  33274.76  0.8570757  21110.48
  22  33296.28  0.8574598  21162.57
  23  33397.05  0.8571769  21252.19
  24  33421.53  0.8575578  21282.18
  25  33461.60  0.8574388  21320.31

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 12.

# plot cross validation results
ggplot(knn_fit2)