Predicting price of diamonds

Business Analytics (23/24)

In this document, we will perform a numeric prediction job using the diamonds dataset from ggplot2. As ggplot2 is included in tidymodels, we will need only the tidymodels package.

Code 
library(tidymodels)
Warning: package 'tidymodels' was built under R version 4.3.3
Warning: package 'dials' was built under R version 4.3.3
Warning: package 'ggplot2' was built under R version 4.3.3
Warning: package 'infer' was built under R version 4.3.3
Warning: package 'modeldata' was built under R version 4.3.3
Warning: package 'parsnip' was built under R version 4.3.3
Warning: package 'recipes' was built under R version 4.3.3
Warning: package 'rsample' was built under R version 4.3.3
Warning: package 'tune' was built under R version 4.3.3
Warning: package 'workflows' was built under R version 4.3.3
Warning: package 'workflowsets' was built under R version 4.3.3
Warning: package 'yardstick' was built under R version 4.3.3

1 The diamonds Dataset

The dataset contains the prices and other attributes of almost 54,000 diamonds.

Code 
data("diamonds")
diamonds
# A tibble: 53,940 × 10
   carat cut       color clarity depth table price     x     y     z
   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
 1  0.23 Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
 2  0.21 Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
 3  0.23 Good      E     VS1      56.9    65   327  4.05  4.07  2.31
 4  0.29 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
 5  0.31 Good      J     SI2      63.3    58   335  4.34  4.35  2.75
 6  0.24 Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48
 7  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47
 8  0.26 Very Good H     SI1      61.9    55   337  4.07  4.11  2.53
 9  0.22 Fair      E     VS2      65.1    61   337  3.87  3.78  2.49
10  0.23 Very Good H     VS1      59.4    61   338  4     4.05  2.39
# ℹ 53,930 more rows

1.1 The target variable

Our job will be predicting the price variable, thus it is a regression or numerical prediction job. It is always a good idea to examine the target variable.

Code 
diamonds |>
  ggplot(aes(price)) +
  geom_histogram(bins = 20) +
  theme_minimal()

We can see that the variable is highly right-skewed, meaning that there are many samples with a price higher than the expected if the variable followed a normal law.

In those cases, it is a good idea to use the logarithm of the variable as a predictor. The transformations of the target variable must be done outside the prediction workflow, so we define log_price as:

Code 
diamonds <- diamonds |>
  mutate(log_price = log(price))

The distribution of log_price is:

Code 
diamonds |>
  ggplot(aes(log_price)) +
  geom_histogram(bins = 20) +
  theme_minimal()

We observe that log_price has a more adequate distribution for predictive modelling.

1.2 The Features

The features of the dataset are some diamond properties. Features from cut to clarity are set as ordered factors, meaning that factor levels are ordinal variables. This means that we can turn them into ordinal variables, rather than a set of dummies.

The other features are numeric, and have to do with diamond size. carat is a measure of diamond weight, and x,y, z, depth and table are measures of diamond size.

As larger diamonds will also be larger, it is not surprising to find high values of correlation between those variables. Let’s use the corrr package to see them.

Code 
library(corrr)
Warning: package 'corrr' was built under R version 4.3.3
Code 
diamonds |>
  select(where(is.numeric)) |>
  correlate() |>
  rearrange() |>
  shave() |>
  rplot()

We observe that carat, and x to z are highly correlated, so some of them will be filtered by a step_corr() recipe.

Let’s see how the geometrical x variable relates with log_price.

Code 
diamonds |>
  ggplot(aes(x, log_price)) +
  geom_point() +
  geom_smooth() +
  theme_minimal()
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

There are some diamonds of geometry zero.

Code 
diamonds |>
  filter(x == 0 | y == 0 | z == 0)
# A tibble: 20 × 11
   carat cut       color clarity depth table price     x     y     z log_price
   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>     <dbl>
 1  1    Premium   G     SI2      59.1    59  3142  6.55  6.48     0      8.05
 2  1.01 Premium   H     I1       58.1    59  3167  6.66  6.6      0      8.06
 3  1.1  Premium   G     SI2      63      59  3696  6.5   6.47     0      8.22
 4  1.01 Premium   F     SI2      59.2    58  3837  6.5   6.47     0      8.25
 5  1.5  Good      G     I1       64      61  4731  7.15  7.04     0      8.46
 6  1.07 Ideal     F     SI2      61.6    56  4954  0     6.62     0      8.51
 7  1    Very Good H     VS2      63.3    53  5139  0     0        0      8.54
 8  1.15 Ideal     G     VS2      59.2    56  5564  6.88  6.83     0      8.62
 9  1.14 Fair      G     VS1      57.5    67  6381  0     0        0      8.76
10  2.18 Premium   H     SI2      59.4    61 12631  8.49  8.45     0      9.44
11  1.56 Ideal     G     VS2      62.2    54 12800  0     0        0      9.46
12  2.25 Premium   I     SI1      61.3    58 15397  8.52  8.42     0      9.64
13  1.2  Premium   D     VVS1     62.1    59 15686  0     0        0      9.66
14  2.2  Premium   H     SI1      61.2    59 17265  8.42  8.37     0      9.76
15  2.25 Premium   H     SI2      62.8    59 18034  0     0        0      9.80
16  2.02 Premium   H     VS2      62.7    53 18207  8.02  7.95     0      9.81
17  2.8  Good      G     SI2      63.8    58 18788  8.9   8.85     0      9.84
18  0.71 Good      F     SI2      64.1    60  2130  0     0        0      7.66
19  0.71 Good      F     SI2      64.1    60  2130  0     0        0      7.66
20  1.12 Premium   G     I1       60.4    59  2383  6.71  6.67     0      7.78

Let’s see how carat works:

Code 
diamonds |>
  ggplot(aes(carat, log_price)) +
  geom_point() +
  geom_smooth() +
  theme_minimal()
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

We observe a strong, nonlinear relationship between carat and price. carat can be a better feature than x:z variables, as it has no anomalous zero values, so we will skip them and keep carat as feature.

2 Predicting Diamond Prices

Let’s define the elements of a workflow to predict diamond prices.

2.1 Initial Split

Performing an adequate split of the dataset into train and test sets. Keeping 90% of data for the training set.

Code 
set.seed(44)
d_split <- initial_split(diamonds, prop = 0.9)

2.2 Preprocessing

Some useful steps for preprocessing:

Performing preprocessing steps in a recipe.

Code 
d_rec <- recipe(log_price ~ ., training(d_split)) |>
  update_role(price, new_role = "original predictor") |>
  step_ordinalscore(all_nominal_predictors()) |>
   step_rm(x:z) |>
  step_sqrt(carat) |>
  step_nzv(all_predictors())

We can see the preprocessed dataset doing:

Code 
d_rec |>
  prep() |>
  bake(new_data = NULL)
# A tibble: 48,546 × 8
   carat   cut color clarity depth table price log_price
   <dbl> <int> <int>   <int> <dbl> <dbl> <int>     <dbl>
 1 1.24      4     7       4  62.9    59  7554      8.93
 2 0.616     5     2       5  62.1    56  1000      6.91
 3 1.01      4     3       2  61.4    57  4586      8.43
 4 0.548     2     2       5  63.2    57   568      6.34
 5 0.566     3     3       6  63.8    55   701      6.55
 6 0.469     1     2       4  65.1    61   337      5.82
 7 0.566     4     1       3  60.9    60   589      6.38
 8 1.00      5     4       3  62.9    56  5107      8.54
 9 0.548     3     1       4  62.5    57   670      6.51
10 1.01      5     4       8  62.5    57  8162      9.01
# ℹ 48,536 more rows

And we can see the roles of variables:

Code 
d_rec |>
  prep() |> 
  summary()
# A tibble: 8 × 4
  variable  type      role               source  
  <chr>     <list>    <chr>              <chr>   
1 carat     <chr [2]> predictor          original
2 cut       <chr [2]> predictor          original
3 color     <chr [2]> predictor          original
4 clarity   <chr [2]> predictor          original
5 depth     <chr [2]> predictor          original
6 table     <chr [2]> predictor          original
7 price     <chr [2]> original predictor original
8 log_price <chr [2]> outcome            original

2.3 Models

Define=ing two predictive models:

  • A ensemble-based model, like boosted trees.
  • A regression-based model, like regularized regression.

I have chosen a regularized regression and a boosted tree model.

Code 
rr <- linear_reg(mode = "regression", engine = "glmnet", penalty = 0, mixture = 1)

bt <- boost_tree(mode = "regression") |>
  set_engine("xgboost")

And the workflows:

Code 
d_rr_wf <- workflow() |>
  add_recipe(d_rec) |>
  add_model(rr)

d_bt_wf <- workflow() |>
  add_recipe(d_rec) |>
  add_model(bt)

Let’s store both workflows into a list.

Code 
d_wf <- list(d_rr_wf, d_bt_wf)
names(d_wf) <- c("reg_regression", "boosted_trees")

2.4 Cross Validation

Defining a set of ten folders for the training set, and testing the two models with cross validation using as metrics mean absolute error mae, root of mean squared errors rmse and r squared rsq.

The folds:

Code 
set.seed(11)
folds <- vfold_cv(training(d_split), v = 10)

The metrics:

Code 
reg_metrics <- metric_set(mae, rmse, rsq)

Let’s test the two models at the same time with lapply().

Code 
d_cv <- lapply(d_wf, \(m) m |> fit_resamples(folds, metrics = reg_metrics) |> collect_metrics())
Warning: package 'glmnet' was built under R version 4.3.3
Warning: package 'xgboost' was built under R version 4.3.3

2.5 Model Decision

Let’s see the results:

Code 
d_cv
$reg_regression
# A tibble: 3 × 6
  .metric .estimator  mean     n  std_err .config             
  <chr>   <chr>      <dbl> <int>    <dbl> <chr>               
1 mae     standard   0.166    10 0.000860 Preprocessor1_Model1
2 rmse    standard   0.215    10 0.00157  Preprocessor1_Model1
3 rsq     standard   0.955    10 0.000533 Preprocessor1_Model1

$boosted_trees
# A tibble: 3 × 6
  .metric .estimator   mean     n  std_err .config             
  <chr>   <chr>       <dbl> <int>    <dbl> <chr>               
1 mae     standard   0.0858    10 0.000314 Preprocessor1_Model1
2 rmse    standard   0.110     10 0.000568 Preprocessor1_Model1
3 rsq     standard   0.990     10 0.000114 Preprocessor1_Model1

We observe that the boosted trees model is the one that performs best.

3 Final Model

3.1 Training the Model

Training a model for diamond pricing on the whole training set.

Code 
model <- d_bt_wf |>
  fit(training(d_split))

3.2 Performance of log_price

Evaluating the performance of the model in the test set.

Code 
model |>
  predict(testing(d_split)) |>
  bind_cols(testing(d_split)) |>
  reg_metrics(truth = log_price, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 mae     standard      0.0878
2 rmse    standard      0.113 
3 rsq     standard      0.989

We can plot target versus prediction:

Code 
model |>
  predict(testing(d_split)) |>
  bind_cols(testing(d_split)) |>
  ggplot(aes(log_price, .pred)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()

3.3 Performance with price

Evaluating the performance of the model on the test set using the original price variable and the exponential transformation of the prediction.

We need to obtain the real price using the reverse function of log() that is exp().

Code 
model |>
  predict(testing(d_split)) |>
  bind_cols(testing(d_split)) |>
  mutate(price_pred = exp(.pred)) |>
  reg_metrics(truth = price, estimate = price_pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 mae     standard     340.   
2 rmse    standard     632.   
3 rsq     standard       0.979

The target versus prediction plot shows high variability in prediction for high prices.

Code 
model |>
  predict(testing(d_split)) |>
  bind_cols(testing(d_split)) |>
  mutate(price_pred = exp(.pred)) |>
  ggplot(aes(price, price_pred)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()