18.2 为什么质量差的钻石更贵

在前面的章节中,我们已经发现了钻石质量与价格间这种令人惊讶的关系:质量差的钻石(切工差、颜色差、纯净度低)具有更高的价格:

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(modelr)
options(na.action = na.warn)
ggplot(diamonds, aes(cut, price)) + geom_boxplot()

ggplot(diamonds, aes(color, price)) + geom_boxplot()

ggplot(diamonds, aes(clarity, price)) + geom_boxplot()

注意,最差的钻石颜色是 J(微黄),最差的纯净度是 I1(肉眼可见内含物)。

18.2.1 价格与重量

质量差的钻石似乎价格更高,造成这一现象的原因是一个重要的混淆变量:钻石的重量(carat)。重量是确定钻石价格的单一因素中最重要的一个,而质量差的钻石往往更重一些:

ggplot(diamonds, aes(carat, price)) +
 geom_hex(bins = 50)

通过拟合一个模型来分离出 carat 变量的作用,我们可以更容易看到钻石的其他特性对price的影响。但是,我们需要先对钻石数据集进行一些调整,以便其更容易处理。

  1. 重点关注小于 2.5 克拉的那些钻石(全部数据的 99.7%)。

  2. 对重量和价格变量进行对数转换。

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

这两个调整可以让我们更轻松地看到 carat 和 price 之间的关系:

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

对数转换在这个示例中非常有用,因为它可以让模式变为线性的,而线性模式是最容易处理的。现在我们进行下一步,从数据中去除这种强烈的线性模式。我们通过拟合一个模型让这种模式成为显式的:

mod_diamond <- lm(lprice ~ lcarat, data = diamonds2)

接着我们检查模型,看看它能够反映出数据中的哪些信息。注意,因为我们对预测值进行了反向变换,还原了对数转换,所以可以将预测值覆盖在原始数据上:

grid <- diamonds2 %>%
 data_grid(carat = seq_range(carat, 20)) %>%
 mutate(lcarat = log2(carat)) %>%
 add_predictions(mod_diamond, "lprice") %>%
 mutate(price = 2 ^ lprice)
ggplot(diamonds2, aes(carat, price)) +
 geom_hex(bins = 50) +
 geom_line(data = grid, color = "red", size = 1)

这张图可以告诉我们关于这份数据的一些有趣信息。如果我们相信这个模型,那么大钻石要比预料中便宜得多。这可能是因为数据集中没有价格超过 $19 000 的钻石。

现在我们可以检查一下残差,它可以用来验证我们是否成功移除了强烈的线性模式:

diamonds2 <- diamonds2 %>%
 add_residuals(mod_diamond, "lresid")
ggplot(diamonds2, aes(lcarat, lresid)) +
 geom_hex(bins = 50)

重要的是,我们现在可以使用残差代替 price 来重新绘图了:

ggplot(diamonds2, aes(cut, lresid)) + geom_boxplot()

ggplot(diamonds2, aes(color, lresid)) + geom_boxplot()

ggplot(diamonds2, aes(clarity, lresid)) + geom_boxplot()

现在我们可以看到期望中的关系了:当钻石的质量下降时,其相应价格也随之下降。为了解释y轴,我们需要思考一下残差的意义及其使用的标度。残差为 -1 表示 lprice 比仅使用重量进行估计的预测值少一个单位。2^-1 就是 1/2,因此值为 -1的点的价格为预计价格的一半,残差为 1 时,价格则是预计价格的 2 倍。

18.2.2 一个更复杂的模型

如果愿意的话,我们可以继续构建模型,用模型明确表示观察到的效果。例如,我们可以在模型中包括 color、cut 和 clarity 变量,以将这 3 个分类变量的效果明确表示出来:

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

现在模型中包括了 4 个预测变量,因此更加难以进行可视化。好在这些变量还是彼此独立的,这意味着我们可以在 4张图中分别绘制出它们。为了让这个过程更简单一些,我们在data_grid() 函数中使用 .model 参数:

grid <- diamonds2 %>%
 data_grid(cut, .model = mod_diamond2) %>%
 add_predictions(mod_diamond2)
grid
## # A tibble: 5 x 5
##   cut       lcarat color clarity  pred
##   <ord>      <dbl> <chr> <chr>   <dbl>
## 1 Fair      -0.515 G     VS2      11.2
## 2 Good      -0.515 G     VS2      11.3
## 3 Very Good -0.515 G     VS2      11.4
## 4 Premium   -0.515 G     VS2      11.4
## 5 Ideal     -0.515 G     VS2      11.4
ggplot(grid, aes(cut, pred)) +
 geom_point()

如果模型需要你还没有明确提供的变量,data_grid()函数会自动使用“典型”值来填充它们。对于连续变量,模型使用中位数;对于分类变量,模型使用最常见的值(或多个值,如果有同样数量的多个值的话):

diamonds2 <- diamonds2 %>%
 add_residuals(mod_diamond2, "lresid2")
ggplot(diamonds2, aes(lcarat, lresid2)) +
 geom_hex(bins = 50)

这张图说明一些钻石有非常大的残差。记住,残差为 2 表示钻石的价格是预计价格的 4倍。通常还应该检查一下异常值:

diamonds2 %>%
 filter(abs(lresid2) > 1) %>%
 add_predictions(mod_diamond2) %>%
 mutate(pred = round(2 ^ pred)) %>%
 select(price, pred, carat:table, x:z) %>%
 arrange(price)
## # A tibble: 16 x 11
##    price  pred carat cut       color clarity depth table     x     y     z
##    <int> <dbl> <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  1013   264  0.25 Fair      F     SI2      54.4    64  4.3   4.23  2.32
##  2  1186   284  0.25 Premium   G     SI2      59      60  5.33  5.28  3.12
##  3  1186   284  0.25 Premium   G     SI2      58.8    60  5.33  5.28  3.12
##  4  1262  2644  1.03 Fair      E     I1       78.2    54  5.72  5.59  4.42
##  5  1415   639  0.35 Fair      G     VS2      65.9    54  5.57  5.53  3.66
##  6  1415   639  0.35 Fair      G     VS2      65.9    54  5.57  5.53  3.66
##  7  1715   576  0.32 Fair      F     VS2      59.6    60  4.42  4.34  2.61
##  8  1776   412  0.29 Fair      F     SI1      55.8    60  4.48  4.41  2.48
##  9  2160   314  0.34 Fair      F     I1       55.8    62  4.72  4.6   2.6 
## 10  2366   774  0.3  Very Good D     VVS2     60.6    58  4.33  4.35  2.63
## 11  3360  1373  0.51 Premium   F     SI1      62.7    62  5.09  4.96  3.15
## 12  3807  1540  0.61 Good      F     SI2      62.5    65  5.36  5.29  3.33
## 13  3920  1705  0.51 Fair      F     VVS2     65.4    60  4.98  4.9   3.23
## 14  4368  1705  0.51 Fair      F     VVS2     60.7    66  5.21  5.11  3.13
## 15 10011  4048  1.01 Fair      D     SI2      64.6    58  6.25  6.2   4.02
## 16 10470 23622  2.46 Premium   E     SI2      59.7    59  8.82  8.76  5.25

这个结果没什么大价值,但或许我们可以花点时间思考一下出现异常值是因为模型有问题,还是数据中有错误。如果是数据中的错误,那么我们就有机会买到那些错误定了低价的钻石。