★ FEATURED ARTICLE 简介围绕R语言钻石价格分析的学习文档面向希望掌握数据清洗、可视化与回归建模的初学者及数据分析爱好者借助ggplot2自带、约5.4万条记录的diamonds数据集完整走通从数据概览、缺失值与重复值检查、异常值识别剔除到价格分布直方图、箱线图和小提琴图探索再到简单线性回归与多元回归建模、价格预测的分析链路。文档重点讨论克拉重量、切工、颜色、净度等变量对价格的影响呈现价格分布的右偏特征与千美元附近峰值并通过逐步筛选梳理主要影响因素兼顾方法理解与结果解读。资源包仅含1个docx文件约2.4MB篇幅集中便于按分析步骤阅读和复盘适合课堂示例、课程作业或自学练手时参考。已有284人学习可作为R语言回归分析与机器学习入门实践的配套材料。1. 反直觉开场Ideal 切工的钻石价格为何不占顶先抛一个我拆这份 R 语言钻石价格分析时没预料到的结果把 ggplot2 自带 diamonds 数据集的 price 画成小提琴图切工等级标到 Ideal理想的那一组价格中位数并没有甩开 Premium甚至分布下沿还压着一大片低价钻石。这跟切工越好越贵的常识是拧着的。这份基于 R 语言的数据分析项目做的就是把这层拧巴拆开——用 53940 条钻石记录、10 个变量从缺失值清洗、异常值剔除、描述性可视化一路走到多元线性回归、交互项筛选和回归诊断最后回答一个问题到底哪些因素在真正推高钻石价格。适合手里有 R 环境、想练一遍完整回归分析链路的人也适合拿它当课程设计或数据岗面试作品底稿的开发者。2. 数据底座与工程骨架54000 行 diamonds 怎么摆进一个能跑的 R 项目diamonds 是 ggplot2 包自带的数据集省掉了找数据、爬数据、对齐字段这些前置脏活。但数据白送不等于项目能跑很多人第一次打开 RStudio 就是一行data(diamonds)然后在一个 .R 文件里从头写到尾写到第 300 行想改个异常值阈值发现前面所有图都得重跑一遍。我习惯先把脚本按职责切开再谈代码。2.1 目录结构把一次性脚本拆成可重跑的流水线一个能反复迭代的分析项目目录不需要复杂但必须让清洗和建模物理隔离。清洗规则一改只需要重跑01_clean.R模型脚本读的是落盘的中间结果不用连带重画十几张图。diamond-price-analysis/ ├── R/ │ ├── 00_setup.R # 依赖包安装与全局参数对比度、随机种子 │ ├── 01_clean.R # 结构检查、缺失值、重复值、异常值、特征构造 │ ├── 02_eda.R # 直方图、核密度、小提琴图、箱线图 │ ├── 03_model.R # 一元回归、多元回归、交互项、逐步回归 │ └── 04_diagnose.R # 残差图、QQ 图、VIF、显著性检验 ├── data/ │ └── diamonds_clean.rds # 清洗后的中间结果建模脚本只认这个 ├── output/ │ ├── fig_price_hist.png │ └── model_summary.txt └── report.Rmd # 汇总成可 knit 的报告data/diamonds_clean.rds是整个结构的枢纽。用saveRDS()而不是write.csv()是因为因子变量的 levels 顺序会随 CSV 读写丢失而 levels 顺序直接决定回归里谁是基准组——这个坑后面单独讲。2.2 十个变量的类型判定哪些能直接进模型哪些要先转diamonds 的 10 个变量里只有 price 是目标变量其余 9 个是特征。R 读进来时 cut、color、clarity 是 character 类型不转 factor 直接丢进lm()模型会把它当连续变量算系数含义全错。变量名含义取值范围原始类型建模前处理carat克拉重量0.2 - 5.01numeric直接入模cut切工质量Fair / Good / Very Good / Premium / Idealcharacter转有序因子color颜色等级J最差到 D最好character转有序因子注意方向clarity净度I1 到 IF最差到最好character转有序因子depth总深度百分比43 - 79numeric直接入模table台面宽度比43 - 95numeric直接入模price美元价格326 - 18823numeric目标变量常取对数x / y / z长 / 宽 / 深毫米x: 0-10.74, y: 0-58.9, z: 0-31.8numeric需先处理 0 值再入模表里 color 这一行要特别留心J 是最差、D 是最好的但字符排序里 D 排在 J 前面。如果直接as.factor()R 会按字母序把 D 设成基准组回归系数就变成相对最好的颜色其他颜色便宜多少方向恰好反了。这是这类项目里最常见的解释性翻车点。2.3 依赖与初始化四行代码定全局# 00_setup.R pkgs - c(ggplot2, dplyr, mice, car, lmtest) for (p in pkgs) if (!requireNamespace(p, quietly TRUE)) install.packages(p) library(ggplot2) library(dplyr) # 关键把无序对比度设为 treatment保证因子系数是相对基准组的解释 options(contrasts c(contr.treatment, contr.poly)) set.seed(2024) # 读数据diamonds 随 ggplot2 一起加载不需要外部文件 data(diamonds, package ggplot2)options(contrasts ...)这一行是整份脚本里性价比最高的一行。R 对有序因子ordered factor默认用contr.poly多项式对比输出的是线性、二次、三次趋势项系数名叫.L、.Q、.C业务方看了完全不知道在说什么。改成contr.treatment后系数就是该水平相对基准组的价格差异能直接写进报告。set.seed(2024)在这个项目里主要给后续的抽样检验用——正态性检验对全量 53940 行跑不动必须抽样抽样就得固定种子否则每次 knit 报告结论都在抖。3. 清洗链路缺失值、重复行与箱线图异常值的处理顺序清洗顺序是有讲究的先看结构再查缺失再查重复最后处理异常值。顺序反了会出问题——比如先去重再查缺失你会误以为数据本来就没有缺失其实缺失行已经被去重顺手删掉了。3.1 结构、缺失值与重复行的三连检查# 01_clean.R str(diamonds) # 看每个变量的类型和前几个取值 summary(diamonds) # 看数值变量的五数概括和因子水平计数 head(diamonds, 6) # 前六行 tail(diamonds, 6) # 后六行 colSums(is.na(diamonds)) # 每列缺失值个数 sum(duplicated(diamonds)) # 完全重复的行数str()输出里能看到num和chr两种类型标注这是判断要不要转因子的第一手依据。summary()更有信息量如果某个数值变量的 Max 明显偏离中位数几个数量级基本就能锁定异常值区间diamonds 里 y 的最大值 58.9 毫米对一颗最大 5 克拉的钻石来说显然不合理。colSums(is.na())返回全 0说明这个数据集没有缺失值原文里用mice::md.pattern()画缺失模式图也是同样的结论。没有缺失值就不需要插补但这一步不能省——真实业务数据里缺失率 30% 以上是常态养成先查的习惯比记住这个数据集没缺失有用得多。sum(duplicated())返回 0说明没有完全重复行。注意这里是完全重复——同一颗钻石被录入两次、但价格写错了duplicated()是发现不了的那是业务主键层面的去重跟统计意义上的重复行不是一回事。3.2 箱线图异常值剔还是不剔先想清楚# 画价格箱线图肉眼确认异常点规模 boxplot(diamonds$price, main 钻石价格箱线图, ylab price (USD)) # 用 IQR 规则定位异常区间 q1 - quantile(diamonds$price, 0.25) q3 - quantile(diamonds$price, 0.75) iqr - q3 - q1 lower - q1 - 1.5 * iqr upper - q3 1.5 * iqr cat(下界:, lower, 上界:, upper, \n) # 剔除并落盘 diamonds_clean - diamonds %% filter(price lower price upper) saveRDS(diamonds_clean, data/diamonds_clean.rds)箱线图跑出来你会看到中位数紧贴下四分位数上须以外拖出一长串点——这就是原文说的存在偏离三个标准差之外的点。IQR 规则1.5 倍四分位距和 3σ 规则在这里结论不一样价格分布严重右偏3σ 的边界会被少数极高价拉得很宽剔掉的点反而少IQR 更严格。我一般先用 IQR 看规模再决定是用 1.5 倍还是 3 倍。这里必须提醒一句在价格分析里剔异常值是个需要犹豫的动作。价格 18823 美元那颗钻石不是录入错误它是一颗真实的大克拉高净度钻石。按 IQR 一刀切下去高端样本几乎被清空模型学到的是价格上限大概 5000 美元拿去预测大钻会系统性低估。原文的做法是剔了课程作业里没问题落到真实业务上我一般改成winsorize 缩尾或者只在诊断阶段标记、不参与训练。3.3 变量类型转换与衍生特征构造diamonds_clean - diamonds_clean %% mutate( # 有序因子levels 从差到好基准组自动落在第一个 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), # 衍生特征体积近似、克拉单价、深度比值 volume x * y * z, price_per_carat price / carat, depth_ratio 2 * z / (x y) )levels的顺序决定了回归输出里谁是基准Intercept 吸收的那一组。cut 的基准是 Faircolor 的基准是 Jclarity 的基准是 I1——都是最差档这样每个系数都是比最差档贵多少正负号一看就懂。三个衍生变量各有用途volume用 x·y·z 近似体积比单看 carat 更能反映看起来大不大price_per_carat是克拉单价做分档对比时比总价更公平depth_ratio用 2z/(xy) 复算深度百分比可以跟原始 depth 列做一致性校验差得多的行往往是尺寸录错了。注意volume x * y * z这一行有隐患diamonds 里存在 x 或 y 或 z 为 0 的记录乘出来 volume 就是 0后面取对数直接变 -Inf。处理办法在避坑那章展开。3.4 价格分布直方图、核密度与条形图# 直方图 核密度曲线 正态参考线 ggplot(diamonds_clean, aes(x price)) geom_histogram(aes(y after_stat(density)), bins 50, fill steelblue, color white, alpha 0.7) geom_density(color darkred, linewidth 1) stat_function(fun dnorm, args list(mean mean(diamonds_clean$price), sd sd(diamonds_clean$price)), linetype dashed) labs(title 钻石价格分布, x 价格美元, y 密度) theme_minimal()after_stat(density)是把直方图的纵轴从计数换成密度否则直方图和核密度曲线的量纲对不上画在一起没法看。bins 50是按 53940 条样本量估的样本越多分箱可以越细。跑出来三个特征很明显整体右偏、1000 美元附近有个尖峰、2000 / 4000 / 8000 美元附近出现断层。右偏符合奢侈品定价规律——低价钻石数量多高价钻石稀缺。1000 美元那个峰是市场心理价位聚集price_per_carat在这个区间成交最密集。stat_function(fun dnorm, ...)叠的正态参考线在 8000 美元之前跟核密度曲线几乎不重合之后两条线才贴近——这说明价格不能用正态分布假设去建模后面做显著性检验时选非参数方法根子就在这张图上。3.5 分类型变量与价格的关系图# 小提琴图切工 vs 价格价格取 log 便于观察低值区 ggplot(diamonds_clean, aes(x cut, y price, fill cut)) geom_violin(trim FALSE, alpha 0.6) scale_y_log10() labs(title 不同切工等级的价格分布, x 切工, y 价格log10) theme_minimal() theme(legend.position none) # 条形图各切工等级的数量 ggplot(diamonds_clean, aes(x cut, fill cut)) geom_bar() geom_text(stat count, aes(label after_stat(count)), vjust -0.5) theme_minimal() theme(legend.position none)scale_y_log10()是关键。价格跨度从 326 到 18823线性轴上低价区全被压成一条线小提琴的形体根本看不出来取对数后分布形状才展开。原文说此处价格取开值log便于观察就是这个意思。小提琴图给出的结论跟条形图要合起来读Ideal 切工数量 21488 颗、Premium 13748 颗、Very Good 12069 颗、Good 4891 颗、Fair 1598 颗Ideal Premium 占比过半。数量最多的 Ideal 组价格分布最宽、下沿最低说明好切工在这个数据集里不是稀缺品而是标配——正因为供给充足它反而不构成价格溢价的主要来源。这个结论要等多元回归把 carat 控制住之后才能确认。color 和 clarity 的箱线图更规整颜色从 J 到 D、净度从 I1 到 IF价格中位数基本单调上升但 clarity 的 I1 组中位数反而比 SI2 高——样本量小加上个别大钻拉高了中位数这种小样本组的图不要过度解读。4. 建模主线从一元回归到带交互项的多元回归与回归诊断清洗完的数据落在data/diamonds_clean.rds建模脚本只读这个文件。整条建模主线是递进的先钉死单个最强特征再堆全变量看 R² 天花板再加交互项抠剩余方差最后用诊断图判断模型假设成不成立。4.1 一元线性回归先把克拉和价格的关系钉死# 03_model.R diamonds_clean - readRDS(data/diamonds_clean.rds) fit1 - lm(price ~ carat, data diamonds_clean) summary(fit1) confint(fit1, level 0.95) # 系数的 95% 置信区间summary()输出里 carat 的系数约为 7760p 值远小于 0.001含义是克拉重量每增加 1 克拉钻石价格平均上涨约 7760 美元在 95% 置信区间内。这是全项目里解释力最强的一个单变量——Multiple R-squared大概在 0.85 上下一颗克拉数就吃掉了大部分方差。但一元回归的残差图会暴露问题残差随 fitted value 增大呈喇叭口张开典型的异方差。原因很直白——1 克拉的钻石价格波动几千美元3 克拉的钻石波动几万美元误差尺度本身就跟价格量级绑在一起。这也是后面要给 price 取对数的动机。confint()这行别省。只报点估计在业务沟通里是耍流氓区间宽窄直接反映这个系数可不可信宽度超过系数本身一半的时候我会回去检查样本量或者共线性。4.2 哑变量与多元回归R²0.88 之后怎么办fit2 - lm(price ~ carat cut color clarity depth table x y z, data diamonds_clean) summary(fit2)因子变量进lm()时 R 自动展开成哑变量7 个水平的 clarity 展开成 7 列基准组 I1 被 Intercept 吸收。这一步不需要手动model.matrix()手动展开反而容易漏掉基准组导致完全共线。fit2的R-squared大约 0.8798。从 0.85 到 0.88加了 8 个变量只换来 3 个百分点说明两个事一是 carat 已经吃掉了主要信号二是剩余方差里可能藏着非线性结构线性加变量加不出来。这时候最容易被误用的判断是R² 0.88 已经很高了模型可用。R² 高只说明拟合样本好不说明残差满足同方差、独立性假设。0.88 的模型如果残差图是喇叭形预测区间仍然不可信——这是后面 4.5 节要处理的事。另外注意 x、y、z 和 carat、volume 之间的共线性。carat 和 volume 的相关系数通常在 0.95 以上同时入模会让两个系数的标准误同时膨胀VIF 会直接报警。我在fit2里保留原始 x/y/z 而不放 volume就是为了避开这个。4.3 交互项与逐步回归AIC 怎么用# 构造带一阶交互项的候选模型 fit3 - lm(price ~ (carat cut color clarity)^2 depth table x y z, data diamonds_clean) # 向后逐步回归按 AIC 删变量 fit3_step - step(fit3, direction backward, trace 1) summary(fit3_step) AIC(fit2, fit3_step) # 两个模型横向比较(carat cut color clarity)^2是 R 的公式语法糖展开后等于主效应加上两两交互项不包含平方项。这里只让四个核心变量交互depth、table、x/y/z 保持主效应——全量交互会产生上百个项跑得慢还容易过拟合。direction backward是从全模型开始逐个删trace 1会把每一步删了哪个变量、AIC 降了多少打印出来。原文说所有变量都被选入了模型说明在这个数据规模下交互项确实带来了拟合提升AIC 没有给出删除信号。AIC(fit2, fit3_step)这行是横向对比的唯一正确姿势。AIC 只能比同一份数据、同一个因变量上的模型fit2和fit3_step满足条件可以比但如果把lm(price ~ ...)和lm(log(price) ~ ...)放一起比 AIC那是无效比较——因变量变了似然不可比。4.4 显著性检验非参数检验在这类数据上更稳# 不同切工类型之间价格是否有显著差异 kruskal.test(price ~ cut, data diamonds_clean) # 不同颜色、不同净度同理 kruskal.test(price ~ color, data diamonds_clean) kruskal.test(price ~ clarity, data diamonds_clean) # 正态性检验前先抽样全量会直接报错 set.seed(2024) sample_price - sample(diamonds_clean$price, 5000) shapiro.test(sample_price)选kruskal.test()而不是aov()理由是 3.4 节那张直方图已经证明了价格分布严重右偏方差分析的 normality 方差齐性两个前提都不满足。Kruskal-Wallis 只要求各组独立对分布形态没要求在这类偏态数据上是更稳的选择。三个检验的 p 值都远小于 0.01结论一致不同切工、不同颜色、不同净度组之间的价格差异是统计显著的。但显著不等于重要——样本量五万多的时候任何微小差异都会被检验成显著效应量比如各组中位数差才是判断业务意义的东西这一点在写报告时要说清楚。shapiro.test()那一行的 5000 不是随便写的这个函数的硬限制是样本量不能超过 5000超了直接抛错。原文里对全量分组做正态性检验如果没抽样是跑不通的。4.5 回归诊断残差图、QQ 图与 VIF# 04_diagnose.R par(mfrow c(2, 2)) plot(fit3_step) # 多重共线性 car::vif(fit3_step) # 残差自相关按行号顺序的 Durbin-Watson lmtest::dwtest(fit3_step)四张诊断图各有分工Residuals vs Fitted 看异方差——点如果呈喇叭形张开说明方差随拟合值变化Normal Q-Q 看残差正态性——点偏离对角线两端翘起说明尾部厚Scale-Location 是标准化残差的开方比第一张更容易看出趋势Residuals vs Leverage 找强影响点——Cooks distance 超过 1 的行要单独看。vif()输出里如果 x、y、z 三个的 VIF 都在 5 以上说明尺寸变量之间高度共线处理办法是只留一个通常留 volume 或者只留 carat。dwtest()的 p 值在这个数据集上通常不显著因为 diamonds 的行序跟价格没有系统性关联不构成时间自相关——但这个检验在真实业务数据里必做一旦数据是按时间导出的残差自相关会让标准误被严重低估。诊断结论直接决定下一步异方差显著就给 price 加 log 变换重跑fit2共线性显著就砍变量强影响点过多就回去看 3.2 的异常值处理是不是剔得太狠。5. 避坑排查这份钻石回归脚本最容易翻车的 5 个地方下面这 5 条都是我在跑这份分析时真真切切卡住过的按现象 → 原因 → 解决记下来照着对能省不少时间。5.1 现象shapiro.test 一跑就报错 sample size must be between 3 and 5000原因shapiro.test()内部实现基于 Shapiro-Wilk 统计量样本量硬上限 5000diamonds 全量 53940 行必然超限。解决按组抽样每组控制在 5000 以内再检验同时固定set.seed()保证可复现set.seed(2024) diamonds_clean %% group_by(cut) %% slice_sample(n 5000, replace FALSE) %% group_modify(~ broom::tidy(shapiro.test(.x$price)))如果抽样后仍然拒绝正态不要硬转参数检验直接用 4.4 节的kruskal.test()走非参数路线。5.2 现象cut 的回归系数输出成.L、.Q、.C方向完全看不懂原因cut被定义成ordered TRUE的有序因子R 对有序因子默认使用contr.poly多项式对比系数表示的是线性趋势、二次趋势、三次趋势不是相对某个基准组的差异。解决两种改法按需选一种。要趋势解释就保留多项式对比但必须在报告里说明.L代表切工等级每升一档的线性效应要基准组差异解释就在00_setup.R里加options(contrasts c(contr.treatment, contr.poly))或者对单个变量用C(cut, contr.treatment)显式指定。改完之后记得重跑一遍summary()基准组变了所有系数的含义都会跟着变。5.3 现象异常值剔完之后模型对高价钻石的预测系统性偏低原因3.2 节用 IQR 规则按 price 剔异常价格分布右偏严重上界通常落在 5000 美元上下19000 美元的大钻全被当成离群点删了。模型训练集里最大价格只剩 5000自然学不会预测 15000。解决分两种情况处理。如果分析目标是解释价格影响因素保留全部样本异常值只在诊断阶段用plot(fit3_step)的第 4 张图标记出来不参与删除如果确实是录入错误比如 y 58.9 毫米那种物理上不可能的尺寸按行删删完要记录删了多少行、删的比例是多少。# 只按物理尺寸合理性剔不按价格剔 diamonds_clean - diamonds %% filter(y 0, y 20, x 0, z 0)这样剔掉的是明显录错的尺寸高价大钻保留在训练集里模型才有机会学到价格上限。5.4 现象volume 算出 0 和 NaN取 log 之后整列变成 -Inf原因diamonds 里有 x 或 y 或 z 为 0 的记录三个尺寸相乘得 0log(0)是 -Inf如果后面还有除法比如price / volume会出现0/0 NaN。这些值一旦进模型lm()会整行丢弃并抛 warning但很多人不看 warning最后发现样本量对不上。解决构造衍生特征前先处理 0 值要么过滤如上一条要么加极小值兜底diamonds_clean - diamonds_clean %% mutate( volume pmax(x * y * z, 1e-6), # 兜底避免 log(0) price_per_carat price / pmax(carat, 1e-6) )pmax()比ifelse()快而且不用写两遍条件。加完兜底之后再跑一遍summary(diamonds_clean$volume)确认最小值不再是 0。5.5 现象step() 跑了很久没出结果或者结果每次不一样原因候选模型变量太多全量交互项轻松上百个逐步回归每一步都要重算 AIC耗时随变量数是乘法级增长另一个原因是没设set.seed()如果候选模型里存在 AIC 相同的平局step()的选择顺序不保证稳定。解决先降维再逐步。把明显共线的变量x/y/z 只留一个和业务上无关的变量比如 depth 和 table 对价格解释力本来就弱先手工剔除候选模型控制在 30 个变量以内再跑step()。同时固定种子并把trace 1打开每一步删了什么、AIC 变化多少都能看到出问题能定位到具体哪一步。6. 二次开发把一次性分析脚本改成可复用的建模流程这份脚本跑通一次之后最有价值的动作不是交报告是把它改成能换数据复用的流程。我一般从三个方向改。第一把建模逻辑封成函数入参是数据框和公式出参是拟合对象加指标表。这样换一份带同样字段的钻石数据改一行调用就行不用复制粘贴整段代码build_price_model - function(df, formula price ~ carat cut color clarity) { fit - lm(formula, data df) tibble( r2 summary(fit)$r.squared, adj_r2 summary(fit)$adj.r.squared, aic AIC(fit), n nobs(fit) ) } # 对比三种变量组合不用写三遍 lm map_dfr( list(price ~ carat, price ~ carat cut color clarity, price ~ (carat cut color clarity)^2), ~ build_price_model(diamonds_clean, .x), .id model )第二补一层留出验证。现在所有 R² 都是样本内拟合值加个 5 折交叉验证看样本外误差才能判断交互项是不是过拟合了改造点当前状态建议改法验证方式模型评估只看样本内 R²加 5 折 CVcaret::train()看 RMSE 波动因变量原始 pricelog(price)对比残差图喇叭口是否收敛共线性靠 VIF 事后发现建模型前先算相关矩阵相关系数 0.9 的变量二选一分类型变量手动设 levels封装成统一转换函数新增变量时不用重写第三把report.Rmd里的数字全部改成引用变量不要手写。报告里写R² 0.8798这种硬编码下次数据一更新就全错。用 inline R 代码r round(summary(fit2)$r.squared, 4)让数字跟着模型走knit 一次自动刷新。改完这三处这个项目就从一份课程作业变成了一个能接新数据的分析模板。我踩过最贵的一次坑是在一个类似的价格预测项目里样本内 R² 做到 0.92 就交付了上线后样本外误差翻了三倍——因为交互项把训练集的噪声也拟合进去了。从那以后我每次交付前都强制走一遍交叉验证样本内指标再漂亮也要看一眼 CV 的 RMSE 标准差。希望这份拆解能帮到你少走一遍我走过的弯路。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
