★ FEATURED ARTICLE 简介围绕钻石价格预测的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一堆用不上的东西包承担职责什么时候必须用ggplot2数据源 全套图形语法数据加载和所有可视化dplyr过滤、变形、分组汇总清洗和派生变量tidyr长宽表转换、缺失值整理分组统计后整理成图用格式micemd.pattern缺失值模式可视化缺失值排查阶段carVIF 共线性诊断、方差分析模型诊断阶段MASSstepAIC逐步回归变量筛选阶段broomtidy/glance把模型结果转数据框需要批量对比多个模型时机器学习里最朴素的回归分析其实就靠这几块拼起来。真正容易翻车的地方不是包不够是因子水平顺序、异常值判定口径这些细节后面单独用一章讲。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())返回全 0md.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)系数约 7760p 值远小于 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参数设成了 2AIC对五万行的大样本来说惩罚太轻容易把噪声变量留在模型里另一个常见原因是交互项写得太宽(.)^2一次性生成上百个交互逐步筛选根本扛不住。解决样本量大时改用 BICk log(nrow(train))交互项只挑业务上说得通的三个carat 和三个品质变量其余不生成。变量多不等于模型好这是我在这个数据集上吃过的最实在的教训。5.4 Ideal 切工均价反而低被质疑结论写反了现象分组统计出来 Ideal 切工的平均价格低于 Premium和切工越好越贵的常识冲突评审直接质疑分析有误。原因这是典型的混杂效应。carat 和 cut 之间存在相关性Ideal 组里集中了大量小克拉钻石均价被拉低单独看 cut 的边际效应是负的控制 carat 之后系数符号会翻过来。解决所有关于单变量的结论必须先做分组控制或者直接读多元回归里的偏效应系数。图表标题里写清未控制克拉重量能省掉很多解释成本。5.5 残差 QQ 图两端翘R² 怎么调都上不去现象模型 R² 卡在 0.88QQ 图两端明显偏离对角线预测高价钻石时系统性偏低。原因price 是右偏分布线性模型假设残差正态右尾的长尾数据把拟合线往下拽。解决对 price 取自然对数后建模预测时再用exp()还原。副作用是模型变成预测价格的对数误差解释要改成相对误差百分比汇报口径得同步调整别拿 RMSE 直接和原尺度模型比。6. 二次开发改造把一次性脚本升级成可复用的定价组件把这份分析从交作业改造成能给别人用我走过几条路按投入产出比排一下改造方向具体做法收益工作量函数化封装把清洗和预测各包成一个函数换数据集只改参数半天交叉验证caret或tidymodels做 5 折 CV指标更可信半天模型对比随机森林、GBM 与线性模型同台找到非线性上限一天交互查询Shiny 输入克拉/切工实时出价演示效果最好一天报告自动化Quarto 参数化渲染每月数据更新一键出报告半天最值得先做的是函数化。把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()看量纲和极值、箱线图看离群点、残差图看模型假设三张图都没问题才允许动模型。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
提前堵住了它)
