R语言钻石价格预测：多元线性回归与数据清洗实战

简介：围绕ggplot2自带diamonds数据集展开的钻石价格分析实战文档，面向已掌握R基础语法、希望系统练习数据清洗与多元线性回归的学习者。内容从54000余条钻石记录入手，覆盖缺失值与重复值检查、箱线图异常值识别与剔除、价格分布直方图及核密度解读，以及切工、颜色、净度、克拉重量等变量与价格关系的可视化探索，并逐步过渡到变量类型转换、新指标构造和多元线性回归建模预测。资源包共1个docx文件，约2.4MB，以图文混排方式呈现分析过程、代码结果与结论解读，便于对照复现和整理笔记。已有284人学习，适合作为回归分析与机器学习入门的完整案例，也可为课程作业、数据分析报告撰写提供参考。

1. 54000 颗钻石的价格账本：R 数据分析能问到什么程度

钻石价格分析这件事，用 R 语言做，最顺手的入口是 ggplot2 包里自带的 diamonds 数据集：53940 条记录、10 个变量，carat、cut、color、clarity、depth、table、price 加 x、y、z 三个尺寸，一条不缺。它干净、够大，又保留着真实业务里的脏东西——尺寸为 0 的错行、被 IQR 判成离群的高价样本、ordered factor 默认跑偏的对比矩阵，这些都会在建模时一个接一个冒出来。这份分析要回答两个问题：哪些因素真正在推动价格，以及给定一颗钻石的 4C 参数能把价格估到什么精度。适合刚学完线性回归想找完整项目落地的人，也适合要给业务方解释「凭什么这颗贵」的人。

2. 数据体检：缺失值、重复值、异常值三道关怎么过

拿到一份 5 万多行的数据集，很多人第一反应是直接lm()一把梭。我一般会先花二十分钟做体检，因为回归分析里最难补的后悔药，就是一开始没看清数据形状。这一章把加载、缺失、重复、异常值、分布形态拆成四步，每步都给你能直接抄的代码和判断标准。

2.1 装包与加载：三个包就够跑通全流程

这套流程不需要重型依赖。ggplot2 提供数据本身和作图能力，dplyr 负责分组汇总，mice 用来把缺失值模式画成图——即便最后确认没有缺失值，这张图也值得出一遍，它是写进报告里的证据。

# ggplot2 既是数据源也是绘图主力，mice 只用来做缺失模式可视化 install.packages(c("ggplot2", "dplyr", "mice", "car")) library(ggplot2) library(dplyr) library(mice) # diamonds 是 ggplot2 的内置数据集，data() 加载后直接可用 data(diamonds, package = "ggplot2") dim(diamonds) # 期望 53940 行 10 列 str(diamonds) # 看每一列的类型：三个 factor，七个 numeric

dim()返回的行列数是第一道校验，如果数字对不上，说明你加载的不是原版数据集，后面所有结论都要打问号。str()的重点是看类型：cut、color、clarity应该是 Factor，其余七个是 num。如果carat或price显示成 chr，多半是被某次读入操作污染了，得先as.numeric()转回来再往下走。

2.2 缺失值与重复值：先证明它干净，再动手

缺失值检查不要只写is.na()就完事，得把结果摊开看。逐列求和能告诉你哪一列有洞，md.pattern()能告诉你洞和洞之间是不是同时出现——这两件事在后续填补策略上完全不同。

# 逐列统计缺失数量，输出应全为 0 colSums(is.na(diamonds)) # 缺失模式矩阵图：行是观测，列是变量，红色方块代表缺失 md.pattern(diamonds, plot = TRUE) # 重复行检查：按整行完全一致来判定 sum(duplicated(diamonds))

colSums()的输出全 0，说明没有缺失；md.pattern()会把每一列画成一条，左侧数字是该模式下的观测数，右侧是该模式含缺失的变量个数。这套数据里通常只会出现一行「53940 0」，即所有行都完整。duplicated()返回 FALSE 的布尔向量，求和为 0 就代表没有完全重复的记录。

这里有个容易踩的细节：duplicated()判的是整行完全相同，而钻石数据里有大量「carat 价格参数都接近但不完全一样」的记录，它们不是重复，是不同实物。不要因为看着像就拿distinct()去重，那是把真实样本删掉。

2.3 异常值：1.5 倍 IQR 不能直接套在价格上

箱线图那套Q1 - 1.5*IQR到Q3 + 1.5*IQR的规则，教科书上讲得很顺，套到这份数据上会翻车。先看代码，再说为什么。

# 通用 IQR 边界函数，k 默认 1.5，可调成 3 做保守剔除 iqr_bound <- function(v, k = 1.5) { q <- quantile(v, c(0.25, 0.75), na.rm = TRUE) r <- q[2] - q[1] c(low = q[1] - k * r, high = q[2] + k * r) } # 分别看价格和三个尺寸的边界 sapply(diamonds[, c("price", "x", "y", "z")], iqr_bound) # 用 boxplot.stats 直接捞出被判为离群的价格 out_price <- boxplot.stats(diamonds$price)$out length(out_price) range(out_price)

sapply()会把四组上下界打成矩阵，boxplot.stats()$out直接返回被判定为离群的具体数值。跑完你会发现，价格这一列被判出局的数量以千计，而它们的下限不过一千多美元——这些恰恰是最贵的、最有商业价值的那批钻石。

所以在价格上按 IQR 剔除，等于亲手把模型的右尾砍掉，训练出来的模型在 1.5 克拉以上会系统性低估。我的做法是：价格一律不按 IQR 删，改用对数变换把右偏压平；真正要删的是物理上不可能的尺寸记录。

# 尺寸为 0 的记录：物理上不存在，且会让 depth 计算出 Inf bad_dim <- subset(diamonds, x == 0 | y == 0 | z == 0) nrow(bad_dim) summary(bad_dim[, c("carat", "price", "x", "y", "z")]) # 只保留三个尺寸都大于 0 的记录 df <- subset(diamonds, x > 0 & y > 0 & z > 0) nrow(df)

bad_dim通常只有个位数行，但破坏力不小：depth的定义是2*z/(x+y)*100，分母为 0 时直接产生 Inf 或 NaN，回归会整段报错。用subset()三条件同时过滤，删完再做一次nrow()确认剩余量级仍在 53900 以上，说明只动了该动的部分。

提示：删异常值前先想清楚「这条记录是录入错误，还是真实的极端业务样本」。前者删，后者留。

2.4 价格分布：右偏、双峰和 log 变换的取舍

价格的直方图是整份分析里信息密度最高的一张图。它同时告诉你三件事：分布形状、峰值位置、以及该不该做变换。

ggplot(df, aes(x = price)) + geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "grey75", colour = "white") + geom_density(colour = "firebrick", linewidth = 0.9) + geom_vline(xintercept = mean(df$price), linetype = "dashed") + geom_vline(xintercept = median(df$price), linetype = "dotted") + labs(x = "价格（美元）", y = "密度")

after_stat(density)是新版 ggplot2 的写法，老教程里的..density..会弹弃用警告，虽然还能跑，但别在新代码里用。bins = 60是经验值，5 万多条数据给到 60 组左右，既能看出峰值又不至于碎成噪声。两条竖线一条均值一条中位数，两者差距越大，右偏越重。

跑完你会看到均值明显在中位数右侧，分布长尾拖到一万多美元。这个形状直接决定了后面两件事：第一，残差正态性假设很难满足，QQ 图尾部一定翘；第二，如果做对数变换，log(price)会把这个形状压得接近钟形，模型残差表现会好很多。

但 log 变换不是免费的。系数从「每增加一单位 X，价格变动多少美元」变成「价格变动百分之多少」，给业务方讲的时候要多绕一道弯；并且反变换回原尺度时，exp()之后的预测区间并不对称，直接取中位数会对不上均值。我的习惯是两套模型都建，用原始尺度的模型讲系数，用对数尺度的模型做预测和区间，最后在测试集上用同一个指标比。

到这里，数据体检的四道关就过完了。剩下的 53900 多条记录，类型干净、无缺失、尺寸合法，可以进入变量工程。

3. 变量工程与可视化：把 cut、color、clarity 变成模型吃得下的形状

回归模型只认识数值，不认「理想切工」这种字面值。把三个定性变量编码进模型这一步，看着简单，实际藏着这份分析里最容易翻车的一个坑。这一章讲清楚因子水平的顺序怎么定、哑变量的参考组怎么选、分组对比图怎么看，以及怎么用非参数检验给「不同切工价格有显著差异」这句话提供证据。

3.1 有序因子：排序和建模是两回事

三个定性变量都有天然顺序。切工从 Fair 到 Ideal 逐级变好，颜色从 J 到 D 逐级变白，净度从 I1 到 IF 逐级变纯。用factor()加上levels和ordered = TRUE，能把顺序固化下来，图表的横轴就会按业务逻辑排，不会按字母序乱跳。

df <- df %>% mutate( cut = factor(cut, levels = c("Fair", "Good", "Very Good", "Premium", "Ideal"), ordered = TRUE), color = factor(color, levels = c("J", "I", "H", "G", "F", "E", "D"), ordered = TRUE), clarity = factor(clarity, levels = c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF"), ordered = TRUE) ) # 看一眼有序因子在建模时会被编码成什么样 contrasts(df$cut)

最后一行contrasts()是这套流程里我建议每个人都跑一次的检查。有序因子默认使用多项式对比（contr.poly），输出会是.L、.Q、.C这样的列名——线性、二次、三次趋势。这意味着模型给你的不是「Ideal 比 Fair 贵多少」，而是「切工等级每升一级的线性趋势贡献」。

对做预测影响不大，对讲结论就是灾难：业务方问「Premium 和 Ideal 差多少」，你没法从.L的系数里直接读出来。所以我的处理是：有序因子只用于排序和画图，建模前统一转回无序因子，让 R 按 treatment 对比生成哑变量。

# 保留 levels 顺序，但取消 ordered 属性，让 lm 生成哑变量 df_model <- df %>% mutate(across(c(cut, color, clarity), ~ factor(.x, levels = levels(.x), ordered = FALSE))) contrasts(df_model$cut) # 想以 Ideal 为参照组，显式 relevel，否则默认拿第一水平当基准 df_model <- df_model %>% mutate(cut = relevel(cut, ref = "Ideal"))

across()配合levels(.x)的写法，能保住业务顺序的同时换掉对比矩阵。relevel()那一步常被忽略：R 默认拿第一水平当参考组，也就是 Fair。结果就是所有切工系数都在讲「比 Fair 贵多少」，而 Fair 在全量数据里只占 3% 左右，这个基准既不稳也不好讲。改成 Ideal 之后，系数直接变成「比理想切工便宜多少」，读起来顺得多。

3.2 分组对比图：小提琴图比箱线图多告诉你什么

价格和切工的关系有个反直觉现象，第一次做这份分析的人几乎都会被绊一下：理想切工的钻石，价格中位数反而低于优质切工。

ggplot(df_model, aes(x = cut, y = price, fill = cut)) + geom_violin(trim = FALSE, alpha = 0.55) + geom_boxplot(width = 0.12, outlier.shape = NA, fill = "white") + scale_y_log10() + labs(x = "切工", y = "价格（美元，对数轴）") + theme(legend.position = "none")

trim = FALSE让小提琴的两端保留完整尾部，能看出分布厚度；scale_y_log10()把右偏压平，否则箱子全挤在底部看不清；outlier.shape = NA关掉箱线图自带的离群点，避免和小提琴的点重复。

图上能读出的信息比箱线图多一层：不只是中位数高低，还有各组的分布形状。理想切工那组的分布明显在低价比其他组更厚，说明大量小克拉钻石被切成了理想切工——这本身是合理的商业行为，小石头切坏了损失小，切好了溢价高。

要解释「中位数反而低」这个矛盾，直接把平均克拉数拉出来对照：

df_model %>% group_by(cut) %>% summarise( n = n(), median_price = median(price), mean_carat = round(mean(carat), 3), mean_x = round(mean(x), 3) )

跑完会看到 Ideal 组的mean_carat是最小的。这就把因果链说清楚了：不是理想切工不值钱，是理想切工这批样本整体克拉数偏小。这是典型的混杂变量效应，也是回归分析存在的理由——单变量对比会骗人，控制住 carat 之后切工的真实溢价才会显现。

同样的画法换到 color 和 clarity 上，方向就符合直觉了：颜色从 J 到 D，价格中位数稳定上移；净度从 I1 到 IF 也整体上移，但到了 VVS1 之后曲线会变平，说明顶级净度的边际溢价开始递减。这三张图摆在一起，就是后面建模时「为什么 clarity 要保留、为什么要加交互项」的直观依据。

3.3 显著性检验：为什么这里不该用方差分析

「不同切工的价格有显著差异」这句话得有统计证据。教科书默认走 ANOVA，但 ANOVA 的两个前提在这份数据上都不成立：正态性和方差齐性。

# 正态性：样本量 5 万时 shapiro 几乎必然拒绝，参考价值有限 tapply(df_model$price, df_model$cut, function(v) shapiro.test(v[1:5000])$p.value) # 方差齐性 bartlett.test(price ~ cut, data = df_model)

tapply()分组跑正态性检验，注意每组只取前 5000 个——shapiro.test()的样本上限就是 5000，超了直接报错。结果多半是 p 值小到打印成 0，但这里要清醒：样本量这么大，任何微小偏离都会被判定为显著，正态性检验在这个规模下几乎失去判别力。真正该看的是 QQ 图和偏度。

bartlett.test()对价格做方差齐性检验，右偏分布下各组方差差异通常也是显著的。两个前提一破，F 检验的 p 值就不可信了。换成不依赖分布假设的 Kruskal-Wallis：

kruskal.test(price ~ cut, data = df_model) # 显著之后做两两比较，bonferroni 控制多重比较带来的假阳性 pairwise.wilcox.test(df_model$price, df_model$cut, p.adjust.method = "bonferroni")

kruskal.test()是秩和检验的多组版本，只比较中位数位置，对分布形状不敏感。pairwise.wilcox.test()会输出一张下三角 p 值矩阵，p.adjust.method参数负责多重比较校正——五组两两比要跑 10 次，不校正的话假阳性率会从 5% 膨胀到 40% 左右。bonferroni 最保守，也可以换"BH"控制 FDR，样本量大时 BH 更实用。

这两步跑完，你手上有了一张完整的证据链：切工、颜色、净度三组变量，组间差异显著，且差异方向与业务认知一致。这才是后面敢把三个变量都塞进回归模型的前提。

3.4 共线性初筛：depth 其实是 x、y、z 的派生变量

建模前还有一件事必须做：看数值变量之间的相关结构。这份数据集里埋了一个恒等式，不注意的话会让模型系数彻底没法解释。

num_vars <- c("carat", "x", "y", "z", "depth", "table", "price") round(cor(df_model[, num_vars]), 3) # 验证 depth 的定义：depth = 2*z/(x+y)*100 df_model %>% summarise(max_err = max(abs(depth - 100 * 2 * z / (x + y))))

cor()输出的矩阵里，carat 和 x 的相关系数通常在 0.97 以上，x 和 y 也在 0.95 以上——这是物理必然，克拉数就是体积换算出来的。第二段代码验证的是depth的公式恒等性，max_err会非常接近 0，说明 depth 完全由 x、y、z 决定，没有独立信息。

这个发现直接影响建模策略：carat、x、y、z、depth 五个变量同时进模型，等于把同一份信息喂了五遍，VIF 会爆表，系数符号甚至可能反过来。我的取舍是主模型用carat + cut + color + clarity + depth + table，把 x、y、z 拿掉；如果确实想利用尺寸信息，就单独建一个含尺寸的版本做对比，两个模型的 R² 和残差摆在一起看，别混着解释。

到这里变量工程收尾，进入真正的建模环节。

4. 多元线性回归：从基准模型到交互项的四步迭代

建模这部分最怕两种写法：一种是一上来就price ~ .全变量梭哈，跑出个 R² 就收工；另一种是死磕 R² 数字，反复调变量直到 0.99 还觉得自己赢了。这一章按「基准 → 加交互 → 逐步回归 → 诊断 → 对数模型」的顺序推，每一步都说清楚为什么这么走，以及哪一步该停下来。

4.1 基准模型：R² 0.88 到底算不算好

先建一个不加任何交互项的模型，只用 4C 加 depth、table。

fit_base <- lm(price ~ carat + cut + color + clarity + depth + table, data = df_model) summary(fit_base) # 把三个尺寸变量也放进去，看 R² 能提多少 fit_size <- lm(price ~ carat + cut + color + clarity + depth + table + x + y + z, data = df_model) summary(fit_size) AIC(fit_base, fit_size)

summary()输出里重点看三处：Estimate列的系数符号和量级、Pr(>|t|)列的显著性、以及底部的Multiple R-squared和Adjusted R-squared。第一个模型的 R² 大致落在 0.85 上下，加上 x、y、z 之后能到 0.88 附近——这和原文记录的 0.8798 是一个量级。

数字本身说明不了太多。R² 的含义是「价格的总变异里，被这组变量解释掉的比例」，0.88 意味着还剩 12% 的变异没被覆盖。放在钻石这个场景里，那 12% 里装的是品牌溢价、证书机构、销售渠道、切割比例的美学评价——这些变量这份数据集里根本没有，所以 0.88 基本就是这套变量组合的天花板，再往上抠就是过拟合。

判断模型好不好，我的标准不是 R² 的绝对值，而是三件事：残差有没有明显结构（4.3 节展开）、调整 R² 和 R² 差多少（差太多说明变量冗余）、以及交叉验证后的误差是否稳定（第 6 章展开）。AIC()那行是给两个嵌套模型做比较用的，值越小越好，但差值小于 2 时基本可以认为没有实质差别。

Adjusted R-squared这个指标值得单独说一句：它会对变量个数做惩罚，所以在往模型里塞变量时，Adj R² 不升反降就是个明确信号——新加的变量没带来信息，只带来了复杂度。

4.2 交互项与逐步回归：carat 的斜率不是常数

基准模型有个隐含假设：克拉数每增加一单位，价格增加固定的美元数。这个假设在钻石上明显不成立——1 克拉到 1.1 克拉的涨价幅度，和 3 克拉到 3.1 克拉完全不是一个量级。

# 三个交互项：克拉数对价格的影响，随切工/颜色/净度而变 fit_int <- lm(price ~ carat * cut + carat * color + carat * clarity + depth + table + x, data = df_model) # 向后逐步回归，按 AIC 准则逐个剔除贡献不足的项 fit_step <- step(fit_int, direction = "backward", trace = 0) summary(fit_step) AIC(fit_int, fit_step)

carat * cut这种写法等价于carat + cut + carat:cut，即主效应加交互项一起进，不要在公式里只写carat:cut——那样会漏掉主效应，模型设定就不完整了。三个交互项同时进，模型参数会明显变多，这也是为什么需要有下一步。

step()的三个参数各有讲究：direction = "backward"从全模型开始逐步删项，比 forward 更稳，因为起点就包含了所有候选；trace = 0关掉每一步的打印，否则控制台会刷出几十行；如果不设scope，它只在当前模型的项里做减法，不会凭空造出新交互。跑完拿AIC()对比两个模型，值小的一方胜出。

这里要泼一盆冷水：AIC 逐步回归不是万能钥匙。它做的是单路径贪心搜索，换一批数据、甚至换个随机种子（如果涉及抽样），选出来的变量集可能就变了。所以我一般把step()的结果当参考，不当终审——真正拍板的是交叉验证的误差，以及模型在业务上说不说得通。

原文的做法是三个交互项全部保留、所有变量都被选入。这在样本量 5 万多的前提下可以接受，样本量撑得住参数膨胀。如果你手上的数据只有几千行，同样的操作大概率过拟合，该砍就砍。

4.3 回归诊断：四张图分别该看什么

模型跑完不看诊断图，等于开车不看后视镜。plot()一次出四张，每张对应一个假设。

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1)) plot(fit_step, which = 1:4) par(mfrow = c(1, 1))

par(mfrow = c(2, 2))把画布切成两行两列，which = 1:4指定输出哪四张。mar参数调边距，默认边距下轴标签容易被裁掉。跑完记得把画布还原成 1×1，否则后面所有图都会被挤在小格子里——这个坑我踩过不止一次。

第一张 Residuals vs Fitted 看线性假设和方差齐性。理想状态是点云均匀铺在 0 线两侧，没有喇叭口、没有弯曲趋势。价格这种右偏变量，这里通常会看到一个明显的漏斗——残差随拟合值增大而扩散。漏斗形状说明存在异方差，普通最小二乘的置信区间会偏窄，系数显著性被高估。

第二张 Normal Q-Q 看残差正态性。点贴着对角线走是理想，两端翘起来是厚尾。5 万条数据下尾部翘几乎是必然的，别指望它完全贴合。真正会伤到推断的是中间段也明显偏离，那说明分布形状不是「近似正态」而是「压根不是正态」。

第三张 Scale-Location 是残差绝对值的平方根对拟合值作图，比第一张更容易看出方差的趋势线是否水平。

第四张 Residuals vs Leverage 看杠杆值和 Cook 距离。右上角和右下角的点是有影响力的观测——它们本身残差不一定大，但杠杆高，删掉一个能让整条回归线转动。虚线标出的 Cook 距离等值线之外的点值得逐条回看，尤其是价格极高、克拉数极大的那几条。

除了四张图，多重共线性还得单独查：

library(car) vif(fit_step)

vif()输出每个变量的方差膨胀因子，经验阈值是 10，超过就说明该变量能被其余变量线性组合出来。含 x、y、z 的模型里，这几个变量和 carat 的 VIF 会明显偏高，这也是 3.4 节建议主模型去掉它们的理由。

4.4 对数模型与预测：反变换那一步最容易错

价格右偏这么重，对数模型往往比原始尺度表现更好。

fit_log <- lm(log(price) ~ log(carat) + cut + color + clarity + depth + table, data = df_model) summary(fit_log)

这里有两个变换：因变量取log(price)，自变量取log(carat)。前者把右偏压平，让残差更接近正态；后者把克拉数和价格之间的幂律关系拉直——钻石价格大致正比于克拉数的某个幂次，取对数之后斜率就是这个幂次指数，直接可读。系数含义也随之变化：log(carat)的系数是弹性，即克拉数每增加 1%，价格平均变化百分之几。

预测的时候，predict()有三个参数必须搞明白：

newd <- data.frame( carat = 0.9, cut = factor("Premium", levels = levels(df_model$cut)), color = factor("G", levels = levels(df_model$color)), clarity = factor("VS2", levels = levels(df_model$clarity)), depth = 61.5, table = 58 ) # 对数尺度上的预测与区间 pred_log <- predict(fit_log, newdata = newd, interval = "prediction", level = 0.95) pred_log # 反变换回美元尺度 exp(pred_log)

newdata必须是数据框，列名和训练数据完全一致。interval = "prediction"给的是预测区间（覆盖单个新观测），不是置信区间（覆盖期望值）；前者一定比后者宽，报告里给哪个取决于你要回答的问题——估算某颗具体钻石的合理成交价，用 prediction。level控制置信水平，默认 0.95，做价格区间报价时可以放到 0.90 收窄一些。

最后exp()那一步是整套流程里最容易翻车的地方。对数尺度上的预测区间是左右对称的，exp()之后变成右偏，中点不再是中位数。严格意义上，要从对数正态反变换得到无偏的均值预测，得乘上一个exp(sigma^2/2)的修正因子。做点估计时这个差别可能只有几百美元，但做区间报价时不能忽略。

5. 踩坑排查：这份分析里最容易翻车的五个地方

这份数据集因为太常见，教程满天飞，但真正跑起来会遇到的坑，教程里基本不写。下面五条是我在不同场合反复见到的翻车现场，每条按「现象 → 原因 → 解决」写，遇到对得上的直接抄解决方案。

5.1 现象：summary 里冒出一堆 cut.L、cut.Q

跑完summary()发现切工的系数不是四行，而是cut.L、cut.Q、cut.C、cut^4这种名字，系数值还特别大。

原因是factor(..., ordered = TRUE)默认使用多项式对比矩阵，模型拟合的是「切工等级的线性趋势、二次趋势、三次趋势」，而不是各水平相对基准的偏移量。这是统计上合法的做法，但业务上没法解释。

解决办法是在建模前把 ordered 属性去掉，或者直接显式指定对比方式：

# 方案一：取消 ordered，保留 levels 顺序 df_model <- df %>% mutate(cut = factor(cut, levels = levels(cut), ordered = FALSE)) # 方案二：显式指定 treatment 对比 contrasts(df_model$cut) <- contr.treatment(n = 5, base = 5)

方案一更干净，推荐。方案二里base = 5表示以第五个水平（Ideal）为参照，等价于relevel()。

5.2 现象：剔完异常值，R² 漂亮了，但大克拉钻石全被低估

按 1.5 倍 IQR 把 price 的离群点删掉，模型 R² 从 0.86 涨到 0.92，看着很美。但拿几颗 2 克拉以上的真实报价去测，预测值普遍低几百到上千美元。

原因是价格分布的右尾不是噪声，是信号。IQR 规则假设数据近似对称，而价格是典型的重尾右偏分布，被判为「离群」的那些点恰恰承载了高端市场的信息。删掉它们，模型就只在低中价位段学过，外推到高端必然偏低。

解决办法是价格不做 IQR 剔除，改用两种手段处理尾部：一是对 price 取对数，把尾部压缩；二是改用分位数回归（quantreg::rq()），直接对中位数或 90 分位数建模，不受极端值拉扯。如果确实要删，只删尺寸为 0 或明显录入错误的记录。

5.3 现象：模型报错，或者出现负价格预测

lm()直接报NA/NaN/Inf in 'y'，或者模型能跑通，但预测出来的价格是负数。

原因是数据里有 x、y、z 为 0 的记录。depth的定义是2*z/(x+y)*100，分母为 0 时产生 Inf；如果某行 z 为 0 而 x、y 正常，depth 算出 0，这个值在物理上不可能，会把 depth 的系数往错误方向拉。

解决办法是先过滤再做任何计算：

df <- subset(diamonds, x > 0 & y > 0 & z > 0) # 顺手验证一遍 depth 的恒等式，误差应该接近 0 max(abs(df$depth - 100 * 2 * df$z / (df$x + df$y)))

第二行是个便宜又好用的自检。如果误差超过 0.1，说明 depth 列被改过，或者单位不一致，得回头查数据来源。

5.4 现象：predict 报错「factor has new levels」

训练集跑得好好的，拿新数据predict()直接抛错，提示因子出现了新水平。

原因是newdata里的因子水平和训练时不一致。哪怕字面值一样，只要水平集合不同（少一个水平、或者顺序不同），R 就认为这是另一个因子。手写data.frame(cut = "Ideal")尤其容易出这个问题，因为默认会按字母序生成水平。

解决办法是把训练集的 levels 显式传给新数据：

newd <- data.frame( cut = factor("Ideal", levels = levels(df_model$cut)), color = factor("G", levels = levels(df_model$color)) )

养成这个习惯，所有进predict()的因子都走一遍factor(x, levels = levels(train$x))，能省掉大量排查时间。

5.5 现象：样本内 R² 高得离谱，上线后误差上千美元

模型在训练集上 R² 达到 0.95 以上，残差标准差只有几百，信心满满拿去预测新钻石，误差动辄上千。

原因通常是两个叠加：一是没有留出测试集，所有评估都在样本内做；二是模型里塞了 x、y、z 这类和 carat 高度共线的变量，它们在训练集上能靠微小的数值差异「记住」样本，泛化能力极差。样本内 R² 衡量的是拟合能力，不是预测能力，这两件事在变量冗余时差距会被放大。

解决办法是划分训练测试集，并且评估指标换到原始尺度上算：

set.seed(2024) idx <- sample(nrow(df_model), size = floor(0.7 * nrow(df_model))) train <- df_model[idx, ] test <- df_model[-idx, ] fit_tr <- lm(log(price) ~ log(carat) + cut + color + clarity + depth + table, data = train) pred <- exp(predict(fit_tr, newdata = test)) rmse <- sqrt(mean((pred - test$price)^2)) mape <- mean(abs(pred - test$price) / test$price)

注意exp()的位置：模型在 log 尺度上训练，预测出来必须反变换回美元，再和真实价格比。如果忘了反变换，RMSE 会小得离谱，因为你在拿 log 单位和美元比。

6. 进阶技巧：一克拉关口、交叉验证和误差的真实含义

这份数据里有个现象，第一次看到会以为是数据错误：把 carat 在 1.0 附近的点单独画出来，价格在 1.00 处有一个明显的跳变，0.99 克拉和 1.01 克拉之间能差出几百美元。

ggplot(subset(df_model, carat > 0.85 & carat < 1.15), aes(x = carat, y = price)) + geom_point(alpha = 0.12, size = 0.7) + geom_vline(xintercept = 1, linetype = "dashed", colour = "firebrick") + scale_y_log10() + labs(x = "克拉数", y = "价格（美元，对数轴）")

这不是数据问题，是消费心理：一克拉是个整数关口，跨过去就是另一个价格带。可惜线性模型捕捉不到这个跳变，carat的系数是一个平滑的斜率，1.00 和 1.01 在模型眼里几乎没有区别。

想把这个效应吃进模型，常见做法是加一个「是否跨过一克拉」的指示变量，或者按克拉数分箱建分段模型：

df_model <- df_model %>% mutate( carat_bin = cut(carat, breaks = c(0, 0.5, 1, 1.5, 2, 5), include.lowest = TRUE), over_one = as.integer(carat >= 1) ) fit_seg <- lm(log(price) ~ carat_bin + cut + color + clarity + depth + table, data = df_model) summary(fit_seg)

cut()在这里是 base R 的分箱函数，和 dplyr 的cut变量名撞了，写的时候要留意作用域——这是 R 里一个经典的命名陷阱，dplyr::cut并不存在，但变量名cut会遮蔽函数。稳妥写法是base::cut(carat, ...)。

over_one这个指示变量的系数，直接读出来就是「跨过一克拉的平均溢价比例」，比从平滑斜率里反推直观得多。carat_bin版本则更彻底，每个区间一个截距，代价是参数变多、边界处不连续。

交叉验证这块，我一般跑 5 折，比单次划分稳：

library(caret) set.seed(2024) ctrl <- trainControl(method = "cv", number = 5) fit_cv <- train( log(price) ~ log(carat) + cut + color + clarity + depth + table, data = df_model, method = "lm", trControl = ctrl ) fit_cv$results

trainControl(method = "cv", number = 5)指定 5 折交叉验证，method = "lm"表示内部还是线性回归，不做正则化。fit_cv$results会给出 RMSE、R²、MAE 在不同折上的均值和标准差——标准差这一项特别值得看，它反映模型在不同数据子集上的稳定性。标准差大的模型，哪怕均值好看，上线也得谨慎。

最后说一个我踩过的坑：评估指标一定要在原尺度上算。在对数尺度上算出来的 RMSE，单位是「log 美元」，数值看着很小很漂亮，但没有任何业务含义，也没法和实际报价误差对照。我现在的习惯是，所有汇报给别人的误差数字，一律先exp()回到美元，再用 RMSE 和 MAPE 两个指标交叉看——RMSE 会被高价钻石主导，MAPE 对低价钻石更敏感，两个一起看才知道模型在哪个价位段表现差。

从那以后我每次跑这类价格回归，都强制走一遍「划分训练测试 → log 尺度训练 → 反变换 → 原尺度算误差」的流程，中间任何一步偷懒都会让最后的结论失真。希望帮到你。

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