简介:《钻石价格分析-基于R语言》是一份面向R语言初学者与数据分析入门者的实战型文档,围绕ggplot2内置diamonds数据集,完整呈现从数据清洗到回归建模的分析路径。内容涵盖缺失值与重复值检查、异常值箱线图识别与剔除、价格分布直方图与核密度解读,以及切工、颜色、净度、克拉重量等变量与价格关系的可视化探索,并进一步完成定性变量因子化、新变量构造、简单线性回归和多元回归预测,保留建模流程、图表解读与结果解释的完整脉络。资源包共1个docx文件,约2.4MB,便于集中阅读和按章节复现实验。已有285人学习下载,适合希望掌握R语言数据处理、回归分析与机器学习基础流程的读者,可据此理解钻石定价的主要影响因素,积累一份可迁移到其他数据集的分析模板。
1. 一颗 0.9 克拉的钻石凭什么标到 8000 美元
同一颗 0.9 克拉的石头,切工评级降一档,报价能差出两千美元;颜色再掉两级,价格又被砍掉一截。想弄清楚这笔钱到底花在哪个属性上,最省事的入口就是 ggplot2 包里那份自带的 diamonds 数据——53940 行、10 个字段。这套 R 语言钻石价格分析项目,做的就是把这份数据从装载、清洗、描述性统计一路推到多元线性回归:先看价格怎么分布,再用克拉、切工、颜色、净度去解释价格,最后落到预测和显著性检验。
它不需要爬虫、不需要数据库,一台装了 R 的笔记本就能从头跑到尾。刚接触 R 数据分析、想完整走一遍回归分析流程的人,拿它练手最合适;已经会写lm()但说不清哑变量基准组、交互项和逐步回归怎么落地的人,这份脚本也能当对照物看。
2. 数据装载与清洗:把 diamonds 从 ggplot2 里挖出来
2.1 依赖清单与脚本骨架
这份分析对包的依赖其实很轻,全部围绕「数据整形 + 画图 + 线性模型 + 诊断」四件事展开。ggplot2 既是数据源又是绘图主力,dplyr 负责分组聚合,mice 用来看缺失值模式,car 提供 VIF,MASS 提供stepAIC,lmtest 做异方差检验,gridExtra 把多张图拼到一起。一次性装齐,后面就不用反复中断流程。
# 依赖安装:先算差集,装了的不再重复装 pkgs <- c("ggplot2", "dplyr", "tidyr", "mice", "car", "MASS", "lmtest", "gridExtra") new_pkgs <- setdiff(pkgs, rownames(installed.packages())) if (length(new_pkgs) > 0) install.packages(new_pkgs) invisible(lapply(pkgs, library, character.only = TRUE))setdiff那一行是给自己省时间的,installed.packages()会返回当前环境里所有已装包的矩阵,取行名做差集即可。invisible包一层是为了让lapply加载包时不要往控制台糊一堆 “Attaching package” 的提示。
我习惯把整条链路拆成四个脚本,中间结果落盘成 rds,改模型的时候不用每次从头清洗。目录大概长这样:
diamond-price-analysis/ ├── 01_load.R # 装载 diamonds、体检、清洗、落盘 ├── 02_eda.R # 描述性统计与全套图形 ├── 03_model.R # 简单回归 / 全模型 / 交互项 / 逐步回归 ├── 04_diagnose.R # 假设检验与残差诊断 ├── data/clean_diamonds.rds └── output/figs/拆脚本的好处是调试成本低。改一个交互项写法,只要重跑03_model.R,清洗那一步读 rds 就行,省下来的时间够你多试三组变量组合。
2.2 数据装载与三轮体检
diamonds 是 ggplot2 的内置数据集,data()一句就出来了,不用read.csv,也没有编码坑。装载完的第一件事不是画图,是把数据的形状摸清楚:多少行、哪些列是数值、哪些是字符,各列的量纲差多少。
data(diamonds, package = "ggplot2") dim(diamonds) # 53940 行 10 列 str(diamonds) # carat/price 等为 num,cut/color/clarity 为 Ord.factor summary(diamonds) # 一眼看清每列的最小值、中位数、最大值str()那一步会暴露出一个重要信息:cut、color、clarity 三列在原始数据里已经是Ord.factor(有序因子)。这一点后面建模会专门处理,先记住。
三轮体检分别是缺失值、重复值、取值越界。
# 第一轮:缺失值 colSums(is.na(diamonds)) mice::md.pattern(diamonds, plot = FALSE) # 第二轮:完全重复行 sum(duplicated(diamonds)) # 第三轮:取值范围是否离谱 sapply(diamonds[, c("carat", "depth", "table", "x", "y", "z", "price")], range)colSums(is.na())返回全 0,说明这份数据本身是完整的,没有缺失值需要插补,md.pattern只是再确认一遍缺失模式。duplicated返回 0,也没有完全重复的行——但要注意,没有重复行不等于没有「近重复」,同一组参数不同价格的记录在这份数据里是正常存在的,别手贱去按除 price 之外的列去重,那会把真实的市场波动删掉。
真正值得警惕的是range()的输出。x、y、z 的最小值全是 0。一颗钻石的长度、宽度、深度不可能为 0,这是录入或测量环节留下的脏记录。原始数据里有 20 条左右这样的行,数量不多,但足以在算体积、画散点的时候把图拉变形。
2.3 异常值清理与特征构造
异常值处理我一般分两个维度做:价格维度按分布尾部裁,尺寸维度按物理意义裁。两者不能混在一起写,出问题的时候不好排查。
# 维度一:价格尾部。用 3 倍 IQR 规则做粗筛 q <- quantile(diamonds$price, probs = c(0.25, 0.75), na.rm = TRUE) iqr <- IQR(diamonds$price) d1 <- subset(diamonds, price >= q[1] - 3 * iqr & price <= q[3] + 3 * iqr) # 维度二:物理尺寸必须为正,否则算不出体积 d2 <- subset(d1, x > 0 & y > 0 & z > 0) nrow(diamonds); nrow(d1); nrow(d2) # 53940 / 53840 上下 / 53820 上下quantile取上下四分位,IQR是四分位距,3 倍 IQR 比 1.5 倍宽松,适合这种本来就右偏的价格分布——用 1.5 倍会砍掉太多真实的高价钻石。你也可以按原文的思路用「三个标准差之外」来剔,效果接近,但对右偏数据会剔得更狠。
清洗完就该造新变量了。原始字段里 x、y、z 是三个独立尺寸,单独放进模型很难解释,合成体积更有物理含义;另外价格和克拉的比值能反映「单位克拉的溢价」,做分群对比时很好用。
diamonds_fe <- d2 %>% mutate( # 体积:比单看 carat 更贴近"看起来大不大" volume = x * y * z, # 单位克拉价格,用于分群对比 price_per_carat = price / carat, # 显式指定因子水平顺序,基准组落在最差档 cut = factor(cut, levels = c("Fair", "Good", "Very Good", "Premium", "Ideal")), color = factor(color, levels = c("J", "I", "H", "G", "F", "E", "D")), clarity = factor(clarity, levels = c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF")) ) saveRDS(diamonds_fe, "data/clean_diamonds.rds")这里有个关键选择:我没有加ordered = TRUE。原因在后面的避坑章会详细说,简单讲就是有序因子在lm()里默认走正交多项式对比,系数会变成一组你没法直接读的线性组合。显式给levels赋值就够了,它决定了谁是基准组——现在基准组是 Fair、J、I1,也就是三个维度里最差的那一档,回归系数解释成「相对于最差档的提升」,最符合直觉。
color 这一列要特别注意:原始 levels 是 D 到 J,D 最好。我这里倒过来写成 J 到 D,是为了让因子水平顺序和「好坏」单调一致,后面画图排序和建有序逻辑都不容易搞混。
3. 探索性分析:价格分布与离散变量的分组对比
3.1 价格分布:直方图叠核密度与正态曲线
描述性分析的第一步是看目标变量的形状,因为形状决定了后面用不用对数变换。把直方图、核密度和正态曲线叠在一张图上,三条线一对比,偏态有多重立刻清楚。
library(ggplot2) ggplot(diamonds_fe, aes(x = price)) + # 直方图纵轴换成密度,才能和密度曲线共用一个坐标系 geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "#4C78A8", colour = "white", alpha = 0.85) + # 核密度:不假设分布形式 geom_density(colour = "#E45756", linewidth = 1) + # 同均值同标准差的正态曲线,作为对照基准 stat_function(fun = dnorm, args = list(mean = mean(diamonds_fe$price), sd = sd(diamonds_fe$price)), colour = "black", linetype = "dashed") + labs(x = "价格(美元)", y = "密度") + theme_minimal()after_stat(density)是 ggplot2 3.4 之后的写法,老版本要写..density..,如果这句报错就是版本问题,升级包或者改写法都行。bins = 60是试出来的,53940 条记录用 30 个箱太粗、100 个箱太碎,60 个刚好能看出断层。
跑出来的图会告诉你三件事:第一,价格分布强烈右偏,绝大多数钻石挤在 0 到 2000 美元之间;第二,核密度曲线在 1000 美元附近有一个明显的峰,说明这个价位是市场最集中的成交区间;第三,和正态曲线比,8000 美元之前两条曲线几乎不重合,8000 之后反而一起贴着横轴走——这对应的是一个残酷的行业事实:高价钻石数量极少,价格越往上,样本越稀疏。
顺带一句,原始 price 的最大值是 18823,按 3 倍 IQR 裁完大致落在 11886 附近,最低 326。这个区间记住,后面解释回归截距为负的时候要用。
3.2 离散变量的分组对比:箱线图、小提琴图与中位数排序
cut、color、clarity 是三个分级变量,看它们和价格的关系,箱线图最直观。但价格跨度大,直接画箱子会被高价样本压扁,所以我习惯对纵轴取 log10。
p_cut <- ggplot(diamonds_fe, aes(x = cut, y = price, fill = cut)) + geom_boxplot(outlier.alpha = 0.12, outlier.size = 0.5) + scale_y_log10() + # 压住长尾,别让箱子挤成一条线 labs(x = "切工", y = "价格(log10 美元)") + theme_minimal() + theme(legend.position = "none") # 分组中位数,比看图更准 diamonds_fe %>% group_by(cut) %>% summarise(n = n(), median = median(price), mean = mean(price)) %>% arrange(desc(median))这里会撞上一个非常反直觉的结果:按中位数排,Fair 组反而排在最前面,Ideal 组排在最后。如果只看这张图就下「切工越差越贵」的结论,就翻车了。
原因不复杂——Fair 切工的钻石里,大克拉的占比高得多;而 Ideal 切工里有大量 0.3 到 0.5 克拉的小石头。克拉这个混杂变量没控制住,任何分组比较都是假的。正确做法是固定克拉区间再比,比如只取 0.7 到 0.9 克拉这一段重新画一遍,结论会立刻翻过来。这也是为什么后面必须上多元回归:单变量对比只能提假设,控制变量之后才谈得上结论。
同理,color 从 J 到 D 画一遍小提琴图,能看到价格中位数随颜色等级提升基本单调上升;clarity 从 I1 到 IF 也是上升趋势,但 I1 到 SI2 这一段不单调,说明净度到优质区间之后,决定价格的维度就不止一个了。
ggplot(diamonds_fe, aes(x = color, y = price, fill = color)) + geom_violin(scale = "count", trim = FALSE) + # scale 控制各组的宽度权重 geom_boxplot(width = 0.12, fill = "white", outlier.shape = NA) + scale_y_log10() + labs(x = "颜色(J 最差 → D 最好)", y = "价格(log10 美元)") + theme_minimal() + theme(legend.position = "none")scale = "count"让各组小提琴的宽度正比于样本量,一眼就能看出 Ideal 和 Premium 占了绝大多数、Fair 只有一千多条。trim = FALSE保留尾部,不然极端值会被削掉,看不到真实分布。
3.3 克拉与价格的散点图:线性之外的弯折
carat 是连续变量里对价格影响最大的一个,值得单独画。散点要设很小的透明度,否则 5 万多个点会糊成一团黑。
ggplot(diamonds_fe, aes(x = carat, y = price)) + geom_point(alpha = 0.04, size = 0.5) + geom_smooth(method = "lm", se = FALSE, colour = "#E45756", linewidth = 1) + geom_smooth(method = "loess", se = FALSE, colour = "#54A24B", linewidth = 1) + scale_x_continuous(breaks = seq(0, 5.5, 0.5)) + theme_minimal()红线是线性拟合,绿线是局部加权回归。两条线在低克拉区几乎重合,超过 1 克拉之后绿线明显上翘,价格增长快于线性。这说明price ~ carat这个模型在低克拉段还行,越往高克拉误差越大。
另外散点里能看到几条竖直的「条带」,那是整数克拉(0.5、1.0、1.5、2.0)的聚集效应——消费者对整数有偏好,同克拉下这些点的价格往往更高。这个现象在建模时没有专门的变量去捕捉,是残差里一块系统性的结构,后面诊断残差图的时候会再遇到。
如果想让线性假设更站得住,两个常规做法:对 price 取对数,或者给 carat 加平方项。我一般先试 log,因为它同时能缓解异方差,一举两得。
4. 建模链路:从简单回归到含交互项的多元模型
4.1 简单线性回归:先把 carat 单独拎出来
不要一上来就堆全变量模型。先用单个变量跑通,确认系数方向和量级合理,再往上加,出问题的时候才知道是哪一步引入的。
fit_carat <- lm(price ~ carat, data = diamonds_fe) summary(fit_carat) confint(fit_carat) # 系数的 95% 置信区间输出里三个数要看:carat 的系数、p 值、R²。系数大致在 7760 上下,意思是克拉数每增加 1,价格预期增加约 7760 美元;p 值远小于 0.001,说明这个关系不是偶然;R² 在 0.85 附近,单靠一个变量就解释了大头,符合「克拉是定价核心」的常识。
截距是负的,大概 -2256。这个数没有现实意义——不存在 0 克拉还倒贴钱的钻石——它只是线性外推的产物。解释模型的时候千万别把这个截距当成「基础价」,这是新手最容易犯的错。
这条模型的价值不在预测精度,而在于给后面的多变量模型立一个基线。加了切工、颜色、净度之后,R² 从 0.85 涨到 0.88,说明这三个定性变量合起来只贡献了 3 个百分点——不是它们不重要,而是它们的作用在简单回归里被 carat 吸收了一部分,得靠交互项才能拆出来。
4.2 定性变量转因子与哑变量全模型
cut、color、clarity 都是分级变量,直接塞进lm()会被当成数值处理,隐含「Fair=1、Good=2……Ideal=5,且相邻档的差距相等」这种强假设,显然不成立。正确做法是让 R 自动生成哑变量。
fit_full <- lm(price ~ carat + cut + color + clarity + depth + table + x + y + z, data = diamonds_fe) summary(fit_full) summary(fit_full)$r.squared # 0.88 上下 car::vif(fit_full) # 多重共线性,重点看 x/y/zR 的处理逻辑是「n 个水平生成 n-1 个哑变量」,被省掉的那个就是基准组。因为我们前面把 levels 定成了从差到好,基准组落在 Fair、J、I1。所以cutIdeal的系数要读作「在克拉数、颜色、净度等其他条件相同的前提下,Ideal 切工比 Fair 切工贵多少美元」。这个「其他条件相同」是关键,也是它能和 3.2 节的分组中位数得出相反结论的原因。
跑完vif()就会看到麻烦:x、y、z 三个变量的 VIF 直接上千。这不奇怪,体积是三者乘出来的,长度和宽度本身相关系数就有 0.95 以上。系数符号会因此变得不可解释,甚至出现 carat 系数为负这种荒唐结果。处理办法是二选一——要么保留 x/y/z 三个原始尺寸,要么保留合成出来的 volume,绝不混着放。我一般留 volume,维度低、解释也方便。
4.3 交互项与逐步回归筛选
全模型有个隐含假设:克拉对价格的边际影响,在所有切工、颜色、净度档次上是同一个常数。这明显不对——一颗 2 克拉的 Ideal 切工钻石和一颗 2 克拉的 Fair 切工钻石,价差应该比 0.3 克拉这一档大得多。交互项就是用来放开这个假设的。
fit_int <- lm(price ~ carat + cut + color + clarity + carat:cut + carat:color + carat:clarity + depth + table + volume, data = diamonds_fe) # 向后逐步回归,按 AIC 删变量 fit_step <- MASS::stepAIC(fit_int, direction = "backward", trace = FALSE) summary(fit_step)三个交互项的含义分别是:carat:cut捕捉「切工溢价随克拉变化」,carat:color捕捉「颜色溢价随克拉变化」,carat:clarity同理。跑出来大概率会看到carat:cutIdeal这类项系数为负,解释是:大克拉时 Ideal 相对 Fair 的溢价被摊薄了,因为大克拉钻石的稀缺性本身就足以支撑高价,切工的边际贡献反而下降。这个结论和 3.2 节的分组现象能对上,互相印证。
stepAIC默认用 AIC(k=2),direction = "backward"表示从全模型开始逐个删。这里trace = FALSE只是关掉每一步的打印,想看删了哪些变量就把它设成 TRUE。如果你的目标是更简约的模型,把准则换成 BIC 更狠——stepAIC(fit_int, direction = "backward", k = log(nrow(diamonds_fe))),惩罚项从 2 变成 log(n),在 5 万多条样本下约为 10.9,会砍掉更多边缘显著的变量。
跑完你会发现所有交互项基本都被留下了,这说明它们确实有信息量。但这不代表模型就没问题了,交互项一多,共线性又会抬头,得回头看 VIF。
4.4 显著性检验与回归诊断
在信模型之前,得先回答两个问题:分组之间的价格差异是不是真的显著?模型的残差符不符合线性回归的前提假设?
第一个问题用非参数检验,因为价格分布右偏严重,不满足方差分析的 normality 假设。
# 分组中位数,先看一眼 diamonds_fe %>% group_by(color) %>% summarise(n = n(), median = median(price)) %>% arrange(desc(median)) # Kruskal-Wallis:不要求正态,比较各组分布位置 kruskal.test(price ~ cut, data = diamonds_fe) kruskal.test(price ~ color, data = diamonds_fe) kruskal.test(price ~ clarity, data = diamonds_fe)三个检验的 p 值都会远小于 0.05,结论是不同切工、颜色、净度之间的价格分布确实存在显著差异。注意 Kruskal-Wallis 检验的是「至少有一组的位置不同」,不是「每一组两两都不同」,想知道具体哪两组有差异还得上事后两两比较。
正态性检验这一步有个硬限制:shapiro.test的样本量必须在 3 到 5000 之间,5 万多条直接报错。
# 随机抽 5000 条做正态性检验 set.seed(42) shapiro.test(sample(diamonds_fe$price, 5000))结果 p 值极小,拒绝正态。但别把这个结论太当回事——样本量到 5000 这个级别,任何轻微偏离都会被判成显著。对大样本来说,QQ 图比 p 值靠谱得多。所以真正的诊断还得看残差图:
par(mfrow = c(2, 2)) plot(fit_full, which = 1:4) par(mfrow = c(1, 1)) # 异方差检验 lmtest::bptest(fit_full)四张图分别看不同问题。Residuals vs Fitted 如果出现喇叭口,就是异方差;QQ 图两端翘起,说明残差厚尾;Scale-Location 进一步确认方差是否恒定;Residuals vs Leverage 看有没有 Cook 距离过大的强影响点,经验阈值是 4/n,5 万条样本下这个线非常低,偶尔几个点超线属于正常。
bptest的 p 值大概率极小,异方差坐实。解决办法是对 price 取对数:
fit_log <- lm(log(price) ~ carat + cut + color + clarity + depth + table + volume, data = diamonds_fe) summary(fit_log)$r.squared # 接近 0.97对数模型的 R² 会明显高于原始尺度,而且系数含义变成「百分比变化」,解释起来更顺:carat 系数乘以 100,大致就是克拉数每增加 1 带来的价格涨幅百分比。如果一定要在原始尺度上做推断,那就别改模型,改用稳健标准误(lmtest::coeftest配合sandwich::vcovHC),系数不变,标准误修正。
5. 排查手册:跑这套分析时最容易踩的五个坑
5.1 数据层的三个坑
坑一:x/y/z 里的 0 值没清,体积特征变成 0。
现象是散点图上出现一条贴在 y 轴的竖线,残差图里对应位置的点杠杆值异常高,volume的分布直方图在 0 处有一根孤零零的柱子。
原因很直接,这部分记录的尺寸字段是 0,x * y * z算出来就是 0。原始数据里这样的行有 20 条左右,is.na()检查不出来,因为它是 0 不是 NA。光看缺失值报告会觉得数据很干净,实际上一脚踩空。
解决就是在构造 volume 之前加一道subset(d2, x > 0 & y > 0 & z > 0),顺手把carat为 0 的也过滤掉。检查习惯是每次造完派生变量,summary()看一眼最小值,非负变量最小值是 0 就该警觉。
坑二:因子水平顺序没定,基准组选错,系数全反着读。
现象是summary(fit_full)里 cut 的系数全是负的,看起来像「切工越好越便宜」,跟 3.2 节的箱线图完全矛盾。
原因是原始 diamonds 里 cut 的 levels 是 Fair 到 Ideal,基准组是 Fair,系数表示「相对 Fair 的差值」,方向本该是正;但如果中途用relevel或factor()重新排过序没注意,基准组可能变成 Ideal,那所有系数就都翻过来了。
解决办法是建模前显式跑一遍levels(diamonds_fe$cut),把它打印出来确认。我现在固定把最差档放第一位,写死在factor(levels = ...)里,不靠默认行为。
坑三:给因子加了ordered = TRUE,系数变成正交多项式。
现象是模型跑完,cut 那一栏的系数名变成了cut.L、cut.Q、cut.C,数值看不懂,也没法解释成「贵多少钱」。
原因是 R 对有序因子默认使用contr.poly对比,生成的是线性、二次、三次正交多项式项,它是为趋势检验设计的,不是为了读单个水平的效应。
解决办法是别加ordered = TRUE,只用factor(levels = ...)定顺序。如果确实需要有序语义(比如做趋势检验),那就得显式指定contrasts = contr.treatment,或者全局设options(contrasts = c("contr.treatment", "contr.poly"))。
5.2 建模层的两个坑
坑四:x/y/z 和 volume 同时进模型,VIF 爆表、系数符号翻转。
现象是car::vif()输出里 x、y、z 的 VIF 上千,carat 的系数变成负数,和 4.1 节简单回归里 7760 的正系数完全对不上。
原因是这四个变量本质上是同一个东西的不同表达,相关系数 0.95 以上,模型无法稳定分配它们各自的贡献,系数估计的方差被放大到失去意义。
解决办法是二选一:要么保留 x/y/z 三个原始尺寸,要么保留合成的 volume,绝不混放。判断标准很简单,跑完模型第一件事就是看 VIF,超过 10 就得处理,超过 100 那模型基本没法用。
坑五:全量数据跑shapiro.test,直接报错中断。
现象是脚本跑到正态性检验那一步抛出sample size must be between 3 and 5000,前面几十分钟的清洗和建模结果全白等。
原因是 Shapiro-Wilk 检验的算法复杂度随样本量增长很快,函数里硬编码了 5000 的上限。
解决办法是抽样,sample(x, 5000)配set.seed()保证可复现。更重要的是心态调整:大样本下正态性检验几乎必然拒绝,与其纠结 p 值,不如直接看 QQ 图和残差-拟合图,那才是判断假设是否成立的实际依据。
6. 二次开发:把一次性脚本改成能复用的小工具
这套脚本跑通之后,最自然的下一步是把它收成一个函数。原始代码是流水账式的,改一个参数要动好几处,封装成函数之后,换数据、换公式、换是否取对数,都只是一个入参的事。
#' 钻石价格回归分析一键流程 #' @param df 原始数据框,需含 carat/price/cut/color/clarity 等列 #' @param formula 回归公式,默认含三个分级变量 #' @param log_y 是否对 price 取对数,默认 TRUE #' @param seed 抽样种子,保证正态性检验可复现 #' @return 列表:清洗后数据、模型对象、R²、VIF 最大值 analyze_diamonds <- function(df, formula = log(price) ~ carat + cut + color + clarity, log_y = TRUE, seed = 42) { set.seed(seed) # 1) 清洗:价格尾部 + 尺寸为正 q <- quantile(df$price, c(0.25, 0.75), na.rm = TRUE) iqr <- IQR(df$price, na.rm = TRUE) clean <- subset(df, price >= q[1] - 3 * iqr & price <= q[3] + 3 * iqr & x > 0 & y > 0 & z > 0) # 2) 特征与因子顺序 clean <- transform(clean, volume = x * y * z, cut = factor(cut, levels = c("Fair","Good","Very Good","Premium","Ideal")), color = factor(color, levels = c("J","I","H","G","F","E","D")), clarity = factor(clarity, levels = c("I1","SI2","SI1","VS2","VS1","VVS2","VVS1","IF")) ) f <- if (log_y) formula else update(formula, price ~ .) fit <- lm(f, data = clean) list( data = clean, model = fit, r2 = summary(fit)$r.squared, max_vif = max(car::vif(fit), na.rm = TRUE) ) }sapply那句清洗逻辑是原有代码的直接搬运,transform比mutate更适合这种不依赖 dplyr 的独立函数。update(formula, price ~ .)是动态换响应变量的写法,避免在函数体里拼字符串。返回值里带上max_vif,是为了让调用方一眼看出这次有没有掉进共线性的坑。
三个模型放在一起对比,差异就清楚了:
| 模型写法 | 响应变量 | 是否含交互项 | R² 量级 | 主要问题 |
|---|---|---|---|---|
price ~ carat | price | 否 | 0.85 | 截距无意义,高克拉段系统性低估 |
price ~ carat + cut + color + clarity + x + y + z | price | 否 | 0.88 | x/y/z 共线性,残差异方差 |
log(price) ~ carat + cut + color + clarity + volume | log(price) | 否 | 0.97 | 系数是百分比语义,绝对值需回代 |
如果还想往前走一步,两条路比较实用。一条是接rmarkdown,把02_eda.R和03_model.R里的图和表直接嵌进参数化报告,换个数据集就能出一份新的 PDF,省掉手工拼图的时间。另一条是接shiny,把predict()包一层,输入克拉、切工、颜色、净度,直接给出价格区间。做这个的时候记得把因子水平的取值范围限制在训练数据覆盖的区间内,否则用户输一个 8 克拉,模型会给出一个荒谬的外推值。
我自己在这类项目上吃过一次亏:早期图省事,把所有变量一股脑塞进全模型,看到 R² 涨了就以为万事大吉,结果交付讨论时被问到「carat 系数为什么是负的」当场卡住。后来才发现是 x/y/z 和 carat 共线,删掉三个尺寸变量,系数立刻回到正轨。从那以后,我每次跑完lm()都强制先走一遍vif()和残差图,R² 再高也不跳过这一步。希望帮到你。
本文还有配套的精品资源,点击获取