The goal is to use machine learning techniques to obtain the predicted diamond price with the minimum errors from the actual price, which involved a complete process of exploratory data analysis, feature engineering, feature selection and predictive modeling analysis.

We assume that the actual price in the historical diamond dataset reflects the current value without other factors such as inflation.

set.seed(503)
library(tidyverse)
diamonds_test <- diamonds %>% mutate(diamond_id = row_number()) %>% 
    group_by(cut, color, clarity) %>% sample_frac(0.2) %>% ungroup()

diamonds_train <- anti_join(diamonds %>% mutate(diamond_id = row_number()), 
    diamonds_test, by = "diamond_id")
set.seed(503)
diamonds_test <- diamonds %>% mutate(diamond_id = row_number()) %>% 
    group_by(cut, color, clarity) %>% sample_frac(0.2) %>% ungroup()
diamonds_train <- anti_join(diamonds %>% mutate(diamond_id = row_number()), 
    diamonds_test, by = "diamond_id")
diamonds_all <- bind_rows(diamonds_train,diamonds_test)

Below is a sample from dataset.

  head(diamonds_train)%>%
  kbl(caption = "Sample from Dataset") %>%
  kable_minimal(full_width = F)
Sample from Dataset
carat cut color clarity depth table price x y z diamond_id
0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31 2
0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31 3
0.29 Premium I VS2 62.4 58 334 4.20 4.23 2.63 4
0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75 5
0.24 Very Good I VVS1 62.3 57 336 3.95 3.98 2.47 7
0.26 Very Good H SI1 61.9 55 337 4.07 4.11 2.53 8

Below is statistical information about our data.

summary(diamonds)
##      carat               cut        color        clarity          depth      
##  Min.   :0.2000   Fair     : 1610   D: 6775   SI1    :13065   Min.   :43.00  
##  1st Qu.:0.4000   Good     : 4906   E: 9797   VS2    :12258   1st Qu.:61.00  
##  Median :0.7000   Very Good:12082   F: 9542   SI2    : 9194   Median :61.80  
##  Mean   :0.7979   Premium  :13791   G:11292   VS1    : 8171   Mean   :61.75  
##  3rd Qu.:1.0400   Ideal    :21551   H: 8304   VVS2   : 5066   3rd Qu.:62.50  
##  Max.   :5.0100                     I: 5422   VVS1   : 3655   Max.   :79.00  
##                                     J: 2808   (Other): 2531                  
##      table           price             x                y         
##  Min.   :43.00   Min.   :  326   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:56.00   1st Qu.:  950   1st Qu.: 4.710   1st Qu.: 4.720  
##  Median :57.00   Median : 2401   Median : 5.700   Median : 5.710  
##  Mean   :57.46   Mean   : 3933   Mean   : 5.731   Mean   : 5.735  
##  3rd Qu.:59.00   3rd Qu.: 5324   3rd Qu.: 6.540   3rd Qu.: 6.540  
##  Max.   :95.00   Max.   :18823   Max.   :10.740   Max.   :58.900  
##                                                                   
##        z         
##  Min.   : 0.000  
##  1st Qu.: 2.910  
##  Median : 3.530  
##  Mean   : 3.539  
##  3rd Qu.: 4.040  
##  Max.   :31.800  
## 
set.seed(20022012)
diamond_samp <- diamonds_train[sample(1:length(diamonds_train$price), 10000), ]
ggpairs(diamond_samp, outlier.shape = I('.')) # params = c(shape = I('.'),

1 Exploring the data

Correlation of variable with Price If we examine the correlation of variables with price. We can see that there is strong positive correlation with x,y,z and carat. Correlation with table and depth variable is low.

with(diamonds_train,
     data.frame(cor_x_price = cor(x, price),
                cor_y_price = cor(y, price),
                cor_z_price = cor(z, price),
                cor_depth_price = cor(depth, price),
                cor_table_price2 = cor(table, price),
                cor_carat_price3 = cor(carat, price)
     )
)
##   cor_x_price cor_y_price cor_z_price cor_depth_price cor_table_price2
## 1   0.8852482   0.8613732   0.8580964     -0.01327053         0.128216
##   cor_carat_price3
## 1        0.9217041

Since we have a few factor variables, I will use one hot encoding to convert them to numeric type.

#*** one hot Encoding

ohe_features<-c("cut","color","clarity")
dummies<-dummyVars(~cut + color + clarity ,data=diamonds_train)

data02_ohe<-as.data.frame(predict(dummies,newdata=diamonds_train))
data03_combined<-cbind(diamonds_train,data02_ohe)

data04<-data03_combined[,!names(data03_combined)%in%ohe_features]

rm(data03_combined)
rm(data02_ohe)
## Correlation plot
y.label <-as.numeric(data04$price)
corrplot(cor(cbind(data04,Price=y.label)),type="upper")

As the carat of diamond increases price also increase as stated in correlation also. Price range is very big even at the same carat. This shows us that there are other variables that effect the price also.

ggplot(diamonds_train) + 
  geom_point(aes(x = carat, y = price), alpha = .05, color = "blue") + 
  labs(title = "Diamond Price vs. Length", x = "Length of diamond (mm)", y = "Price of diamond (USD)")

Let’s add the clarity information to plot.

ggplot(aes(x = carat, y = price), data = diamonds_train) + 
  geom_point(alpha = 0.5, size = 1, position = 'jitter',aes(color=clarity)) +
  scale_color_brewer(type = 'div',
                     guide = guide_legend(title = 'Clarity', reverse = T,
                                          override.aes = list(alpha = 1, size = 2))) +                         
  ggtitle('Price by Carat and Clarity')

qplot(price, data=diamonds_train ,geom="density", fill=cut, alpha=I(.5), 
      main="Distribution of Carat", xlab="Different kinds of cut", 
      ylab="Density") + theme_minimal()

qplot(price, data=diamonds_train, geom="density", fill=color, alpha=I(.5), 
      main="Distribution of Carat", xlab="Different Colors", 
      ylab="Density") + theme_minimal()

qplot(price, data=diamonds_train, geom="density", fill=clarity, alpha=I(.5), 
      main="Distribution of Carat", xlab="Different clarity parameters", 
      ylab="Density") + theme_minimal()

2 CART modelling

First we prepare the dataset we use to train the model then test the model.

diamond_model <- rpart(price ~ x+y+z+carat+cut+color+clarity, data=diamonds_train)

summary(diamond_model)
## Call:
## rpart(formula = price ~ x + y + z + carat + cut + color + clarity, 
##     data = diamonds_train)
##   n= 43143 
## 
##           CP nsplit rel error    xerror        xstd
## 1 0.60913716      0 1.0000000 1.0000574 0.009817577
## 2 0.18517972      1 0.3908628 0.3908901 0.004372515
## 3 0.03361554      2 0.2056831 0.2057352 0.002295471
## 4 0.02661853      3 0.1720676 0.1722087 0.002292179
## 5 0.02583190      4 0.1454491 0.1413036 0.001913812
## 6 0.01000000      5 0.1196172 0.1202350 0.001729446
## 
## Variable importance
##   carat       y       x       z clarity   color 
##      25      24      24      23       3       1 
## 
## Node number 1: 43143 observations,    complexity param=0.6091372
##   mean=3938.805, MSE=1.597919e+07 
##   left son=2 (27894 obs) right son=3 (15249 obs)
##   Primary splits:
##       carat < 0.995 to the left,  improve=0.60913720, (0 missing)
##       y     < 6.345 to the left,  improve=0.60805950, (0 missing)
##       x     < 6.315 to the left,  improve=0.60428840, (0 missing)
##       z     < 3.915 to the left,  improve=0.59942370, (0 missing)
##       color splits as  LLLLRRR,   improve=0.02261306, (0 missing)
##   Surrogate splits:
##       x       < 6.265 to the left,  agree=0.983, adj=0.953, (0 split)
##       y       < 6.285 to the left,  agree=0.981, adj=0.946, (0 split)
##       z       < 3.895 to the left,  agree=0.978, adj=0.937, (0 split)
##       clarity splits as  RRLLLLLL,  agree=0.679, adj=0.091, (0 split)
##       color   splits as  LLLLLRR,   agree=0.661, adj=0.040, (0 split)
## 
## Node number 2: 27894 observations,    complexity param=0.03361554
##   mean=1632.056, MSE=1248233 
##   left son=4 (19820 obs) right son=5 (8074 obs)
##   Primary splits:
##       carat   < 0.625 to the left,  improve=0.66557750, (0 missing)
##       y       < 5.535 to the left,  improve=0.66520910, (0 missing)
##       x       < 5.495 to the left,  improve=0.66238230, (0 missing)
##       z       < 3.375 to the left,  improve=0.66219560, (0 missing)
##       clarity splits as  RRRLLLLL,  improve=0.01216575, (0 missing)
##   Surrogate splits:
##       x       < 5.485 to the left,  agree=0.993, adj=0.976, (0 split)
##       y       < 5.515 to the left,  agree=0.991, adj=0.969, (0 split)
##       z       < 3.385 to the left,  agree=0.991, adj=0.969, (0 split)
##       clarity splits as  RRLLLLLL,  agree=0.725, adj=0.051, (0 split)
##       cut     splits as  RLLLL,     agree=0.719, adj=0.028, (0 split)
## 
## Node number 3: 15249 observations,    complexity param=0.1851797
##   mean=8158.389, MSE=1.538716e+07 
##   left son=6 (10296 obs) right son=7 (4953 obs)
##   Primary splits:
##       y       < 7.195 to the left,  improve=0.54407490, (0 missing)
##       x       < 7.195 to the left,  improve=0.53738530, (0 missing)
##       carat   < 1.495 to the left,  improve=0.53615420, (0 missing)
##       z       < 4.435 to the left,  improve=0.52631830, (0 missing)
##       clarity splits as  LLLRRRRR,  improve=0.05470891, (0 missing)
##   Surrogate splits:
##       x     < 7.185 to the left,  agree=0.984, adj=0.950, (0 split)
##       carat < 1.455 to the left,  agree=0.980, adj=0.937, (0 split)
##       z     < 4.435 to the left,  agree=0.964, adj=0.890, (0 split)
##       color splits as  LLLLLLR,   agree=0.679, adj=0.012, (0 split)
## 
## Node number 4: 19820 observations
##   mean=1050.302, MSE=264735.3 
## 
## Node number 5: 8074 observations
##   mean=3060.143, MSE=792289.4 
## 
## Node number 6: 10296 observations,    complexity param=0.02661853
##   mean=6151.569, MSE=4751812 
##   left son=12 (7824 obs) right son=13 (2472 obs)
##   Primary splits:
##       clarity splits as  LLLLRRRR,  improve=0.3750779, (0 missing)
##       y       < 6.785 to the left,  improve=0.1210156, (0 missing)
##       carat   < 1.195 to the left,  improve=0.1066457, (0 missing)
##       color   splits as  RRRRLLL,   improve=0.1054268, (0 missing)
##       x       < 6.805 to the left,  improve=0.1042365, (0 missing)
##   Surrogate splits:
##       y < 4.175 to the right, agree=0.76, adj=0.001, (0 split)
## 
## Node number 7: 4953 observations,    complexity param=0.0258319
##   mean=12330.05, MSE=1.172079e+07 
##   left son=14 (3223 obs) right son=15 (1730 obs)
##   Primary splits:
##       y       < 7.855 to the left,  improve=0.30675820, (0 missing)
##       x       < 7.845 to the left,  improve=0.30102780, (0 missing)
##       carat   < 1.915 to the left,  improve=0.29420540, (0 missing)
##       z       < 4.805 to the left,  improve=0.27774230, (0 missing)
##       clarity splits as  LRRRRRRR,  improve=0.07274662, (0 missing)
##   Surrogate splits:
##       x       < 7.905 to the left,  agree=0.983, adj=0.952, (0 split)
##       carat   < 1.895 to the left,  agree=0.974, adj=0.924, (0 split)
##       z       < 4.825 to the left,  agree=0.952, adj=0.862, (0 split)
##       clarity splits as  RRLLLLLL,  agree=0.675, adj=0.068, (0 split)
## 
## Node number 12: 7824 observations
##   mean=5401.156, MSE=2073671 
## 
## Node number 13: 2472 observations
##   mean=8526.661, MSE=5804893 
## 
## Node number 14: 3223 observations
##   mean=10940.83, MSE=8720276 
## 
## Node number 15: 1730 observations
##   mean=14918.16, MSE=7016977
fancyRpartPlot(diamond_model)

With decision tree we can test our model with test data.

pred_Diamond_test <- predict(diamond_model, newdata = diamonds_test)
head(pred_Diamond_test,10)
##         1         2         3         4         5         6         7         8 
## 10940.832  1050.302  5401.156  5401.156  3060.143  3060.143  5401.156  5401.156 
##         9        10 
##  1050.302  5401.156

Before we prune the tree, display CP table look for the lowest cross-validation error(xerror). Lowest xerror is 0.12015 at CP value of 0.01

printcp(diamond_model)
## 
## Regression tree:
## rpart(formula = price ~ x + y + z + carat + cut + color + clarity, 
##     data = diamonds_train)
## 
## Variables actually used in tree construction:
## [1] carat   clarity y      
## 
## Root node error: 6.8939e+11/43143 = 15979193
## 
## n= 43143 
## 
##         CP nsplit rel error  xerror      xstd
## 1 0.609137      0   1.00000 1.00006 0.0098176
## 2 0.185180      1   0.39086 0.39089 0.0043725
## 3 0.033616      2   0.20568 0.20574 0.0022955
## 4 0.026619      3   0.17207 0.17221 0.0022922
## 5 0.025832      4   0.14545 0.14130 0.0019138
## 6 0.010000      5   0.11962 0.12023 0.0017294
#Get the lowest CP value from CP table
min.xerror <- diamond_model$cptable[which.min(diamond_model$cptable[,"xerror"]),"CP"]

min.xerror
## [1] 0.01

Next, we prune the tree based on this value of CP:

# Prune the tree
diamond_model.pruned <- prune(diamond_model, cp = min.xerror) 

# Draw the prune tree
fancyRpartPlot(diamond_model.pruned)

Then use this prune tree to evaluate the our test data

pred_Diamond_test.pruned <- predict(diamond_model.pruned, newdata = diamonds_test)

Obtain the pseudo R2 - a correlation.

fitcorr <- format(cor(diamonds_test$price, pred_Diamond_test.pruned)^2, digits=4)

fitcorr
## [1] "0.8786"

Result

# Create a data frame with the predictions for each method
all.predictions <- data.frame(actual = diamonds_test$price,
                              full.tree = pred_Diamond_test,
                              pruned.tree = pred_Diamond_test.pruned)

#For each actual create model and predictions row
all.predictions <- gather(all.predictions, key = model, value = predictions, 2:3)

# Plot "Predicted vs. actual, by model""
ggplot(data = all.predictions, aes(x = actual, y = predictions)) + 
  geom_point(colour = "blue") + 
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  facet_wrap(~ model, ncol = 2) + 
  ggtitle("Predicted vs. Actual, by model")