R语言钻石价格分析：从数据清洗到多元回归建模

简介：围绕钻石价格预测的R语言数据分析案例文档，面向R语言、回归分析与机器学习入门者，也适合需要完成课程设计或数据分析报告的高校学生。压缩包内共1个docx文件，约2.4MB，内容围绕ggplot2自带的diamonds数据集展开，涵盖数据整体检查、缺失值与重复值排查、箱线图异常值识别与剔除、价格直方图分布观察，以及切工、颜色、净度、克拉重量等变量与价格关系的可视化分析。文档还涉及定性变量因子化、新指标构造、简单线性回归和多元线性回归建模，并尝试用模型解释价格影响因素与预测走势。已有284人学习，可作为R语言数据清洗、探索性分析和回归建模的完整练习素材，帮助读者理解从数据准备到模型解释的基本流程。

1. 从两颗钻的报价差额说起：用 R 语言把 diamonds 数据集跑成一份定价参考

朋友拿了两份报价单过来，一颗 1 克拉 H 色 SI1、一颗 0.9 克拉 D 色 VS2，价差三千多美元，问我哪颗更值。我打开 R，把 ggplot2 自带的 diamonds 数据集跑了一遍钻石价格分析：五万四千多条记录，carat、cut、color、clarity 四个维度铺开，先看分布，再上多元线性回归把每个属性的价格权重估出来。这套流程不复杂，却是数据分析里最典型的回归分析练手场——名义变量和连续变量混在一起，异常值和右偏分布全占齐了。适合刚学完 R 基础语法、想找个真实数据集把清洗、可视化和建模串起来的人，也适合拿它当模板去套自己业务数据的人。后面从工程结构、清洗脚本、建模链路一路讲到怎么改造成能复用的组件。

2. 项目骨架与技术栈：一个能复现的 R 分析工程该有几层

2.1 为什么选 diamonds，而不是自己爬一份数据

diamonds 最大的价值是"干净到可以专注方法"。54000 条记录、10 个字段，三种名义变量（cut、color、clarity）、七个连续变量（carat、depth、table、x、y、z、price），变量类型分布刚好覆盖了回归分析最常遇到的组合：既有需要转哑变量的，也有需要做标准化和变换的。换个业务数据集，光对齐字段口径就能耗掉两天。

代价是它太理想化了，没有时间维度，没有地区维度，也没有真实市场上的促销、品牌溢价。所以这份代码的定位是"方法论模板"，不是"钻石行情预测器"。把这一点想明白，后面读系数、解读图表时才不会被带偏——数据里 Ideal 切工的钻石平均价格并不是最高的，如果不是提前知道这是数据本身的分布问题，很容易写出一个反常识的结论。

我一般会在项目根目录放一个README，第一段就写清数据来源和它的局限，避免半年后自己都忘了这份 R² 到底代表什么。

2.2 技术栈拆解：每个包负责哪一段

整条链路上的包分工明确，不要一上来就library一堆用不上的东西：

机器学习里最朴素的回归分析，其实就靠这几块拼起来。真正容易翻车的地方不是包不够，是因子水平顺序、异常值判定口径这些细节，后面单独用一章讲。

2.3 目录结构和分层思路

一次性脚本谁都会写，但撑不过第二次复用。按"数据 → 清洗 → 探索 → 建模 → 报告"切层，每层一个文件，函数之间只通过数据框传递：

diamond-analysis/ ├── R/ │ ├── 01_load.R # 加载 diamonds，统一字段类型 │ ├── 02_clean.R # 缺失值、重复值、异常值处理 │ ├── 03_eda.R # 描述统计与全部图形 │ ├── 04_model.R # 哑变量、回归、逐步筛选、诊断 │ └── 05_eval.R # 训练/测试拆分、RMSE/MAE ├── output/ │ ├── figures/ # 导出的 png，命名带图号 │ └── models/ # saveRDS 保存的模型对象 ├── report/ │ └── analysis.Rmd # 汇总成可重复渲染的报告 └── run_all.R # 按顺序 source 上面五个文件

分层的好处很实际：调参时只重跑04_model.R，不用每次从读数据开始；图要改配色只动03_eda.R；模型对象用saveRDS落盘，下次预测直接readRDS，省掉重新拟合的几分钟。run_all.R里按顺序source()，保证任何人拿到这份代码，一条命令就能从零跑到底。

2.4 依赖安装与数据加载的第一段代码

# 首次运行时装依赖，已装过的会被 requireNamespace 跳过 pkgs <- c("ggplot2", "dplyr", "tidyr", "mice", "car", "MASS", "broom") for (p in pkgs) { if (!requireNamespace(p, quietly = TRUE)) install.packages(p) } library(ggplot2) library(dplyr) # diamonds 是 ggplot2 的内置数据集，不需要下载 csv data("diamonds", package = "ggplot2") # 立刻做一次体检：维度、字段类型、前六行 dim(diamonds) # 53940 x 10 str(diamonds) # 确认 cut/color/clarity 是 chr 或 Ord.factor head(diamonds, 6)

data()里的package参数别省，工作区里如果存在同名的diamonds变量，不加参数会静默加载错的那个，这类问题排查起来极其费时。str()一定要看，diamonds 里的 cut、color、clarity 在不同 R 版本下可能是chr、Ord.factor或ordered，后续lm()对有序因子的处理方式和普通因子完全不同，这是建模阶段第一个要盯的点。

3. 数据清洗与分布画像：缺失值、重复值、异常值三道闸

3.1 缺失值排查：is.na全 0 不等于没有缺失

# 逐列统计 NA 数量 colSums(is.na(diamonds)) # 用 mice 的 md.pattern 看缺失模式矩阵 library(mice) md.pattern(diamonds, rotate.names = TRUE)

colSums(is.na())返回全 0，md.pattern画出来的矩阵里也没有 0 值行，说明这批数据在 NA 意义上的确完整。但这里有个陷阱：x、y、z三个尺寸字段里存在取值为 0 的记录，is.na()是查不出来的。一颗钻石长度为 0 毫米、宽度 58.9 毫米，物理上不成立，本质是录入错误，只是被填成了 0 而不是留空。

我一般把这类值叫"隐性缺失"，处理方式和 NA 一样——直接剔除，但不能用na.omit()一刀切，得显式写条件，方便日后回溯剔了多少条。

3.2 重复值与异常值：箱线图先看，再定剔除规则

# 完全重复行 sum(duplicated(diamonds)) # 结果是 0 行 # 三围字段的合理性过滤，同时干掉 0 尺寸和明显超限的记录 diamonds_clean <- diamonds %>% filter(y > 0, y < 20, z > 0, z < 10, x > 0) # price 用 ±3 标准差之外视为异常 mu <- mean(diamonds_clean$price) sd_ <- sd(diamonds_clean$price) diamonds_clean <- diamonds_clean %>% filter(between(price, mu - 3 * sd_, mu + 3 * sd_)) # 剔除前后对比数量 nrow(diamonds) - nrow(diamonds_clean)

filter(y > 0, y < 20)这两个阈值不是随手拍的，y 字段官方范围标注到 58.9 毫米，超出正常裸钻尺寸的记录基本都是数据录入问题；20 毫米留了余量，能保住大克拉钻又不至于误伤。price 用 ±3σ 是因为价格近似右偏，均值加三倍标准差大约落在两万美元左右，恰好把极少数天价记录排除，同时保留 1 万到 1.8 万这一段真实存在的高价钻。

剔除规则一定要落到代码里写着，别在控制台手动删几行就完事——三个月后别人问"你这 53940 怎么变成 53000 了"，你至少还能答上来。

3.3 描述性统计：先把价格骨架摸清楚

summary(diamonds_clean$price) # 最小 326，最大 18823，均值 3933 附近 quantile(diamonds_clean$price, probs = c(0.25, 0.5, 0.75, 0.9, 0.99)) # 分位数 + 偏度峰度，psych::describe 一步到位 psych::describe(diamonds_clean$price)

看这几个数：中位数明显低于均值，说明右尾被高价钻石拉长；99% 分位在 1.8 万左右，说明极端值很少但仍然存在。这个"均值 > 中位数"的形态是所有奢侈品数据的共性，也是后面必须考虑对 price 取对数的直接原因。

描述统计只做这一步还不够，还得看分组后各组的价格中枢有没有差别，否则你不知道到底该把哪些变量放进模型，也不知道有没有变量纯属噪音。

3.4 直方图与核密度：右偏、峰值、断层一次看全

library(ggplot2) ggplot(diamonds_clean, aes(x = price)) + geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "#4C72B0", colour = "white") + geom_density(colour = "#C44E52", linewidth = 1) + stat_function(fun = dnorm, args = list(mean = mean(diamonds_clean$price), sd = sd(diamonds_clean$price)), colour = "#55A868", linetype = "dashed") + labs(x = "价格（美元）", y = "密度", title = "钻石价格分布：直方图 + 核密度 + 正态参照") + theme_minimal(base_family = "")

after_stat(density)是新一代写法，老代码里的..density..在新版 ggplot2 会报警告，虽然还能跑，但建议统一换掉，不然控制台一堆 warning 会淹没真正的错误信息。stat_function叠的那条虚线是同样均值标准差的正态曲线，两条线一比就能看出偏差：低于 8000 美元这段，核密度曲线明显比正态曲线高而窄；8000 以上两线几乎重合。这个形态说明简单线性回归的残差不会服从正态，后面做诊断时必然会亮红灯。

bins = 60是试出来的——30 太粗看不清断层，100 以上每个柱子只有个位数样本，反而全是噪音。这类参数没有标准答案，出图看一眼再定。

3.5 分组可视化：切工、颜色、净度谁在推价格

# 小提琴图：既有分布形态又有中位数，比箱线图信息量大 ggplot(diamonds_clean, aes(x = cut, y = price, fill = cut)) + geom_violin(trim = FALSE, alpha = 0.7) + geom_boxplot(width = 0.12, outlier.shape = NA, fill = "white") + scale_y_log10() + labs(x = "切工等级", y = "价格（对数刻度）") + theme_minimal() + theme(legend.position = "none")

价格用scale_y_log10()转成对数刻度，不转的话所有分布会被压成贴着 x 轴的一条线，中位数差异完全看不出来。图上能直接读出的结论有两个：一是 Ideal 切工那一组的下四分位明显偏低，说明存在大量小克拉的高切工钻石，切工好不代表贵；二是各组的分布宽度差异远大于中位数差异，意味着 cut 单独解释价格的能力很弱，必须在模型里和 carat 交互才有意义。

颜色和净度同理，把aes(x = color)和aes(x = clarity)换上去各跑一次就行，不用改其他参数。真正值得注意的是 color 从 J 到 D 的均价是一条单调上升的曲线，而 clarity 在 VVS1/VVS2 之后出现回落——高档位的净度差异对价格的边际影响开始变缓，这条结论在后面解读回归系数时能对上。

4. 回归建模全链路：哑变量、交互项与逐步回归怎么搭

4.1 因子编码：有序因子会悄悄改掉你的系数含义

diamonds_model <- diamonds_clean %>% mutate( # 保留原始等级顺序用于画图 cut_lv = factor(cut, levels = c("Fair","Good","Very Good","Premium","Ideal")), color_lv = factor(color, levels = c("J","I","H","G","F","E","D")), clarity_lv = factor(clarity, levels = c("I1","SI2","SI1","VS2","VS1","VVS2","VVS1","IF")), # 建模用的无序因子，对照水平显式指定 cut_f = relevel(factor(cut), ref = "Fair"), color_f = relevel(factor(color), ref = "J"), clarity_f = relevel(factor(clarity), ref = "I1"), # 派生变量：体积、每克拉单价 volume = x * y * z, ppc = price / carat )

关键区别：ordered = TRUE的因子进lm()会走多项式对比（.L、.Q、.C），输出的是线性趋势、二次趋势系数，不是"相对基准组高多少美元"。很多人看系数表里出现cut.L、cut.Q一头雾水，就是这里来的。建模时统一用无序因子加relevel指定基准，基准选最差那一档，系数解读起来最符合直觉——"比 Fair 贵多少"。

volume = x*y*z和ppc = price/carat这两个派生变量不是凑数：volume把三个高度相关的尺寸字段压缩成一个，能显著降低共线性；ppc是业务上最常用的性价比指标，做残差分析时拿它当 y 轴，异常样本一眼就能挑出来。

4.2 一元回归：先看 carat 到底能解释多少

fit_carat <- lm(price ~ carat, data = diamonds_model) summary(fit_carat)

系数约 7760，p 值远小于 0.001，直白地说：每多 1 克拉，价格预期涨 7760 美元。但 R² 只有 0.85 出头，残差图里明显是喇叭口——小克拉时残差集中在 0 附近，大克拉时残差散开到正负几千美元。这说明关系是非线性的，carat 越高，价格离散度越大。

这一步的意义不是拿这个模型去预测，而是给后面的多元回归定一个基准线。如果加了七八个变量，R² 才从 0.85 提到 0.87，那就得想想这些变量到底贡献了什么。

4.3 多元回归与交互项：0.8798 的 R² 卡在哪

# 无交互项主效应模型 fit_main <- lm(price ~ carat + cut_f + color_f + clarity_f + depth + table + volume, data = diamonds_model) summary(fit_main)$r.squared # 约 0.8798 # 加三个业务上说得通的交互项 fit_inter <- lm(price ~ carat * cut_f + carat * color_f + carat * clarity_f + depth + table + volume, data = diamonds_model) # 向后逐步回归，AIC 准则筛变量 fit_step <- step(fit_inter, direction = "backward", k = 2, trace = 0) summary(fit_step)$r.squared

主效应模型 R² 停在 0.88 附近，再往上堆变量提升有限，因为真正的信息藏在"克拉数和品质的交互"里：同样是加 1 克拉，D 色和 J 色的价格增幅完全不同，这层关系用加法模型表达不出来，必须写成carat * color_f。

step()的k = 2对应 AIC，改成k = log(nrow(diamonds_model))就是 BIC，惩罚更重、留下的变量更少。trace = 0关掉每步的打印，不然控制台会刷出几百行。向后逐步回归从全模型出发逐步删变量，比向前法更稳，代价是初始模型变量多、计算慢，五万行的数据大概几十秒能出结果。

4.4 显著性检验与回归诊断：模型能不能信

# 非参数检验：不假设正态，比较不同切工组的price分布 kruskal.test(price ~ cut, data = diamonds_model) # 颜色组间两两比较，BH 校正控制假阳性 pairwise.wilcox.test(diamonds_model$price, diamonds_model$color, p.adjust.method = "BH") # 共线性诊断，VIF 超过 10 的变量必须处理 car::vif(fit_step) # 残差四图：线性、正态、同方差、杠杆值 par(mfrow = c(2, 2)) plot(fit_step)

kruskal.test的 p 值极小，说明不同切工组的价格分布确实存在差异，但注意它只能告诉你"有差异"，不能告诉你"差多少"，具体幅度还得看回归系数。pairwise.wilcox.test做多重比较必须加p.adjust.method，不加的话 21 组两两比较里必然冒出假阳性。

vif()是最容易被跳过的一步。carat和volume相关系数接近 0.95，两个都放进去 VIF 会飙到 15 以上，系数符号甚至会翻转。处理办法通常二选一：删掉volume只留carat，或者改用carat加depth/table的组合。

plot(fit_step)的四张图必须逐张看。残差 vs 拟合值如果呈漏斗形，说明异方差；QQ 图两端翘起，说明残差重尾，直接原因就是 price 没取对数。

4.5 训练测试拆分与误差评估

set.seed(2024) idx <- sample(nrow(diamonds_model), 0.7 * nrow(diamonds_model)) train <- diamonds_model[idx, ] test <- diamonds_model[-idx, ] # 对数变换缓解右偏，模型更贴合残差假设 fit_log <- lm(log(price) ~ carat + cut_f + color_f + clarity_f + volume, data = train) pred_log <- predict(fit_log, newdata = test) pred <- exp(pred_log) # 还原到美元尺度 rmse <- sqrt(mean((pred - test$price)^2)) mae <- mean(abs(pred - test$price)) c(RMSE = rmse, MAE = mae) # 保存模型，供后续预测复用 saveRDS(fit_log, "output/models/fit_log.rds")

set.seed固定随机种子，保证每次拆分一致，否则两次跑出来的 RMSE 能差几百，调参就变成玄学。取对数后模型 R² 通常能上到 0.93 以上，不是因为模型变强了，而是因为 log 尺度下的相对误差被压平了，这一点在汇报时必须说清楚，别把指标提升当成算法改进。

exp(pred_log)是朴素还原，严格来说对数模型的预测值需要加一个 smearing 修正因子（残差均值的指数），量级上和朴素还原差不了几十美元，做业务估算可以忽略，做精确报价就得补上。

5. 避坑与排查：从因子顺序到 log 变换的五条血泪记录

5.1 缺失值全是 0，模型却报错"NA/NaN/Inf in 'y'"

现象：colSums(is.na(diamonds))返回全 0，但跑lm()时报NA/NaN/Inf in 'y'。

原因：x或y字段存在取值为 0 的记录，is.na()判定不出来，但 0 宽度的钻石会让volume = x*y*z变成 0，取对数后就是-Inf。

解决：清洗阶段显式加一行filter(x > 0, y > 0, z > 0)，并且在建模前统一跑一次sum(!is.finite(diamonds_model$volume))作为守门检查。从那以后我每次派生新变量都会顺手查一遍is.finite。

5.2 系数全是 NA，或者出现cut.L、cut.Q

现象：summary()里切工相关的行显示NA，或者变量名变成cut.L、cut.Q、cut.C。

原因：cut 字段被读成了ordered因子，lm()自动套用多项式对比；另一种情况是某个因子的某个水平样本量为 0，形成完全共线。

解决：建模前统一factor(cut)转成无序，再用relevel指定基准组；用table(diamonds$cut)确认每个水平都有足够样本。这段判断现在被我固化成一个小函数，每次建模前先跑。

5.3 逐步回归筛完，R² 反而比主效应模型低

现象：step()砍掉一堆变量后，调整 R² 下滑，预测误差也变大。

原因：k参数设成了 2（AIC），对五万行的大样本来说惩罚太轻，容易把噪声变量留在模型里；另一个常见原因是交互项写得太宽，(.)^2一次性生成上百个交互，逐步筛选根本扛不住。

解决：样本量大时改用 BIC，k = log(nrow(train))；交互项只挑业务上说得通的三个（carat 和三个品质变量），其余不生成。变量多不等于模型好，这是我在这个数据集上吃过的最实在的教训。

5.4 Ideal 切工均价反而低，被质疑结论写反了

现象：分组统计出来 Ideal 切工的平均价格低于 Premium，和"切工越好越贵"的常识冲突，评审直接质疑分析有误。

原因：这是典型的混杂效应。carat 和 cut 之间存在相关性，Ideal 组里集中了大量小克拉钻石，均价被拉低；单独看 cut 的边际效应是负的，控制 carat 之后系数符号会翻过来。

解决：所有关于单变量的结论必须先做分组控制，或者直接读多元回归里的偏效应系数。图表标题里写清"未控制克拉重量"，能省掉很多解释成本。

5.5 残差 QQ 图两端翘，R² 怎么调都上不去

现象：模型 R² 卡在 0.88，QQ 图两端明显偏离对角线，预测高价钻石时系统性偏低。

原因：price 是右偏分布，线性模型假设残差正态，右尾的长尾数据把拟合线往下拽。

解决：对 price 取自然对数后建模，预测时再用exp()还原。副作用是模型变成"预测价格的对数"，误差解释要改成相对误差百分比，汇报口径得同步调整，别拿 RMSE 直接和原尺度模型比。

6. 二次开发改造：把一次性脚本升级成可复用的定价组件

把这份分析从"交作业"改造成"能给别人用"，我走过几条路，按投入产出比排一下：

最值得先做的是函数化。把02_clean.R里的逻辑抽成一个clean_diamonds(df, sigma = 3)，阈值当参数传进去，下次换成自己公司的商品数据，只要字段名对齐就能直接跑：

clean_diamonds <- function(df, sigma = 3, y_max = 20, z_max = 10) { stopifnot(all(c("price", "x", "y", "z") %in% names(df))) df <- df %>% filter(x > 0, y > 0, y < y_max, z > 0, z < z_max) %>% filter(between(price, mean(price) - sigma * sd(price), mean(price) + sigma * sd(price))) attr(df, "removed") <- nrow(df) # 留个痕迹，方便回溯 df } predict_price <- function(model, newdata) { exp(predict(model, newdata = newdata)) # 与训练时的 log 变换配套 }

stopifnot那行是防御性的，字段缺失时直接报错，而不是让filter静默返回零行——零行数据流到后面，报的错会离真正的原因十万八千里。predict_price里把exp()绑死在预测函数内，这样调用方永远不用记住"这个模型是对数尺度的"，接口干净，出错概率也低。

第二个值得做的是交叉验证。现在用的单次 7:3 拆分，RMSE 会随随机种子浮动，做模型对比时根本分不清差异来自模型还是来自拆分。换成rsample::vfold_cv(v = 5)，指标立刻稳定下来，随机森林和线性模型的差距才有可比性。我试过一次，线性模型 RMSE 大约 1100 美元，随机森林能压到 900 上下，但可解释性掉了大半——给消费者做报价参考，能说清"贵在克拉还是贵在净度"比省两百美元误差更重要。

最后一个小心思：把显著性检验的结果也写进报告。很多人只贴 R² 和系数表，但业务方最常问的是"颜色这一档到底有没有影响到价格"，kruskal.test的 p 值和两两比较的结果正好回答这个问题，成本很低，说服力提升明显。

这个数据集我前后跑过七八遍，每次重跑都会发现前一版漏掉的东西。最惨的一次是直接拿原始 diamonds 建模，异常值和 0 尺寸记录一个没清理，交互项跑出来系数符号全反，图表交上去被当场问住，只能回去重来。从那以后我每次接手一份新数据，都强制先过三关：summary()看量纲和极值、箱线图看离群点、残差图看模型假设，三张图都没问题才允许动模型。希望帮到你。

本文还有配套的精品资源，点击获取