在前面的章节中,我们已经发现了钻石质量与价格间这种令人惊讶的关系:质量差的钻石(切工差、颜色差、纯净度低)具有更高的价格:
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(肉眼可见内含物)。
质量差的钻石似乎价格更高,造成这一现象的原因是一个重要的混淆变量:钻石的重量(carat)。重量是确定钻石价格的单一因素中最重要的一个,而质量差的钻石往往更重一些:
ggplot(diamonds, aes(carat, price)) +
geom_hex(bins = 50)
通过拟合一个模型来分离出 carat 变量的作用,我们可以更容易看到钻石的其他特性对price的影响。但是,我们需要先对钻石数据集进行一些调整,以便其更容易处理。
重点关注小于 2.5 克拉的那些钻石(全部数据的 99.7%)。
对重量和价格变量进行对数转换。
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 倍。
如果愿意的话,我们可以继续构建模型,用模型明确表示观察到的效果。例如,我们可以在模型中包括 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
这个结果没什么大价值,但或许我们可以花点时间思考一下出现异常值是因为模型有问题,还是数据中有错误。如果是数据中的错误,那么我们就有机会买到那些错误定了低价的钻石。