Assignment 3: Diamonds Price Estimation

Introduction

This analysis focuses on the relationship between price and quality information of diamonds. The main objective is to create a model which predicts the price of a given diamond. The diamonds dataset is included in the ggplot2 package. Therefore it is necessary to load it.

Loading Libraries

Overview and Manipulation of Data

Summary

Let’s take a quick view of the variables of the diamonds dataset. There are 53940 observations and 10 variables. I also stored the original data in raw_diamonds.

# In case I need the original data
raw_diamonds <- diamonds

# Check the variables and types
str(diamonds)

## tibble [53,940 × 10] (S3: tbl_df/tbl/data.frame)
##  $ carat  : num [1:53940] 0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 2 6 7 7 6 5 2 5 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 2 3 5 4 2 6 7 3 4 5 ...
##  $ depth  : num [1:53940] 61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
##  $ table  : num [1:53940] 55 61 65 58 58 57 57 55 61 61 ...
##  $ price  : int [1:53940] 326 326 327 334 335 336 336 337 337 338 ...
##  $ x      : num [1:53940] 3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
##  $ y      : num [1:53940] 3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
##  $ z      : num [1:53940] 2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

Variables

• carat : weight of the diamond between (0.2 ,5.01)
• cut : quality of the cut in (Fair, Good, Very Good, Premium, Ideal) from worst to best
• color : diamond color in (D, E, F, G, H, I, J) from worst to best
• clarity : clarity of the diamond in (I1, SI2, SI1, VS2, VS1, VVS2, VVS1, IF) from worst to best
• depth : total depth percentage = z / mean(x, y) between (43 ,79).
• table : width of top of diamond relative to widest point between (43 ,95).
• price : price in US dollars between (326 ,18823)
• x : length in mm between (0 ,10.74)
• y : width in mm between (0 ,58.9)
• z : depth in mm between (0 ,31.8).

The variables cut, color and clarity are ordinal categorical variables, which means they are categorical but they can be ordered as factors. Therefore their levels are written from worst to best.
Figure below gives some information on carat, cut, color and clarity.

url <- "https://github.com/pjournal/boun01-ilaydacelenk/raw/gh-pages/Assignment3/cut-color-clarity-carat.png"

Data Exploration

Let’s see some statistics about the dataset diamonds. For each numerical attribute we can compare median and mean. If they are very far away from each other, then they cannot be normally distributed. It is better to have them close. For the categorical attributes, we can see how many diamonds are in each category. We also would want them to close to each other, in order to have a realiable model.

# Median/mean for numerical, count for categorical attributes
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  
## 

Price vs Numerical Attributes

The following plot, shows linear relationships between price and numerical attributes.

# Linearity
diamonds %>% gather(-price, -cut, -color, -clarity, key = "var", value = "value") %>% 
  ggplot(aes(x = value, y = price)) + geom_point() + geom_smooth(method='lm', formula= y~x) + facet_wrap(~ var, scales = "free") +  theme_bw()

Log(Price) vs Numerical Attributes

The following plot, shows linear relationships between log(price) and numerical attributes.

# Linearity
diamonds %>% gather(-price, -cut, -color, -clarity, key = "var", value = "value") %>% 
  ggplot(aes(x = value, y = log(price))) + geom_point() + geom_smooth(method='lm', formula= y~x) + facet_wrap(~ var, scales = "free") +  theme_bw()

Correlation

This is the correlation matrix of numerical variables with significance levels according to the Pearson method. Since our objective is predicting the price, we focus on the relationship with price. For example, from the heatmap we can say that there is a significant relationship between price and carat, but price does not have a significant relationship with table and depth, therefore they are crossed on the map. Carat, x, y and z variables are highly correlated. This was expected since carat is the weight of the diamond and it has a positive relationship with volume which is described by x, y and z, which are length, width and depth. Above 0.8 means high correlation and below -0.8 means high anti-correlation.

# Select numerical attributes
diamonds_numerical <- diamonds %>% select(-cut, -color, -clarity)

# Plot correlation matrix for numerical attributes
corr <- round(cor(diamonds_numerical), 1) #correlation matrix
p_values_matrix <- cor_pmat(corr) #p-values to observe the significance of the variables
ggcorrplot(corr, hc.order = TRUE, lab = TRUE, p.mat = p_values_matrix)

Distribution of Prices

The figure below shows that distribution of prices is right skewed. Since there are a lot of large values, we can take a look at logarithm of th price column. When we do that, we observe log(price) is more suitable since it is closer to being normal. Therefore, replacing price column with log(price) makes sense for the regression model.

# Right skewed
ggplot(diamonds, aes(x=price)) + geom_histogram(binwidth=500)

# More normal
ggplot(diamonds, aes(x=log(price))) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Use log(price) instead of price
diamonds$log.price = log(diamonds$price)
diamonds <- diamonds %>% select(-price)

Cut Proportions

# Count the duplicates
table(diamonds$cut)

## 
##      Fair      Good Very Good   Premium     Ideal 
##      1610      4906     12082     13791     21551

# Get proportions
plot_cut <- diamonds %>% group_by(cut) %>% summarise(Count = n()) %>% mutate(Proportion=round(Count/sum(Count),1))

Out of 53940 different diamonds in diamonds dataset, 21551 are classifies as Ideal, which corresponds to 0.3995365.

# Show proportions on a pie plot
ggplot(plot_cut, aes(x = "", y = Proportion, fill = cut)) + geom_bar(stat="identity", width=1) + coord_polar("y", start=0) + geom_text(aes(label = paste0(round(Proportion*100), "%")), position = position_stack(vjust = 0.5)) + labs(x = NULL, y = NULL, fill = NULL, title = "Proportion of Cut Types") + theme_classic() + theme(axis.line = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(hjust = 0.5))

Distribution of Price for Different Cuts

Even though proportions of cut types are not close to each other, using log(price) instead of price is more preferable since price has normal distribution in each category. Using box plot, we see the distribution of log(price) for different cut types. The plot below shows the importance of cut types. If we plotted price instead of log(price) then the distribution would be very skewed.

diamonds %>% ggplot(aes(x = cut, y = log.price, fill=cut)) + geom_boxplot() + 
  xlab("Cut Types") + ylab("log(price)") + ggtitle("Cut vs log(price)")

Supervised Learning

Since we are trying to predict labeled column values, this is a supervised learning.

Train Test Split

Sometimes machine learning models may memorize instead of learning. This is called overfitting. In order to understand such situation, we split the data into test and train data sets. We fit the model using train, predict and evaluate accuracy using test data. If the train prediction accuracy is much higher from test prediction accuracy then overfitting may have occurred and we may want to change the parameters of the model or use a different model.

For train data, we take a sample of 80% of the data without replacement and the rest 20% will be the test data.

# Split
train_indices <- sample(seq_len(nrow(diamonds)), size = 0.8*nrow(diamonds), replace = FALSE)
train <- diamonds[train_indices, ]
test <- diamonds[-train_indices, ]

Model: Multiple Linear Regression

According to p-values for t-tests and F-test, all of the attributes except have significant effect on log(price).

linear_model1<-lm(log.price ~ ., data=train)
summary(linear_model1)

## 
## Call:
## lm(formula = log.price ~ ., data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2638 -0.0911  0.0015  0.0902  9.0693 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -2.5579352  0.0685251  -37.328  < 2e-16 ***
## carat       -0.6093868  0.0084957  -71.729  < 2e-16 ***
## cut.L        0.1068633  0.0039515   27.044  < 2e-16 ***
## cut.Q       -0.0383897  0.0031574  -12.158  < 2e-16 ***
## cut.C        0.0382466  0.0027078   14.125  < 2e-16 ***
## cut^4        0.0116628  0.0021604    5.399 6.75e-08 ***
## color.L     -0.4535165  0.0030227 -150.039  < 2e-16 ***
## color.Q     -0.1041008  0.0027459  -37.911  < 2e-16 ***
## color.C     -0.0119221  0.0025639   -4.650 3.33e-06 ***
## color^4      0.0170633  0.0023579    7.237 4.68e-13 ***
## color^5     -0.0082932  0.0022247   -3.728 0.000193 ***
## color^6      0.0009673  0.0020218    0.478 0.632350    
## clarity.L    0.8912203  0.0052705  169.096  < 2e-16 ***
## clarity.Q   -0.2586442  0.0049184  -52.587  < 2e-16 ***
## clarity.C    0.1440099  0.0042061   34.238  < 2e-16 ***
## clarity^4   -0.0673837  0.0033572  -20.071  < 2e-16 ***
## clarity^5    0.0274296  0.0027393   10.013  < 2e-16 ***
## clarity^6   -0.0053161  0.0023854   -2.229 0.025842 *  
## clarity^7    0.0284388  0.0021078   13.492  < 2e-16 ***
## depth        0.0524856  0.0007773   67.524  < 2e-16 ***
## table        0.0088021  0.0005084   17.315  < 2e-16 ***
## x            1.1697516  0.0053947  216.832  < 2e-16 ***
## y            0.0292519  0.0030215    9.681  < 2e-16 ***
## z            0.0378274  0.0052685    7.180 7.09e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.176 on 43128 degrees of freedom
## Multiple R-squared:   0.97,  Adjusted R-squared:   0.97 
## F-statistic: 6.056e+04 on 23 and 43128 DF,  p-value: < 2.2e-16

# Remove the one with highest p-value 
linear_model2<-lm(log.price ~ . -color^6, data=train)
summary(linear_model2)

## 
## Call:
## lm(formula = log.price ~ . - color^6, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2822 -0.1232  0.0088  0.1281  9.3911 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.2674415  0.0854524 -26.535  < 2e-16 ***
## carat       -0.7837840  0.0104911 -74.709  < 2e-16 ***
## cut.L        0.1000367  0.0049284  20.298  < 2e-16 ***
## cut.Q       -0.0365415  0.0039388  -9.277  < 2e-16 ***
## cut.C        0.0416451  0.0033762  12.335  < 2e-16 ***
## cut^4        0.0136627  0.0026950   5.070 4.00e-07 ***
## clarity.L    0.8346509  0.0065322 127.774  < 2e-16 ***
## clarity.Q   -0.2557648  0.0061263 -41.749  < 2e-16 ***
## clarity.C    0.1409132  0.0052450  26.866  < 2e-16 ***
## clarity^4   -0.0846360  0.0041839 -20.229  < 2e-16 ***
## clarity^5    0.0353468  0.0034149  10.351  < 2e-16 ***
## clarity^6    0.0006886  0.0029752   0.231    0.817    
## clarity^7    0.0484121  0.0026244  18.447  < 2e-16 ***
## depth        0.0483248  0.0009689  49.877  < 2e-16 ***
## table        0.0082405  0.0006340  12.997  < 2e-16 ***
## x            1.1987336  0.0067195 178.397  < 2e-16 ***
## y            0.0303250  0.0037696   8.045 8.87e-16 ***
## z            0.0405558  0.0065725   6.171 6.87e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2196 on 43134 degrees of freedom
## Multiple R-squared:  0.9532, Adjusted R-squared:  0.9532 
## F-statistic: 5.173e+04 on 17 and 43134 DF,  p-value: < 2.2e-16

# Now all attributes are significant
linear_model3<-lm(log.price ~ . -color^6 -clarity^6, data=train)
summary(linear_model3)

## 
## Call:
## lm(formula = log.price ~ . - color^6 - clarity^6, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1031 -0.1668 -0.0031  0.1603  8.8468 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.9506774  0.1066354  -8.915  < 2e-16 ***
## carat       -0.7022745  0.0130724 -53.722  < 2e-16 ***
## cut.L        0.1830228  0.0061224  29.894  < 2e-16 ***
## cut.Q       -0.0620637  0.0048934 -12.683  < 2e-16 ***
## cut.C        0.0719969  0.0042204  17.059  < 2e-16 ***
## cut^4        0.0289130  0.0033786   8.558  < 2e-16 ***
## depth        0.0370300  0.0012116  30.562  < 2e-16 ***
## table        0.0049572  0.0007955   6.231 4.67e-10 ***
## x            1.0966612  0.0083785 130.891  < 2e-16 ***
## y            0.0336991  0.0047332   7.120 1.10e-12 ***
## z            0.0550543  0.0082517   6.672 2.56e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2758 on 43141 degrees of freedom
## Multiple R-squared:  0.9263, Adjusted R-squared:  0.9262 
## F-statistic: 5.419e+04 on 10 and 43141 DF,  p-value: < 2.2e-16

Prediction

Let’s get the predicted log(price) values for the test set in order to compute residuals.

predict_test = predict(object=linear_model1, newdata=test)

Mean Squared Error

We can compare different models by their mean squared error(MSE). The model’s MSE is 0.0302056.

mse1 <- mean((predict_test-test$log.price)^2)
mse1

## [1] 0.03020564

Possible Problems Related to Assumptions of Regression Analysis

In order to carry out a regression analysis, we need to check the assumptions. Otherwise, the hypothesis tests are not going to be valid.

• Normality of Error Terms
• Linearity
• Independence of Error Terms
• Constant Variance of Error Terms (Homoscedasticity)

Normality can be observed by histogram of the error distribution, for the other three we can carry out a residual analysis.

Normality of Error Terms

As we have seen above, the distribution of prices are right skewed. Let’s check how its logarithm behaves. The figure below is the best we can get in terms of normality.

ggplot(test, aes(x=log.price)) + geom_histogram()

Linearity

The following graph shows residuals’ behavior according to changing log(price) values. By ignoring a couple of outliers, linearity is not violated since we can draw parallel lines where the points between would be evenly distributed.

test$residuals = test$log.price - predict_test
test$residuals = round(test$residuals, digits=2)
ggplot(test) + 
  geom_point(aes(x = log.price, y = residuals)) +
  geom_hline(yintercept=0,col="red") 

Independence of Error Terms

Same graph shows that error terms are not switching signs by different log(price) values and the error terms are independent.

Constant Variance of Error Terms (Homoscedasticity)

According to Residuals vs log(price) graph, error variance is not changing due to logarithm of prices. Therefore, this assumption is also not violated.

Results

• Regression model for diamonds dataset works better for log(price) instead of price.
• Assumptions for regression are shown to be valid by residual analysis.
• Predictors are highly significant.

References

• Assignment page
• Inserting an image
• Train/Test Split

To see all my work related to IE 48A course, visit my Progress Journal website. You can also visit my personal Github Website.