Why are low quality diamonds more expensive?

suppressPackageStartupMessages(library("tidyverse"))
suppressPackageStartupMessages(library("modelr"))
suppressPackageStartupMessages(library("lubridate"))
suppressPackageStartupMessages(library("broom"))
suppressPackageStartupMessages(library("nycflights13"))

The splines package is needed for the ns() function used in one of the solutions.

suppressPackageStartupMessages(library("splines"))
options(na.action = na.warn)

This code appears in the section and is necessary for the exercises.

diamonds2 <- diamonds %>%
  filter(carat <= 2.5) %>%
  mutate(
    lprice = log2(price),
    lcarat = log2(carat)
  )

mod_diamond2 <- lm(lprice ~ lcarat + color + cut + clarity, data = diamonds2)

diamonds2 <- add_residuals(diamonds2, mod_diamond2, "lresid2")

1. In the plot of lcarat vs. lprice, there are some bright vertical strips. What do they represent?

The distribution of diamonds has more diamonds at round or otherwise human-friendly numbers (fractions).

2 If log(price) = a_0 + a_1 * log(carat), what does that say about the relationship between price and carat?

Following the examples in the lecture, I use a base-2 logarithm.

mod_log <- lm(log2(price) ~ log2(carat), data = diamonds)
mod_log

Call:
lm(formula = log2(price) ~ log2(carat), data = diamonds)

Coefficients:
(Intercept)  log2(carat)  
     12.189        1.676  

The estimated relationship between carat and price looks like this.

tibble(carat = seq(0.25, 5, by = 0.25)) %>%
  add_predictions(mod_log) %>%
  ggplot(aes(x = carat, y = 2^pred)) +
  geom_line() +
  labs(x = "carat", y = "price")

The plot shows that the estimated relationship between carat and price is not linear. The exact relationship in this model is if \(x\) increases \(r\) times, then \(y\) increases \(r^{a_1}\) times. For example, a two times increase in carat is associated with the following increase in price:

2^coef(mod_log)[2]

log2(carat) 
   3.195002 

Let’s confirm this relationship by checking it for a few values of the carat variable. Let’s increase carat from 1 to 2.

2^(predict(mod_log, newdata = tibble(carat = 2)) -
  predict(mod_log, newdata = tibble(carat = 1)))

       1 
3.195002 

Note that, since predict() predicts log2(carat) rather than carat, the prediction is exponentiated by 2. Now let’s increase carat from 4 to 2.

2^(predict(mod_log, newdata = tibble(carat = 4)) -
  predict(mod_log, newdata = tibble(carat = 2)))

       1 
3.195002 

Finally, let’s increase carat from 0.5 to 1.

2^(predict(mod_log, newdata = tibble(carat = 1)) -
  predict(mod_log, newdata = tibble(carat = 0.5)))

       1 
3.195002 

All of these examples return the same value, \(2^{a_1}=3.2\).

So why is this? Let’s ignore the names of the variables in this case and consider the equation: \(\log_by = a_0 +a_1\log_bx\)

We want to understand how the difference in \(y\) is related to the difference in \(x\). Now, consider this equation at two different values

\(\log_by_0 = a_0 +a_1\log_bx_0\) \(\log_by_1 = a_0 +a_1\log_bx_1\)

What is the value of the difference, \(\log y_1 - \log y_0\) ?

\(\log_b\left(y_1\right)-\log_b\left(y_0\right) = \left( a_0 + a1\log_bx_1\right)- \left( a_0 + a1\log_bx_0\right)\)

\(=a_1\left( \log_bx_1- \log_bx_0\right),\)

\(\log_b\left(\frac{y_1}{y_0}\right)=\log_b\left(\frac{x_1}{x_0}\right)^{a_1}\)

\(\frac{y_1}{y_0}=\left(\frac{x_1}{x_0}\right)^{a_1}\)

Let \(s = \frac{y_1}{y_0}\) and \(r=\frac{x_1}{x_0}\). Then,

\(s = r^{a_1}\).

In other words, an \(r\) times increase in \(x\), is associated with a \(r^{a_1}\) times increase in \(y\). Note that this relationship does not depend on the base of the logarithm, \(b\).

There is another approximation that is commonly used when logarithms appear in regressions. The first way to show this is using the approximation that \(x\) is small, meaning that \(x\approx0\),

\(\log \left( 1+x \right)\approx x\)

This approximation is the first order Taylor expansion of the function at \(x=0\). Now consider the relationship between the percent change in \(x\) and the percent change in \(y\),

\(\log\left(y+\Delta y\right)− \log y=\left(\alpha+\beta\log\left(x+\Delta x\right)\right)−\left(\alpha+\beta \log x\right)\) \(\log\left(\frac{y+\Delta y}{y}\right)=\beta\log\left(\frac{x+\Delta x}{x}\right)\) \(\log\left(1+\frac{\Delta y}{y}\right)=\beta\log\left(1+\frac{\Delta x}{x}\right)\) \(\frac{\Delta y}{y} \approx \beta\left(\frac{\Delta x}{x}\right)\)

Thus a 1% percentage change in \(x\) is associated with a \(\beta\) percent change in \(y\).

This relationship can also be derived by taking the derivative of \(\log y\) with respect to
\(x\). First, rewrite the equation in terms of \(y\),

\(y=\exp(a_0+a_1 \log(x))\)

Then differentiate \(y\) with respect to \(x\),

\(dy=\exp(a_0+a_1 \log x)\left(\frac{a_1}{x}\right)dx\) \(=a_1y\left(\frac{dx}{x}\right)\) \((dy/y)=a_1(dx/x)\) \(\% \Delta y = a_1 \% \Delta x\)

3. Extract the diamonds that have very high and very low residuals. Is there anything unusual about these diamonds? Are they particularly bad or good, or do you think these are pricing errors?

diamonds2 %>%
  filter(abs(lresid2) > 1) %>%
  add_predictions(mod_diamond2) %>%
  mutate(pred = round(2^pred)) %>%
  select(price, pred, carat:table, x:z) %>%
  arrange(price)

I did not see anything too unusual. Do you?

4. Does the final model, mod_diamonds2, do a good job of predicting diamond prices? Would you trust it to tell you how much to spend if you were buying a diamond?

Plotting the residuals of the model shows that there are some large outliers for small carat sizes. The largest of these residuals are a little over two, which means that the true value was four times lower. Most of the mass of the residuals is between -0.5 and 0.5, which corresponds to about \(\pm 40\). There seems to be a slight downward bias in the residuals as carat size increases.

ggplot(diamonds2, aes(lcarat, lresid2)) +
  geom_hex(bins = 50)

lresid2_summary <- summarise(diamonds2,
  rmse = sqrt(mean(lresid2^2)),
  mae = mean(abs(lresid2)),
  p025 = quantile(lresid2, 0.025),
  p975 = quantile(lresid2, 0.975)
)
lresid2_summary

While in some cases the model can be wrong, overall the model seems to perform well. The root mean squared error is 0.19 meaning that the average error is about -14%. Another summary statistics of errors is the mean absolute error (MAE), which is the mean of the absolute values of the errors. The MAE is 0.15, which is -11%. Finally, 95% of the residuals are between -0.37 and 0.38, which correspond to 23—31.

Whether you think that this is a good model depends on factors outside the statistical model itself. It will depend on the how the model is being used. I have no idea how to price diamonds, so this would be useful to me in order to understand a reasonable price range for a diamond, so I don’t get ripped off. However, if I were buying and selling diamonds as a business, I would probably require a better model.