R语言钻石价格回归分析：从数据清洗到模型诊断的完整实战

简介：围绕钻石价格影响因素展开的R语言数据分析实践文档，面向具备一定R基础、希望将回归分析与可视化结合成一个完整项目的数据分析学习者和统计课程使用者。内容基于ggplot2自带的diamonds数据集，覆盖数据整体信息查看、缺失值与重复值检查、箱线图异常值识别与剔除、价格分布直方图与核密度分析，以及切割、颜色、净度、克拉重量等因素与价格关系的分组可视化，最终进行多元线性回归建模，并给出价格预测与走势判断。资源包共1个docx文件，约2.4MB，以图文排版的实验记录形式呈现，变量说明、代码结果与分析结论相互对应，便于按小节复现，也便于整理成课程报告或项目笔记。已有284人学习下载，适合作为回归分析从入门到综合应用的参照范例。

1. 从一个报价问起：为什么钻石定价要用回归而不是经验

钻石价格分析这个题目，难点从来不在 R 语言写不出回归代码，而在于你凭什么说某个报价“合理”。我见过不少人，加载完 ggplot2 自带的 diamonds 数据集就直接lm(price ~ .)，跑出个 R² 约 0.88 就收工，答辩时被追问三句——异常值凭什么删、交互项为什么加、三克拉以上敢不敢预测——就露了底。这套方案走另一条路：先把约 54000 条记录拆成变量族并定下清洗口径，再做分布与分组探索找出价格偏态和断层，然后用一元回归热身、哑变量建模、交互项配合向后逐步回归来完善，最后给出一套能自证的验收检查点。它适合要交数据分析类毕设课设、需要完整流程和可复现代码的读者，也适合手上有报价数据、需要一套能讲清楚的定价口径的从业者。

2. 数据底座：diamonds 十个变量的分层与三步清洗口径

任何一份价格分析报告，可信度都是从数据底座开始往下掉的。diamonds 这个数据集好处是干净、公开、量级够，坏处也恰恰在这里——太干净会让人放松警惕，把“没有缺失值”当成“不用清洗”。这一章要定下三件事：哪些变量进模型、按什么口径剔除记录、定类和派生变量怎么处理。这三件事定完，后面所有图、所有系数才有可比性，否则换一个人重跑一遍结果就对不上，这种翻车在评审现场最难解释。

2.1 变量分层：三个定类、七个连续，先决定谁进模型

先把十个变量摊开看。因变量只有一个——price，单位美元；自变量按性质分成三族：连续型尺寸（carat、x、y、z）、比率型工艺参数（depth、table）、定类型分级（cut、color、clarity）。

表里最后一行要特别留意：x、y、z 三个尺寸和 carat 高度相关，甚至可以说 carat 就是它们的综合表达。如果目标是“解释哪些因素影响价格”，把它们和 carat 一起塞进模型会让系数互相稀释；如果目标是“预测价格”，加进去能明显压低残差。我一般会准备两个模型：解释版只保留 carat + 三个定类 + depth + table，预测版才把尺寸放进去，答辩时分别讲，逻辑就顺了。

另外，y 的最大值 58.9 毫米，和 x 的 10.74 毫米完全不在一个量级——这种记录在钻石尺寸里不成立，属于典型录入错误，要在清洗阶段处理掉，别指望模型自己消化。

2.2 缺失、重复、异常值：三步清洗的判断口径

第一步查缺失。diamonds 本身没有 NA，但流程要写出来，这是给别人看你的分析是不是走完整套的关键。

library(ggplot2) library(mice) data(diamonds) # 整体结构与前六行 str(diamonds) head(diamonds, 6) # 缺失值模式图：矩阵中每一列代表一个变量，红色块表示缺失位置 md.pattern(diamonds) sum(is.na(diamonds)) # 应为 0

第二步查重复。duplicated()统计的是整行完全一致的记录，钻石数据里存在规格完全相同的样本，这在业务上未必是错误（同批次同参数的两颗钻石很正常），所以判据不是一律删，而是看重复比例——超过 5% 才需要回头核对采集环节。

sum(duplicated(diamonds)) # 重复行数量 diamonds <- dplyr::distinct(diamonds) # 保险起见去重

第三步才是重头戏：异常值。原文用的是三个标准差口径，箱线图上价格那一列的中位数贴近下四分位，说明存在明显右偏，超出一侧 3σ 的点不少。这里给一段可以直接复用的代码，把连续变量统一打标：

num_vars <- c("carat", "depth", "table", "price", "x", "y", "z") # 每个连续变量标准化后取绝对值，超过 3 记为一处越界 flag_mat <- sapply(diamonds[num_vars], function(v) abs(scale(v)[, 1]) > 3) # 一颗钻石只要任一变量越界就剔除；想保守一点可改成 rowSums(flag_mat) >= 2 diamonds_clean <- diamonds[rowSums(flag_mat) == 0, ] # 清洗前后行数对比，是验收时要交代的第一个数字 cat("清洗前:", nrow(diamonds), " 清洗后:", nrow(diamonds_clean), "\n")

参数上有两个选择要想清楚。一是阈值，3σ 适合分布接近对称的变量，对 price 这种右偏变量会误删不少高价真样本，替代方案是分位数截尾（如保留 1% ~ 99%）。二是“任一越界即删”还是“至少两处越界才删”，前者干净但样本损失大，后者保留信息但脏点会稀释系数。我一般先用“任一越界即删”跑一版，如果剔除比例超过 8%，就退回到分位数截尾，理由是“样本大量丢失带来的偏差，比几个极端值更伤模型”。

同时补一条硬规则：x、y、z 任一为 0 的记录直接删。物理尺寸为零的钻石不存在，这是量测缺失被填成 0 的典型情况，留着会严重压低回归线。

2.3 因子化与派生变量：让 depth、table 能进模型

cut、color、clarity 三个变量本质是等级，直接当字符型向量进线性回归，R 会把它们当成一堆互不相干的水平，模型就没法利用“等级递进”这层信息。因子化并且显式声明水平顺序，是把这层信息交代清楚的唯一方式：

# 有序因子：levels 由差到好，ordered = TRUE 让后续对比有明确方向 diamonds_clean$cut <- factor(diamonds_clean$cut, levels = c("Fair", "Good", "Very Good", "Premium", "Ideal"), ordered = TRUE) diamonds_clean$color <- factor(diamonds_clean$color, levels = c("J", "I", "H", "G", "F", "E", "D"), ordered = TRUE) diamonds_clean$clarity <- factor(diamonds_clean$clarity, levels = c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF"), ordered = TRUE) # 派生变量：体积用三边乘积，避免和 carat 完全同义，留给模型做补充解释 diamonds_clean$volume <- with(diamonds_clean, x * y * z) # 校验 depth：官方定义是 2*z/(x+y)*100，算出来和原列做差应接近 0 diamonds_clean$depth_check <- with(diamonds_clean, 2 * z / (x + y) * 100) summary(diamonds_clean$depth - diamonds_clean$depth_check)

这段里有个容易被忽略的点：ordered = TRUE只影响对比编码的方式，lm()默认仍然按处理对比（treatment contrast）生成 k−1 个哑变量，基准水平取第一个 level，也就是 Fair、J、I1。想让基准变成“最理想档”便于解读，用relevel()换一下就行；想让系数表示“相邻等级之间的增量”，改用contr.sum或多项式对比。这个选择不改变模型拟合优度，但直接决定你能不能把系数讲成人话——答辩时被问“你这个 Ideal 系数 −800 是什么意思”，答案完全取决于基准是谁。

还有一点：depth 和 table 都是百分比，depth = 2z/(x+y)这条恒等式意味着它和 x、y、z 存在结构性关联。上面那行depth_check的差值如果集中在 0 附近，说明数据自洽；如果出现成片的偏差，说明有记录的量测方式不一致，这类样本在解释型模型里要单独标注。

3. 探索性分析：价格分布形状与三组分类变量的读图顺序

清洗完直接进模型是新手最常见的跳步。价格这种右偏严重的变量，不做分布诊断就上线性回归，残差图一定难看；不做分组对比就直接把 cut、color、clarity 丢进去，系数符号反常时你连“是数据问题还是模型问题”都判断不了。这一章的三组图，作用不是好看，是给第 4 章的每一个建模决策提前找到依据。

3.1 价格分布：右偏、千元峰值和三处断层

先看整体形状。直方图叠加核密度曲线，再加一条均值线，一张图能读出三件事：偏态方向、峰值位置、长尾厚度。

library(ggplot2) p_summary <- summary(diamonds_clean$price) print(p_summary) # 最小值、四分位、均值、最大值 cat("中位数:", median(diamonds_clean$price), "\n") 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) + geom_vline(xintercept = mean(diamonds_clean$price), linetype = "dashed", colour = "grey30") + labs(x = "价格（美元）", y = "密度", title = "钻石价格分布")

清洗后的典型读数是：最低 326 美元、最高约 11886 美元、均值在 3160 美元上下、中位数明显低于均值。均值大于中位数就是右偏的直接证据，长尾全部落在高价侧。密度曲线在 1000 美元附近出现峰值，说明这个价位是市场接受度最集中的区间；2000、4000、8000 三个位置出现明显断层，通常和克拉整数档（0.5、1.0、2.0 克拉）的定价台阶有关。

这张图给建模的直接结论是：price 不适合直接做线性回归的因变量。要么对 price 取对数，要么在解释时明确“系数是平均边际效应，不适用于高价段”。我在方案里一般两个都做：主模型用原始 price 保证系数可读（美元/克拉），诊断章节再用 log(price) 验证一遍结论方向是否一致。

3.2 分组对比：切工、颜色、净度三张图各自说明什么

单变量看完，接着按三个定类变量分组。这里统一用箱线图或小提琴图，纵轴取 log10 是为了压住长尾、让组间中位数的差异看得清。

# 切工与价格：小提琴图看分布形状，箱线看中位数 ggplot(diamonds_clean, aes(x = cut, y = price, fill = cut)) + geom_violin(alpha = 0.6, trim = FALSE) + geom_boxplot(width = 0.12, outlier.size = 0.4) + scale_y_log10() + labs(x = "切工", y = "价格（美元，log10）") + theme(legend.position = "none") # 颜色与净度：同样结构，换 x 轴变量即可 ggplot(diamonds_clean, aes(x = color, y = price, fill = color)) + geom_boxplot(outlier.size = 0.4) + scale_y_log10() + labs(x = "颜色（J 差 → D 好）", y = "价格（美元，log10）") + theme(legend.position = "none")

三张图各自给出一条不能省的结论。切工那张，Ideal 组的中位数并不是最高的，甚至有一批 Ideal 钻石价格偏低——这不是数据错误，而是小克拉钻石也能做到理想切工，切工等级和克拉数之间存在结构性交叉，正是这个现象让第 4 章的交互项有了依据。颜色那张，从 J 到 D 中位数基本单调上升，方向清晰，说明这个变量进模型后系数应该稳定为正。净度那张最反直觉：VVS1、VVS2 的价格中位数最高，但 I1 并不是最低，SI2 一带反而低——净度到达优质区间之后，区分价格的维度就更多了，单一等级解释不了。

顺带补一张计数图，确认各等级的样本量，避免某个等级只有几十条却拿来下结论：

table(diamonds_clean$cut) table(diamonds_clean$color) table(diamonds_clean$clarity)

计数结果通常显示 Ideal 与 Premium 合计占了一半以上，而 Fair 只有一千多条。样本量悬殊意味着 Fair 组的均值方差都大，做组间检验时要留意，这一点放到 3.4 之前先说清楚。

3.3 尺寸与价格：为什么 carat 一定要做变换或加交互

散点图是最直白的一步。横轴 carat、纵轴 price，直接画出来是一团向右上方发散的扇形——均值随克拉上升，但离散程度上升得更快，这就是异方差的视觉证据。

ggplot(diamonds_clean, aes(x = carat, y = price)) + geom_point(alpha = 0.12, size = 0.6) + geom_smooth(method = "lm", formula = y ~ x, colour = "#C44E52") + labs(x = "克拉重量", y = "价格（美元）") # 双对数坐标：如果近似一条直线，说明幂律关系成立 ggplot(diamonds_clean, aes(x = carat, y = price)) + geom_point(alpha = 0.12, size = 0.6) + geom_smooth(method = "lm", formula = y ~ x, colour = "#4C72B0") + scale_x_log10() + scale_y_log10() + labs(x = "克拉（log10）", y = "价格（美元，log10）")

双对数坐标下散点基本贴合一条直线，斜率通常在 1.7 附近，换算过来就是克拉翻一倍、价格涨约 3.2 倍。这个数字比“每克拉涨多少钱”更贴近真实定价逻辑，也是后面加交互项的理论支点：既然 carat 的边际效应不是常数，那么 carat 与切工、颜色、净度的乘积项就应该显著。

到这里探索部分可以收口了，三条待验证的假设已经列清楚：price 右偏需要变换或加权、三个定类变量的效应方向不完全单调、carat 的边际效应随等级变化。第 4 章的模型结构，就是围绕这三条搭的。

4. 建模路线：从一元回归到交互项加逐步回归的四个台阶

建模最忌讳一步到位。lm(price ~ .)这种写法跑得出来，但每一步在干什么、为什么这么设，自己都说不清。这一章按四个台阶往上走：一元回归建立基线并暴露问题，哑变量模型解决“定类变量怎么进”，交互项模型回答“效应是否随等级变化”，逐步回归负责在解释力和复杂度之间找平衡。四个模型的输出放在一起对比，才是能拿去讲的东西。

4.1 一元回归热身：carat 系数约 7760 美元/克拉怎么读

只用克拉数预测价格，是最小可用的基线模型。它的价值不在预测精度，而在于把“单因素解释力”这个参照点钉住。

fit1 <- lm(price ~ carat, data = diamonds_clean) summary(fit1) confint(fit1, level = 0.95) # 系数的 95% 置信区间

输出里三件事要读。第一，carat 系数约 7760，p 值远小于 0.001，说明克拉每增加 1 单位、价格平均增加约 7760 美元，这个方向没有争议。第二，截距是显著的负值，这意味着 carat = 0 时模型预测出负价格——纯粹的无意义外推，提醒你线性模型只在观测区间内有效。第三，R² 通常在 0.85 上下，看着很高，但残差图会暴露问题：

plot(fit1, which = 1) # 残差 vs 拟合值

残差随拟合值增大而扩散成喇叭口，这是异方差。异方差本身不会让系数估计有偏，但会让标准误失真，p 值和置信区间都不可靠。解法有三条：对 price 取对数、用lm(log(price) ~ carat)；用加权最小二乘，权重取 1/carat²；或者老老实实承认区间预测在高价段不可靠。我一般主模型保留原始尺度便于解释，同时在诊断章节跑一遍对数模型验证结论方向一致。

4.2 哑变量与无交互项模型：R² = 0.8798 卡在哪

把三个定类变量和工艺参数加进来，得到第一个多元模型。定性变量转哑变量这一步不用手工做，lm()会自动按因子水平生成 k−1 列，基准水平是各因子的第一个 level。

fit2 <- lm(price ~ carat + cut + color + clarity + depth + table, data = diamonds_clean) summary(fit2) anova(fit1, fit2) # 嵌套模型比较：加了这些变量到底有没有用

读数要点：R² 落在 0.88 附近，比一元模型提升有限，说明 carat 一个人扛下了大部分解释力；anova()的 F 检验显著，说明三个定类变量整体上是有效的，不能删。逐个看系数，color 从 J 到 D 方向单调（越接近 D 价格越高），clarity 大体同向，cut 的系数则可能符号摇摆——这正是 3.2 里“Ideal 组中位数不高”的模型侧印证。

depth 这个变量要单独说。它的系数经常不显著，甚至出现负号，看着像“深度越大越便宜”，其实是因为 depth 与 carat、x、y、z 存在结构性关联，被其他变量吸收了信息。处理方式不是硬删，而是在报告里单独说明“depth 的最优区间在 61% ~ 62.5%，偏离两侧都会折价”，把它从线性项改成分段或加二次项I(depth^2)，往往能让符号恢复正常。

4.3 交互项与向后逐步回归：AIC 准则下怎么取舍

第四个台阶是给 carat 装上“随等级变化的斜率”。理论上讲得通：一克拉的 Ideal 切工钻石和两克拉的 Fair 切工钻石，克拉溢价的幅度不一样，所以 carat 与 cut、color、clarity 的乘积项应该有解释力。

# 用 (A + B + C)^2 展开，等价于主效应 + 两两交互，避免手写一长串 fit_full <- lm(price ~ (carat + cut + color + clarity)^2 + depth + table, data = diamonds_clean) # 向后逐步回归：k = 2 是 AIC 准则；改成 k = log(nrow(diamonds_clean)) 即 BIC，更保守 fit_step <- step(fit_full, direction = "backward", k = 2) summary(fit_step) AIC(fit1, fit2, fit_step)

step()的逻辑是每一步尝试删掉一个变量，看 AIC 是否下降，下降就删，直到删不动为止。AIC 惩罚项较轻（每多一个参数罚 2），倾向于保留更多项；BIC 惩罚项随样本量增长（log(n)），在五万多条样本下会明显更严格，模型更精简。实际跑下来，交互项模型通常会保留 carat 与三个定类变量的交互，R² 提升到 0.98 附近，看起来很好看，但这一步有个必须提前想清楚的问题：参数数量翻了几倍，每个系数的业务含义还讲得出来吗？

我的处理原则是双模型并存。解释版用 fit2，系数少、方向清楚，报告正文以它为主；预测版用 fit_step，只出现在精度评估和误差分析里，明确标注“用于预测，不用于归因”。这条界线一旦模糊，评审现场很容易被问住。

4.4 模型诊断：残差图、共线性与影响点

模型跑完不看诊断，等于开车不看仪表。四项检查按顺序走：

par(mfrow = c(2, 2)) plot(fit_step) # 残差vs拟合值、QQ图、标准化残差、杠杆值 par(mfrow = c(1, 1)) car::vif(fit2) # 方差膨胀因子，> 10 提示严重共线 cooksd <- cooks.distance(fit_step) thr <- 4 / nrow(diamonds_clean) sum(cooksd > thr) # 强影响点数量，超过 1% 就要回头查

残差 vs 拟合值看异方差是否缓解，QQ 图两端翘说明尾部不服从正态，杠杆值图能定位那几颗把回归线拽偏的样本。VIF 这一项在解释版模型上要重点看，carat 和 x、y、z 一旦同时入模，VIF 很容易冲到 10 以上，此时系数的正负号就不再可靠。影响点用 Cook 距离筛，阈值 4/n 只是起步线，真正要做的不是删掉这些点，而是把它们单独列出来，看是不是高克拉少见的规格——如果是，就在报告里说明“模型对超大克拉样本敏感”，这比假装没问题要专业得多。

5. 避坑与排查：这份分析里五个最容易翻车的环节

5.1 异常值一刀切：样本掉了两成，系数全变

现象：清洗后样本从 53940 掉到 40000 出头，重新拟合时 carat 系数从 7760 变成 9200，颜色系数方向也动了。

原因：abs(scale(v)) > 3对右偏的 price 极不友好，高价真样本被大批误删，剩下的样本高价侧被截断，回归线自然变陡。另外“任一变量越界即删”的规则在十个变量上叠加，误删率会放大好几倍。

解决：先统计每个变量的剔除贡献，找出主要来源；对 price 改用分位数截尾（如 0.5% ~ 99.5%），对尺寸变量用 IQR 规则；剔除比例超过 8% 就退回去用软处理——保留样本，加一列is_outlier标记，在稳健回归MASS::rlm()里让权重自动降下来。

5.2 定类变量当连续变量：J→D 编码成 1~7 的代价

现象：为了“省事”，把 color 手工映射成 1~7、clarity 映射成 1~8 直接入模，R² 看着还行，但残差图出现明显分层，同一颜色档内残差系统性偏正或偏负。

原因：这种编码强行假设“相邻等级之间的价格差相等”。实际上从 D 到 E 的溢价和从 I 到 J 的溢价完全不同，等差假设一旦不成立，系数就是有偏的，而且残差里会残留等级信息。

解决：老老实实因子化，让每个等级有自己的哑变量；如果确实需要降维，用有序因子的多项式对比（contr.poly），它假设的是“趋势”而不是“等距”，解释成本也更低。

5.3 只看 R² 选模型：交互项模型的解释成本

现象：fit_step 的 R² 0.98，fit2 只有 0.88，于是报告全篇用交互项模型下结论，写“理想切工使克拉溢价下降 XXX 美元”，被追问这个数字怎么算出来时答不上来。

原因：R² 衡量的是拟合，不是可解释性。交互项模型的系数是“某个等级相对基准等级的斜率修正量”，单独看某一项没有业务含义，必须组合起来算边际效应。

解决：分模型定用途。解释用无交互项模型，预测用交互项模型；非要用交互项下结论，就用emmeans或marginaleffects算特定组合下的边际效应再讲，别直接念系数。

5.4 先清洗再划分：数据泄漏

现象：报告里写“模型在测试集上 RMSE 很低”，但换一批新数据一跑，误差翻倍。

原因：异常值阈值、分位数截尾的边界都是拿全量数据算出来的，测试集的信息已经通过清洗口径渗进了训练过程。样本量五万多时这种泄漏不明显，一旦做交叉验证就会暴露。

解决：把清洗参数（3σ 的均值方差、分位数的分位点）只在训练集上估计，再原样应用到测试集；用rsample或caret的preProcess把这套流程封装成可复现的预处理对象。

5.5 外推到 3 克拉以上：模型没见过的区间

现象：拿模型预测一颗 4 克拉钻石的价格，结果偏得离谱，甚至出现负数区间。

原因：diamonds 里 3 克拉以上的样本占比极低，回归线在那个区间基本由少数几个点决定，杠杆值很高。加上 4.1 提到的异方差，高价段的预测区间本来就宽得没意义。

解决：先在报告里画出 carat 的样本分布，标出“有效预测区间”；预测超出区间时不给点估计，只给“超出建模范围”的结论。要真做高价段预测，就单独建一个高克拉子模型，或者改用分位数回归看区间而不是看均值。

6. 验收与答辩：五个检查点证明模型站得住

一套分析能不能交，不看图有多漂亮，看能不能通过下面五个检查点。我习惯在报告最后一页放一张自检表，答辩前逐条打钩，哪条没过就回头补，比临场解释省事得多。

第五项最容易被忽略，也最容易被问到。整体 RMSE 是一个被平均值掩盖的数字，必须拆开看：

# 在训练集上拟合，在留出的测试集上评估，避免自评自夸 set.seed(2024) idx <- sample(nrow(diamonds_clean), size = 0.7 * nrow(diamonds_clean)) train <- diamonds_clean[idx, ] test <- diamonds_clean[-idx, ] fit_final <- lm(price ~ carat + cut + color + clarity + depth + table, data = train) pred <- predict(fit_final, newdata = test) # 分段 RMSE：把预测误差按真实价格三分位切开看 brk <- quantile(test$price, probs = c(0, 1/3, 2/3, 1)) test$seg <- cut(test$price, breaks = brk, include.lowest = TRUE) rmse <- function(a, b) sqrt(mean((a - b)^2)) tapply(seq_len(nrow(test)), test$seg, function(i) rmse(test$price[i], pred[i]))

跑出来通常能看到：低价段 RMSE 小，高价段 RMSE 大出好几倍。这不是模型失败，而是必须写进结论的事实——“本模型在 5000 美元以下区间的预测较可靠，高价段仅作趋势参考”。答辩时主动说出这句话，比等着被问出来强太多。

还有一个小技巧：把 3.2 里那张 Ideal 切工价格偏低的小提琴图留在答辩 PPT 里，它是全篇最能体现“数据说话”的一页，也是解释为什么必须加交互项的入口。讲图比讲系数容易让人信服。

从那以后我每次做完回归，都会强制走一遍这套流程：清洗日志留档、三个模型并列、四项诊断齐全、误差分段报出来——少一步就不往报告里写结论。希望帮到你。

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