简介：「钻石价格分析-基于R语言」是一份面向R语言初学者与数据分析入门者的实操型文档，围绕ggplot2内置的diamonds数据集展开，尝试解决从数据清洗到回归建模预测价格这一完整分析链路中的学习需求，可用作统计学课程作业、课程设计或数据分析练手项目的参考范例。资源包内仅含1个docx文档，压缩后约2.4MB，内容集中于单一主题，便于一次性通读并对照复现。文档从数据探索切入，依次覆盖缺失值与重复值检查、用箱线图识别异常值并剔除、借助价格直方图与核密度曲线判断右偏分布，再延伸到把cut、color、clarity等定性变量编码为factor、构造新变量、做简单线性回归与多元回归的变量筛选；同时配有各属性与价格关系的小提琴图、透明度分组对比等可视化结果与解读，最后给出价格走势预测的思路。已有284人学习，适合希望在真实数据集上跑通探索、建模、解释全流程的读者。

1. 54000 颗钻石摆在一起，先看清这份 R 语言分析的真实难点

钻石价格分析这个题目，几乎每本 R 语言教材都拿 ggplot2 自带的 diamonds 数据集练一遍——54000 多条记录、10 个变量，看着像个玩具数据集。真把回归跑完就会发现不对劲：carat 和 price 的关系明显不是直线，cut 的系数跑出来是负的，R² 顶到 0.88 就上不去了。这份资源走的是完整链路，从缺失值、重复行、异常值清洗，到描述性统计、分组可视化，再到简单线性回归、带交互项的多元回归、非参数显著性检验和回归诊断。它适合已经能写几行 R、但没完整跑过一次「清洗—可视化—建模—诊断」的人；也适合想拿一个真实分布（右偏、长尾、混着定性变量）练回归的人。下面按我实际复现的顺序拆，代码能直接抄。

2. 数据装载与清洗：哪些行不该进回归模型

这一步决定后面所有结论的上限。diamonds 数据集本身是干净的，坑不在缺失值，而在因子类型和物理上不可能存在的尺寸记录。前者会让你的回归系数含义整个跑偏，后者会把 x、y、z 的系数拉到一个没法解释的位置。

2.1 装载与类型确认：ordered factor 是第一个陷阱

# diamonds 挂在 ggplot2 包里，先确认包在不在，不要无脑重装 if (!requireNamespace("ggplot2", quietly = TRUE)) { install.packages("ggplot2") } library(ggplot2) # 显式指定 package，避免工作区里有同名对象时拿到旧数据 data("diamonds", package = "ggplot2") dim(diamonds) # 行数、列数 str(diamonds) # 变量类型 + 因子水平 sapply(diamonds, class) # 只看 class，一眼扫完 10 个变量

逻辑说明：data()的第二个参数package是关键。R 的data()会先在工作区找、再去已加载的包里找，如果之前load()过一个叫 diamonds 的 rda，直接data(diamonds)拿到的东西可能不是你以为的那份。str()要看两件事——数值列的 num/int 类型对不对，因子列的水平数和顺序对不对。

参数说明：quietly = TRUE让requireNamespace不打印加载信息，脚本source()的时候输出干净；install.packages()只在包缺失时触发，不会每次都跑。

跑完str()你会看到cut、color、clarity都是Ord.factor。这是第一个必须处理的地方。R 对有序因子（ordered factor）默认使用多项式对比（contr.poly），也就是把 5 个切工等级拆成线性、二次、三次、四次趋势项，回归结果里出现的是cut.L、cut.Q、cut.C这种名字。你想要的「Ideal 相对 Fair 贵多少」根本不在输出里。所以建模前必须显式重建：

dm <- diamonds # 重建为普通因子，levels 按业务从差到好手工排 dm$cut <- factor(dm$cut, levels = c("Fair", "Good", "Very Good", "Premium", "Ideal"), ordered = FALSE) dm$color <- factor(dm$color, levels = c("J", "I", "H", "G", "F", "E", "D"), ordered = FALSE) dm$clarity <- factor(dm$clarity, levels = c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF"), ordered = FALSE) # 确认对比方式变成 treatment 对比 contrasts(dm$cut)

逻辑说明：levels的第一个元素就是回归里的基准组。cut 的基准是 Fair（最差切工），color 的基准是 J（最差颜色），clarity 的基准是 I1（最差净度）——三个都选最差档做基准，好处是所有系数都是正数，读起来就是「比最差档贵多少」，方向和业务直觉一致。

参数说明：ordered = FALSE不能省。哪怕你levels排得再对，只要保留 ordered 属性，lm()照样走多项式对比。我见过不少人的模型里cut.L显著、cut.Q不显著，然后对着输出发半天呆，就是这么来的。

2.2 缺失值与重复行：先量化，再决定删不删

library(mice) # 缺失值模式矩阵：最后一列是每行缺失个数，最后一行是每列缺失个数 md.pattern(dm, rotate.names = TRUE) # 不装 mice 的话用这两行也够 colSums(is.na(dm)) sum(duplicated(dm))

逻辑说明：md.pattern()返回的矩阵里，1 表示有值、0 表示缺失，最后一行统计每个变量的缺失条数，最后一列统计每条记录的缺失变量数。diamonds 这个数据集跑出来全是 1，说明没有缺失值。duplicated()返回逻辑向量，sum()得到重复行数。

参数说明：rotate.names = TRUE把变量名转 90 度，10 个变量的表横着排屏幕塞不下，转过来才读得清。

重复行这里有个判断要做。diamonds 里确实存在属性完全一致的行，但它们可能是两颗独立打磨的钻石——克拉、切工、颜色、净度、尺寸都撞上并不是不可能。这份资源的做法是直接去重，我的做法是先记数量再决定：

n_before <- nrow(dm) dm <- dm[!duplicated(dm), ] cat("去重前:", n_before, "去重后:", nrow(dm), "差异占比:", round((n_before - nrow(dm)) / n_before * 100, 3), "%\n")

逻辑说明：:-号在[` 里配合逻辑向量是标准的「保留非重复行」写法。打出去重前后的行数和差异占比，用来判断这个操作要不要写进最终流程。

参数说明：差异占比在千分之几的量级，删掉对结论没影响；如果超过 1%，就要考虑这些「重复」是不是清洗阶段的副产物（比如你把连续变量四舍五入过），那就不能删。

2.3 异常值判定：箱线图只能告诉你「有」，不能告诉你「删哪些」

# 第一步：物理规则过滤，三维尺寸必须为正 dm <- subset(dm, x > 0 & y > 0 & z > 0) # 第二步：尺寸上限，y 最大 58.9mm 这类是典型的小数点录入错误 dm <- subset(dm, x < 20 & y < 20 & z < 20) # 第三步：对 price 只做标记，先不删 q <- quantile(dm$price, c(0.25, 0.75)) iqr <- q[2] - q[1] dm$flag_price_out <- dm$price < (q[1] - 1.5 * iqr) | dm$price > (q[2] + 1.5 * iqr) table(dm$flag_price_out)

逻辑说明：这里刻意拆成三步而不是一行subset()。第一步去掉尺寸为 0 的记录——0 毫米的钻石不存在，这类行通常是缺测被填成 0 的结果。第二步卡上限，原始数据里 y 的最大值到了 58.9 毫米，而 carat 最大才 5.01，一颗 5 克拉的钻长边通常在 10 毫米量级，58.9 是数量级错误。第三步对价格只打标记不删除，因为 price 是右偏长尾，大克拉高价钻恰恰是这个数据集最有信息量的部分，一刀切会把模型往中间压缩。

参数说明：1.5 * iqr是箱线图触须的标准倍率，对应「离群但不算离谱」的边界；想更严格就换3 * iqr。flag_price_out是一个新列，后面做敏感性分析时用它分层重跑一遍模型，看系数有没有大变化，比直接删掉稳妥得多。

对比一下原文用三倍标准差剔除异常值的做法。在正态假设下 3σ 是合理的，但 price 明显不服成正态分布（下一章会看到），用正态的尺子去量长尾分布，等于系统性地砍掉高价样本，最后模型会低估大克拉钻石的价格。所以我一般对价格类变量用 IQR 标记 + 分层验证，对物理尺寸用硬规则剔除。

2.4 衍生变量与建模底表：把清洗流程固化成一次

# 体积：三个维度的乘积，和 carat 高度相关但不等价 dm$volume <- dm$x * dm$y * dm$z # 每克拉单价，用来观察「越大越贵」是不是线性 dm$price_per_carat <- dm$price / dm$carat # 把清洗结果落盘，后面所有分析都从这份快照读 saveRDS(dm, "diamonds_clean.rds")

逻辑说明：volume和x、y、z之间存在极强的线性关系，建模时这两组变量二选一，同时放进去会让系数标准误爆炸。price_per_carat是个诊断变量，不是建模变量——把它按 carat 分箱画出来，能直接看到每克拉单价在 1 克拉、1.5 克拉附近的变化，比看回归系数直观。

参数说明：saveRDS()存单个对象，读回来用readRDS()，比write.csv()快且不会丢因子水平。清洗脚本从头到尾跑一遍只要几秒，但如果你每次分析都重跑清洗、每次都微调一下过滤阈值，前后结论就对不上了——固化快照是保证可复现的最低成本手段。

这一步别用attach()。attach(diamonds)之后工作区里凭空多出 10 个变量名，price到底是原始的diamonds$price还是你改过的dm$price，你自己三天后也说不清，这是黑匣子级的事故源。

3. 描述性统计与可视化：价格右偏、1000 美元峰值和四张必看图

清洗完先别急着lm()。这个数据集的形状决定了后面模型怎么设：价格右偏到什么程度、峰值在哪、分组之后分布有没有重叠、carat 和 price 到底是直线还是分段。这四件事没搞清楚就上回归，系数解释一定出事。

3.1 价格分布三件套：直方图、核密度与正态曲线叠加

library(ggplot2) ggplot(dm, aes(x = price)) + geom_histogram(aes(y = after_stat(density)), bins = 60, fill = "grey75", colour = "white") + geom_density(colour = "#c0392b", linewidth = 0.8) + stat_function(fun = dnorm, colour = "steelblue", linetype = "dashed", args = list(mean = mean(dm$price), sd = sd(dm$price))) + labs(x = "价格（美元）", y = "密度") + theme_minimal()

逻辑说明：after_stat(density)把直方图的纵轴从计数换成密度，这样核密度曲线和正态曲线才能画在同一根轴上比较。三条线叠在一起的目的很明确——灰柱是真实分布，红线是平滑后的形状，蓝虚线是「如果它服从正态会是什么样」。三者一比，右偏有多严重一目了然。

参数说明：bins = 60是对 5 万多条记录比较合适的粒度，太少看不出断层，太多全是毛刺。linewidth在新版 ggplot2 里替代了旧的size，写size会报警告。stat_function的args必须传list()，直接写mean = ..., sd = ...会被当成未知参数。

跑完能看到三件事。第一，峰值在 1000 美元附近，0 到 2000 美元区间堆了绝大多数样本。第二，2000、4000、8000 美元附近有明显断层——这是定价习惯造成的，整数关口附近的样本密度会突然掉下去。第三，正态曲线和核密度曲线在 8000 美元之前几乎完全不重合，之后才逐渐贴近横轴。这意味着任何假设 price 正态的检验都不成立，后面做显著性检验只能走非参数路线，这也是原文选非参数检验的原因。

3.2 分组比较：小提琴图为什么比箱线图更适合 cut

ggplot(dm, aes(x = cut, y = price, fill = cut)) + geom_violin(trim = FALSE, scale = "width", alpha = 0.7) + geom_boxplot(width = 0.12, outlier.shape = NA, fill = "white") + scale_y_log10() + theme_minimal() + theme(legend.position = "none")

逻辑说明：小提琴图给的是完整的密度形状，箱线图给的是中位数和四分位。两张图叠在一起，既能看出「Ideal 组是不是整体更贵」，也能看出「Ideal 组内部是不是有两坨」。这个数据集里 Ideal 组确实存在一批低价样本，只看中位数会被平均值带跑。

参数说明：trim = FALSE保留尾部，不然小提琴的两端会被截断，看不出长尾；scale = "width"让每组小提琴等宽，默认是按样本量缩放，Fair 组样本少会缩成一条缝，视觉上造成误导；scale_y_log10()是必须的，price 右偏，线性轴上所有箱子都压在底部；箱线图里outlier.shape = NA是因为小提琴已经表达了尾部，再画一遍离群点会糊成一团。

同样套路换color和clarity再跑两遍。颜色从 J 到 D，价格中位数基本单调上升，这个方向符合直觉。净度那边就有意思了——IF、VVS1、VVS2 价格最高，但 I1 的中位数并不是最低的，SI2 反而比 I1 低。原文对这个现象的解释是「净度到达优质状态后，区分钻石价格的维度更多」，说白了就是净度不是唯一变量，低净度区间里 carat 的权重更大。

3.3 克拉与价格：线性坐标会骗人，换对数轴再看

ggplot(dm, aes(x = carat, y = price)) + geom_point(alpha = 0.15, size = 0.6, colour = "grey30") + geom_smooth(method = "lm", formula = y ~ x, colour = "red", linewidth = 0.8) + scale_x_log10() + scale_y_log10() + labs(x = "克拉（对数轴）", y = "价格（对数轴）") + theme_minimal()

逻辑说明：线性坐标下 5 万多个点全挤在左下角，什么都看不出来。双对数轴上能看到三层结构——整体趋势是直线，说明 carat 和 price 之间接近幂律关系；在 0.5、1.0、1.5 克拉处有垂直的「墙」，说明整数克拉有溢价，定价不是连续的；大克拉区域点的离散度明显变大，说明高价钻的定价更依赖其他因素。

参数说明：alpha = 0.15是 5 万个点的合理透明度，再高就成实心块了；size = 0.6让单点足够小，避免重叠成色块；formula = y ~ x显式写给geom_smooth，新版 ggplot2 不写会警告。

这张图直接决定了回归怎么设。如果只用price ~ carat，残差图一定有明显的弧形，因为真实关系是对数级的。后面第 4 章会看到，单变量模型的 R² 已经不低，但残差结构不对，这就是原因。

3.4 相关性矩阵与共线性初判

num_vars <- c("carat", "depth", "table", "x", "y", "z", "price") round(cor(dm[num_vars], use = "complete.obs"), 3) # 只看尺寸和 carat 的共线性 library(car) vif(lm(price ~ carat + depth + table + x + y + z, data = dm))

逻辑说明：cor()给的是两两线性相关，vif()给的是「某个变量被其他变量解释的程度」，后者对多元回归更有意义。x、y、z 三者之间的相关系数普遍在 0.95 以上，VIF 会直接冲到几十。

参数说明：use = "complete.obs"是保险写法，虽然这个数据集没缺失，但脚本换个数据源就跑，加一行省事。VIF 的经验阈值是 10，超过就说明共线严重，系数符号可能反掉。

结果是 x、y、z 的 VIF 全部超标。这不奇怪，一颗钻的长宽深本来就同步变化。处理方式有两种：要么只留volume一个综合维度，要么留carat然后把 x、y、z 全砍掉。我一般选后者，因为 carat 是行业里真正的交易单位，系数解释起来有业务含义，而「体积每增加 1 立方毫米贵多少钱」没人听得懂。

4. 回归建模：从 carat 单变量系数 7760 到带交互项的多元模型

到这一步才开始写lm()。顺序不能乱：先单变量看方向和量级，再全变量看整体解释力，最后加交互项看边际效应是不是随分组变化。每一步都要留残差图，看结构有没有被吃掉。

4.1 简单线性回归：7760 美元/克拉这个数怎么读

fit1 <- lm(price ~ carat, data = dm) summary(fit1) # 系数置信区间，比只看 p 值信息量大 confint(fit1, level = 0.95) # 四张诊断图一次出 par(mfrow = c(2, 2)) plot(fit1) par(mfrow = c(1, 1))

逻辑说明：summary()里重点看三行——carat 的 Estimate、它的 Pr(>|t|)、还有最下面的 R²。这份资源里跑出来 carat 的系数约 7760，也就是「克拉数每增加 1，价格平均增加约 7760 美元」。截距是负的，但在业务上没有意义，因为不存在 0 克拉的钻石，别硬解释。

参数说明：confint()默认 95% 置信区间，可以改成level = 0.99。par(mfrow = c(2,2))把四张图排成 2×2，跑完记得用par(mfrow = c(1,1))复位，不然后面所有图都会被挤到一个小格子里——这是新手最常见的「图怎么画不出来」问题。

四张诊断图里，先看 Residuals vs Fitted。如果点云呈喇叭口（左窄右宽），说明存在异方差，残差的方差随拟合值增大，这时候系数的标准误会算小，p 值会虚高。再看 Normal Q-Q，点偏离对角线两端翘起，说明残差不正态——price 右偏，残差不可能正态。这两条都在提示同一件事：直接对 price 建线性模型，假设层面就不成立，后面才会去考虑对数变换和交互项。

4.2 哑变量处理：基准组选谁，系数含义就变了

# 确认三个因子的对比方式 contrasts(dm$cut) contrasts(dm$color) contrasts(dm$clarity) # 显式指定基准组，写成代码里比靠 levels 顺序更稳 dm$cut <- relevel(dm$cut, ref = "Fair") dm$color <- relevel(dm$color, ref = "J") dm$clarity <- relevel(dm$clarity, ref = "I1")

逻辑说明：contrasts()打印出来的是对比矩阵，行是水平名，列是哑变量名，第一行全 0 的那个水平就是基准组。relevel()的作用是把指定水平挪到第一位，也就是改基准。

参数说明：ref参数只接受已经存在于 levels 里的字符串，写错会直接报错。这一步我建议在建模前单独跑一次relevel()并写死，而不是依赖上游factor()里 levels 的排列顺序——上游改一行，下游所有系数含义全变，这种坑排查起来极其费时间。

为什么基准组这么重要？因为系数是「相对基准组的差值」。如果 color 的基准是 J（最差），那 D 的系数就是正的一大截，「颜色越好越贵」一眼就懂。如果按默认字母序，基准变成 D（最好），所有系数都变成负数，输出看起来像「颜色越差越贵」，实际上只是参照系反了。这个符号问题在下一章的避坑清单里还会再出现一次——它是这份数据集中最容易翻车的地方。

4.3 无交互项多元回归：R²=0.8798 说明什么

fit2 <- lm(price ~ carat + cut + color + clarity + depth + table + x + y + z, data = dm) summary(fit2)

逻辑说明：这是标准的加性模型，每个变量对 price 的边际效应是常数，不随其他变量变化。这份资源里 R² ≈ 0.8798，调整 R² 略低一点。也就是说还有 12% 的价格波动没被解释。

参数说明：lm()里公式写+就是加性，不写交互。注意这里同时放了carat和x、y、z，共线性严重，你会看到x的系数可能是负的——这不是「钻石越长越便宜」，而是「在 carat 和其他尺寸不变的前提下，x 单独变大反而降价」，业务上没意义。正确做法是砍掉 x、y、z，或者用上一章算的 volume 替代。

depth和table这两个变量在这个模型里通常不显著，符号还可能反。depth 是深度百分比，定义是2z/(x+y)，它和 x、y、z 是函数关系，放进来纯属添乱。我的处理是先在模型里保留它们看 p 值，不显著就从最终模型里去掉，别舍不得——每多一个不显著的变量，调整 R² 就掉一点。

R² = 0.8798 这个数怎么评价？单看绝对值不算低，但结合残差图看就不够。残差在低价区和高价区呈现不同的离散度，说明模型对便宜钻和贵钻的拟合能力不一样。这就是这份资源往下走、加交互项的原因。

4.4 交互项与 AIC 逐步回归：三个交互项全进模型之后

# carat 与三个定性变量的交互，主效应会被自动保留 fit3 <- lm(price ~ carat * (cut + color + clarity) + depth + table + x + y + z, data = dm) summary(fit3)$r.squared AIC(fit3) # 向后逐步，按 AIC 删变量 fit4 <- step(fit3, direction = "backward", trace = FALSE) summary(fit4) AIC(fit3, fit4)

逻辑说明：carat * (cut + color + clarity)是展开写法，等价于 carat、cut、color、clarity 四个主效应加上 carat:cut、carat:color、carat:clarity 三个交互项。交互项的含义是「克拉数对价格的影响，在不同切工/颜色/净度下不一样」——这正是加性模型缺的东西。

参数说明：step()的direction可选"backward"、"forward"、"both"，向后是从全模型出发逐步删；trace = FALSE关掉每一步的打印，不然几万行输出刷屏；向前选择时必须给scope参数指定候选变量范围。5 万多行数据、几十个参数，step()会反复拟合几十次模型，总耗时也就一两秒，行数再上一个量级就要考虑换方案了。

原文的结论是「所有变量都被选入模型」，这个结果本身值得琢磨——AIC 逐步回归在 5 万样本上非常保守，因为 AIC 的惩罚项是2k，而样本量一大，任何一个微小改进的似然增益都超过这个惩罚。所以别把「全都留下了」理解成「所有变量都重要」，它只说明在这么大样本下，删掉任何一个都会让 AIC 变差一点点。

交互项进来之后 R² 会往上走一点，但更值得看的是交互项的系数。如果carat:cutIdeal是负的，意思就是「同样增加一克拉，Ideal 切工的钻石涨价比 Fair 少」——这在业务上是说得通的，因为小克拉钻更强调切工，大克拉钻更强调保重。这类结论才是交互项真正的价值，光看 R² 涨了零点零几没什么感觉。

5. 踩坑排查：R 语言钻石回归里最容易翻车的五个地方

下面五条都是我在这份数据上真实撞过的。每条按「跑出来什么样 → 为什么 → 怎么改」写，遇到了直接对号入座。

5.1 cut 的系数是负的，是不是数据错了

现象：模型里cutIdeal的系数是负数而且显著，可你单独按 cut 分组算平均价格，Ideal 组一点也不便宜。两个结论打架，第一反应是数据有问题。

原因有两层。第一层是参照系——如果 cut 还是 ordered factor 没改过来，你看到的根本不是cutIdeal，而是cut.L、cut.Q这些趋势项，符号含义完全不同。第二层是混杂：Ideal 切工大量出现在小克拉钻上（切工好能让小钻显大），而 Fair 切工往往出现在大克拉钻上（保重优先）。当你控制住 carat 之后，Ideal 在同等克拉下反而便宜，这是典型的辛普森悖论式现象，不是数据错。

解决：先contrasts(dm$cut)确认对比方式，再tapply(dm$price, dm$cut, median)看分组中位数，两个都正常说明是混杂，不是错误。解释系数时必须带限定语——「在控制克拉重量、颜色、净度等其他变量不变的前提下」。想进一步验证就加carat:cut交互项，看不同克拉区间里 cut 的效应是不是反过来的。

5.2 对 price 取对数后，回归系数不能直接读成美元

现象：为了处理右偏，把模型改成lm(log(price) ~ log(carat) + ...)，跑出来 carat 的系数是 1.8，然后有人直接写「克拉数每增加 1，价格增加 1.8 美元」。

原因：双对数模型里系数是弹性，1.8 的意思是「carat 增加 1%，price 平均增加 1.8%」，跟美元没有半点关系。如果是半对数模型（只对 y 取 log），系数 β 对应的是「x 增加 1 单位，y 变化(exp(β)-1)×100%」，β 小的时候近似等于 β×100%，β 大的时候这个近似误差很难看。

解决：写个转换函数固定下来，别每次心算。

# 半对数模型系数转百分比 pct_change <- function(beta) (exp(beta) - 1) * 100 # 双对数模型的弹性直接读，但要写成百分比表述 fit_log <- lm(log(price) ~ log(carat) + cut + color + clarity, data = dm) summary(fit_log)$coefficients["log(carat)", "Estimate"]

参数说明：pct_change()传入的是summary()里的 Estimate 列，不是 t 值。双对数模型的弹性解释可以直接用系数值乘以 100 变成百分比，不需要再套exp()。

5.3 逐步回归结果每次跑都不一样

现象：同一份数据，同事那边step()选出来的变量集跟你不一样，或者你自己隔一天重跑，变量少了一个。

原因：step()本身不吃随机数，结果不稳定通常来自三个地方。一是起始模型不同——你从全模型出发，他从只含主效应的模型出发，贪心搜索的路径当然不同。二是样本量变了——上游某一步na.omit()或者过滤条件改过，参与拟合的行数不一样，AIC 没法比。三是工作区里有残留对象，lm()实际用的是旧的dm。

解决：把清洗结果saveRDS()落盘，建模脚本第一行dm <- readRDS("diamonds_clean.rds")，彻底切断对工作区状态的依赖。起始模型写成固定公式，别用~ .这种简写——~ .会把数据框里所有列都拉进来，包括你后面加的flag_price_out、price_per_carat这些诊断列，直接造成数据泄漏。每次跑完把AIC(fit3, fit4)和变量列表一起记在脚本注释里。

5.4 交互项全加进去之后报错或出现 NA 系数

现象：summary()里冒出几个NA的系数，或者直接报对比相关的错误，提示某个因子水平有问题。

原因：某个因子水平的组合样本量为 0。比如你按某个条件过滤之后，cut == "Fair"且color == "D"的记录一条都不剩，交互项在这个格子上就没有数据，设计矩阵秩亏，R 只能给 NA。另一个常见原因是把depth又按定义2*z/(x+y)算了一遍塞进模型，和已有的 x、y、z 构成精确共线。

解决：建模前先查空格子。

# 因子组合的样本量，出现 0 的地方就是问题所在 table(dm$cut, dm$color) # 过滤后重新降水平，防止空水平残留 dm$cut <- droplevels(dm$cut) dm$color <- droplevels(dm$color)

参数说明：droplevels()只删掉样本量为 0 的水平，不影响其他水平。table()输出里出现 0 的格子，要么合并水平，要么在模型里去掉这个交互项，没有第三种办法。

5.5 拿训练集 R² 当预测精度

现象：模型 R² = 0.98，看着很美，拿几颗新钻石的报价去预测，误差两三千美元。

原因：R² 是样本内指标，变量越多它越高，跟预测能力不是一回事。尤其是这个数据集里 carat 和 x、y、z 高度共线，模型相当于用四个变量反复描述同一件事，样本内拟合自然好。另外异常值对 R² 的影响很大，几个极端点就能把它拉上去或拽下来。

解决：划分训练集和测试集，在测试集上算 RMSE 和 MAE，这两个才是能拿出去说的数。同时看残差的分布，如果残差在不同价格区间上系统性偏正或偏负，说明模型有结构性偏差，加变量也救不回来。具体做法放下一章。

6. 模型验证与一个可复用技巧：留出集 RMSE 和 log 反变换的修正

模型跑完到能交付，中间还差两步：在没见过的数据上量一次误差，以及把对数模型的预测值正确换回美元。第二步最容易被忽略，也最容易造成系统性低估。

6.1 留出集验证：RMSE 怎么算才有意义

set.seed(2024) n <- nrow(dm) idx <- sample.int(n, size = floor(0.7 * n)) train <- dm[idx, ] test <- dm[-idx, ] fit_tr <- lm(log(price) ~ log(carat) + cut + color + clarity, data = train) pred_log <- predict(fit_tr, newdata = test) rmse <- function(actual, pred) sqrt(mean((actual - pred)^2)) mae <- function(actual, pred) mean(abs(actual - pred))

逻辑说明：sample.int()抽行号而不是直接抽数据框，后面用-idx取补集时不会出错。70/30 是个常规比例，5 万多条记录里拿 1.6 万条做测试足够稳定。predict()的newdata必须包含模型里用到的所有变量，缺一个就报错——这是好事，能提前发现测试集列名不一致的问题。

参数说明：set.seed(2024)固定划分，换机器重跑得到同一个测试集，否则两次结果没法比。RMSE 的单位和响应变量一致，对 log 模型来说是「对数价格」的单位，不能直接读成美元，得先反变换。

6.2 log 反变换：直接 exp 会系统性低估

朴素做法是exp(pred_log)，但这一步有问题。对数变换是非线性的，E[exp(Y)] ≠ exp(E[Y])，直接取指数得到的是中位数预测而不是均值预测，整体偏低。修正办法叫 smearing estimate，用训练集残差构造一个缩放因子：

# 朴素反变换 pred_naive <- exp(pred_log) # smearing 修正因子：训练集残差取 exp 后的均值 smear <- mean(exp(residuals(fit_tr))) pred_smear <- pred_naive * smear # 两种反变换在测试集上的误差对比 rbind( naive = c(RMSE = rmse(test$price, pred_naive), MAE = mae(test$price, pred_naive)), smear = c(RMSE = rmse(test$price, pred_smear), MAE = mae(test$price, pred_smear)) )

逻辑说明：residuals(fit_tr)是训练集上的对数残差，exp()之后取均值就是那个缩放因子。它衡量的是「对数尺度上的误差在原始尺度上会被放大多少」，价格分布越分散，这个因子离 1 越远。

参数说明：smear必须用训练集残差算，用测试集算等于偷看了答案。这个因子通常在 1.05 到 1.3 之间，价格右偏越厉害越大。如果算出来接近 1，说明对数变换没什么必要，直接建线性模型更省事。

对比表我一般这么看：

6.3 从那以后我每次都先跑一遍反变换对比

吃了两次亏之后我养成一个习惯：任何带 log 变换的模型，在写结论之前先把 naive 和 smear 两列误差并排打出来。第一次踩坑是给一个采购场景做报价预估，模型 R² 很好看，交付之后对方反馈「你们的

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