★ FEATURED ARTICLE 简介面向希望用R语言入门并实践机器学习的数据分析学习者这份资料围绕监督学习中的回归任务展开覆盖简单线性回归、多元线性回归、多项式回归、决策树回归与支持向量回归等模型并配有数据预处理模板帮助读者理解从数据清洗、特征处理到模型训练与比较的基本流程。压缩包共20个文件以8个R脚本和6个CSV示例数据为主另含少量系统缓存及授权说明整体约16KB体量轻便便于按回归方法分目录查阅和运行。已有245人学习适合作为课程练习、项目起步或算法复现的参考。目录按回归模型与数据预处理模块组织脚本与数据对应读者可据此快速搭建实验、对照不同模型的实现差异并在此基础上扩展调参与评估。1. 一个公式加一个函数R 的机器学习凭什么这么“短”在 R 里做机器学习一行lm(y ~ ., data d)就能拟合出一个模型换成别的技术栈从定义占位符到组装计算图代码可能要多出好几倍。很多人把它当成 R 语法糖吃得好但真正让这套写法能同时服务线性回归、广义线性模型、决策树和随机森林的是一套叫 formula 加设计矩阵的中间表示。理解这道中间层就能明白 R 的机器学习算法为什么能“换模型不换写法”也能知道它在极限场景下会从哪里漏出破绽。2.1 terms 对象符号表达式怎样变成一张可执行的计算计划y ~ x * g这行代码在 R 里并不会立刻计算任何东西它首先被解析成一个call对象也就是一段还没被求值的函数调用表达式。真正把它变成计算计划的是terms()它把公式拆成一个带属性的terms对象记录这次建模要哪些变量、它们以什么方式进入、有没有截距、交互项到几阶。这个对象的属性不是装饰而是所有下游函数共同遵守的接口约定。d - data.frame( y c(2.1, 3.4, 3.9, 5.2, 6.1, 7.8), x c(0, 1, 2, 3, 4, 5), g factor(c(a, b, a, b, a, b)) ) tt - terms(y ~ x * g, data d) attr(tt, term.labels) # x g x:g 要进入模型的项 attr(tt, order) # 1 1 2 每项的交互阶数 attr(tt, intercept) # 1 是否含截距 attr(tt, factors) # 变量 x 项 的指示矩阵 # x g x:g # x 1 0 2 # g 0 1 2factors这个矩阵是整套机制的关键行是原始变量列是模型项数字1表示该变量直接构成该项2表示它参与了一个交互项。有了它任何建模函数都不用再去解析公式字符串只需要读这一个矩阵就能知道“第 j 个系数对应哪几个变量的乘积”。variables属性则保存了每个变量的调用表达式所以像log(x)、poly(x, 2)这种写在公式里的变换也能被正确求值。这也解释了一个常见疑问为什么新写的算法只要提供一个接收设计矩阵的fit函数就能直接吃y ~ .的语法因为公式解析这一步已经被terms和model.frame消化掉了。2.2 model.matrix 的编码规则因子、交互项和截距到底怎么排布terms只描述“要什么”model.matrix负责“怎么摆”。它先调用model.frame在数据环境里求值顺带处理subset、na.action、drop.unused.levels然后按对比编码contrasts把因子展开成数值列。默认行为由options(contrasts)决定无序因子用contr.treatment也就是 k 个水平只产生 k-1 列第一个水平被截距吸收成基准有序因子用contr.poly展开成正交多项式列。X - model.matrix(tt, data d) colnames(X) # (Intercept) x gb x:gb dim(X) # 6 4 # 去掉截距后第一个因子会展开成全部 k 列 colnames(model.matrix(~ x g - 1, data d)) # x ga gb这里有个容易被忽略的细节即便是-1去掉截距也只有第一个因子会展开成 k 列后面的因子仍然只用 k-1 列。原因是全部展开会让设计矩阵列和等于截距列直接造成秩亏。model.matrix用这种方式在“每个水平都有独立参数”和“矩阵满秩”之间找平衡。交互项x:gb则是把x列和gb哑变量列逐元素相乘算出来的是一列新特征而不是某种特殊的运算。秩亏在lm里不会报错而是被安静地处理掉d$x2 - d$x * 2 # 与 x 完全共线 coef(lm(y ~ x x2, data d)) # (Intercept) x x2 # 2.0 1.0 NANA不是拟合失败而是 QR 分解在选主元时判定这一列与前面对应列线性相关把它踢出了独立列集合。判定阈值tol 1e-7是相对于原始列范数的不是绝对值所以量纲变化不会轻易改变判定结果。2.3 QR 分解为什么 lm 不去求 (XX) 的逆教科书上的最小二乘解是(XX)^{-1} Xy但 R 的lm.fit走的是另一条路底层 C 代码Cdqrls对设计矩阵做 Householder QR 分解解R b Qy。原因是条件数会被平方。如果X的条件数是 10^8XX的条件数就是 10^16已经超过双精度浮点约 10^16 的有效位数正规方程解出来的系数可能连符号都是错的。q - qr(X) q$rank # 独立列数 q$pivot # 列的重排顺序前 rank 个是独立的 R - qr.R(q) # 上三角因子 all.equal(crossprod(R), crossprod(X)) # TRUE仅限无主元交换时summary.lm计算标准误时用的也是这个Rchol2inv(R)得到(XX)^{-1}再乘以残差方差sum(resid^2) / (n - rank)。这里的分母用的是矩阵的秩而不是列数所以刚才那个共线的例子里自由度是 n-2 而不是 n-3。整套推断链条从系数、标准误到 F 检验全部建立在同一次 QR 分解的中间结果上不需要重复计算。提示判断共线性严重程度看qr(X)$rank和1 / rcond(X)不要凭肉眼看相关系数两个变量相关系数 0.6 也可能因为第三个变量的存在而联合共线。3.1 自己写一遍逻辑回归的负对数似然与解析梯度树模型和线性模型之外R 里大部分“训练”动作都能归结为一个优化问题给定损失函数找一组参数让损失最小。逻辑回归是最适合拿来拆的例子因为它的梯度形式极简X(p - y)一眼就能看出每一步在把预测概率往标签方向拉。# 负对数似然数值稳定写法避免 exp(-z) 溢出成 Inf nll - function(b, X, y) { z - drop(X %*% b) # log(1 exp(-z)) 改写成 pmax(z,0) log1p(exp(-abs(z))) sum(pmax(z, 0) - y * z log1p(exp(-abs(z)))) } # 解析梯度p 是当前参数下的预测概率 grad - function(b, X, y) { p - plogis(drop(X %*% b)) drop(crossprod(X, p - y)) # X(p - y) }那个pmax(z, 0) log1p(exp(-abs(z)))不是炫技。当z是 -800 时exp(800)直接是Inf整个损失变成Inf优化器会立刻停在原地而log1p(exp(-800))是 0结果是精确的 800。很多人自己写的损失函数在中等规模数据上跑得好好的数据一放大就返回NaN罪魁祸首几乎都是这里。3.2 optim 与 glm两种接管方式IRLS 是隐藏的主角R 提供两条优化路径。一条是通用的optim你给损失函数和梯度它用 BFGS、L-BFGS-B 或 Nelder-Mead 去找最小值另一条是glm内部专门实现的迭代重加权最小二乘IRLS它利用了指数族的结构把非线性问题转成一串加权线性最小二乘每一轮都是一次 QR 分解。set.seed(1) n - 2000; p - 5 X - cbind(1, matrix(rnorm(n * p), n, p)) beta - c(-0.5, 1.2, -0.8, 0.4, 0, 0.9) y - rbinom(n, 1, plogis(drop(X %*% beta))) # 平均梯度配合固定步长收敛更稳 b - rep(0, ncol(X)); lr - 0.5 for (i in 1:400) b - b - lr * grad(b, X, y) / n bIRLS 的一轮可以写成先算线性预测eta Xb再算均值mu linkinv(eta)和方差函数V(mu)构造工作响应z eta (y - mu) / mu.eta(eta)和工作权重w mu.eta(eta)^2 / V(mu)最后解一个以sqrt(w)加权的线性最小二乘。收敛判据是偏差的相对变化|dev - dev_old| / (|dev| 0.1) 1e-8默认最多 25 轮。这个结构解释了为什么glm比手写梯度下降快一个量级它每一步都在用二阶信息而梯度下降只用了一阶。optim的默认行为则是有限差分。不传gr时它对每个参数做前向差分步长由ndeps控制默认1e-3。参数有 100 个每轮就要多算 100 次函数值参数上千时这个开销会让本来能跑完的模型直接跑不完。3.3 收敛判据、步长和标准化loss 曲线抖动的三个来源第一是步长。求和梯度和平均梯度差一个n倍学习率自然也要差n倍拿对小样本调好的步长去跑大样本loss 一定发散。第二是特征尺度。年龄是几十收入是几万设计矩阵的条件数被拉到 10^5 以上梯度下降的等高线变成一条细长峡谷步长稍大就来回震荡稍小就爬不动。标准化之后再优化本质上是在降低条件数不是“让数据更好看”。第三是收敛判据只看相对变化。optim的reltol默认是1e-8当目标函数本身量级很大时绝对变化可能还很大就已经满足相对判据了反过来量级很小的时候又会一直迭代到maxit。判断是否真收敛看返回的convergence字段0才代表成功1是到达迭代上限10是 Nelder-Mead 的退化。注意优化前把特征标准化、把目标函数写成平均损失、给optim传解析梯度这三件事能把绝大部分“训练不收敛”的问题挡在门外。4.1 分裂点搜索Gini 增益在 C 代码里怎么被枚举决策树的训练过程比线性模型更依赖高效的底层实现因为要枚举的分裂点数量是变量数乘以每个变量的不同取值数。rpart把这部分放在 C 里流程是对每个预测变量把样本按该变量排序只在相邻的不同取值中点处试探分裂然后计算左右子节点的加权不纯度。分类树默认用 Gini1 - sum(p_k^2)回归树用误差平方和的下降量。gini - function(tab) { p - tab / sum(tab); 1 - sum(p^2) } # 一次分裂的加权不纯度 split_gini - function(y, left) { n - length(y); nl - sum(left) (nl / n) * gini(table(y[left])) ((n - nl) / n) * gini(table(y[!left])) }这个过程有两个隐含代价。一是它只做贪心搜索每一步选当前最优不保证全局最优所以树的结果对数据微小扰动很敏感。二是 Gini 天然偏好取值多的变量一个 ID 列能把每个样本分到单独一类加权不纯度降到 0但完全没有泛化能力。这个选择偏差不是实现 bug而是准则本身的性质条件推断树用置换检验来校正代价是计算量上去了。4.2 cp 与代价复杂度剪枝printcp 那张表怎么读rpart不靠预剪枝控制树的大小而是先长一棵大树再用代价复杂度剪枝往回砍。剪枝准则是最小化R(T) cp * |T|其中R(T)是树的误差|T|是叶子数cp是每增加一个叶子的“价格”。cp越大树越小。library(rpart) fit - rpart(Species ~ ., data iris, method class, control rpart.control(cp 0.001, minsplit 10, minbucket 5, xval 10)) printcp(fit)printcp输出的表里CP是每个候选剪枝点对应的复杂度参数nsplit是分裂次数rel error是训练集上的相对误差xerror是交叉验证误差xstd是它的标准误。选择通常看 1-SE 规则先找xerror最小的那一行再往上找第一棵满足xerror min(xerror) xstd的更小的树。这么做的逻辑是交叉验证误差本身有方差追求最小xerror很可能是在拟合验证集的噪声宁可要一棵小一点、误差在统计上 indistinguishable 的树。rpart.control里还有两个参数值得留意minsplit默认 20节点样本数不到就不尝试分裂minbucket默认是minsplit / 3限制叶子的最小样本量。这两个值在样本量小的数据集上设得太大会让树长不起来看起来像“模型没学到东西”实际是参数卡住了。4.3 bagging、mtry 与 OOB随机森林的随机性到底加在哪一层单棵树的问题是方差大随机森林的思路是用两组随机性把方差平均掉。第一层是 bagging每棵树用一个有放回抽样的自助样本训练样本量仍是 n大约 36.8% 的样本不会被抽中。第二层是特征随机每次分裂时不从全部 p 个变量里选而是先随机抽mtry个再挑最优。library(randomForest) set.seed(42) rf - randomForest(Species ~ ., data iris, ntree 500, mtry 2, importance TRUE, nodesize 1) rf$err.rate[500, ] # OOB 总体与各类错误率 importance(rf, type 2) # MeanDecreaseGinimtry的默认值随任务类型变分类是floor(sqrt(p))回归是floor(p/3)。这个差异的道理是分类树的分裂准则在低维子空间里更容易找到“够用”的分裂降低mtry能进一步降低树间相关性回归对单变量信息量更敏感抽得太少会让每棵树都太弱。nodesize也是分类默认 1、回归默认 5回归叶子强制多几个样本是为了压住单棵树在连续目标上的方差。未抽中的那部分样本顺带提供了免费的验证集每棵树用它的 OOB 样本做预测再汇总成 OOB 误差。这个估计几乎等价于一次 n 折交叉验证却不需要额外训练一遍模型是随机森林在中小数据上特别讨喜的原因。重要性指标有两种基于 Gini 下降的MeanDecreaseGini在训练集上算偏向高基数变量基于置换的MeanDecreaseAccuracy在 OOB 上算更可靠但更慢。5.1 predict 报 “factor has new levels” 的踩坑记录现象训练集上拟合好的模型换一份新数据做预测直接报factor g has new levels c整个批量预测任务中断。原因model.frame在训练时把每个因子的水平集合存进了模型的xlevels预测时会用这份记录去重新编码新数据。新数据里出现了训练时没见过的水平无法映射到任何一个哑变量列只能报错。这其实是保护机制不报错就会静默产出一堆错误的预测值。解决训练前就固定水平集合factor(g, levels c(a, b, c))让训练集里暂时没出现的水平也占一个位置或者在预测前做一次droplevels并把未知水平显式归到“其他”类。批量预测脚本里这一条不处理跑一晚上可能在第 3000 条才崩。5.2 预处理写在交叉验证外面的踩坑记录现象交叉验证准确率 0.95上线后掉到 0.7怎么调参都复现不了那个数字。原因标准化、缺失值填补、特征筛选是在全量数据上做的然后才切分训练折和验证折。验证折的均值和缺失值填补信息已经泄漏进了训练过程模型提前“看过”验证集。特征筛选的泄漏更严重用全量数据的相关性挑出来的变量本身就带着标签信息。解决把所有依赖数据的预处理步骤放进重采样的内层。用caret的preProcess参数或者recipes定义一个recipe再交给工作流保证每一折都是先拟合预处理参数、再应用到验证折。判断标准很简单把整个流程套在一个只含单折数据的子集上能不能独立跑通。5.3 na.omit 静默改行数的踩坑记录现象拟合完模型residuals(fit)的长度比原数据少了几十行想把这些残差拼回原表做分析对不上号。原因na.action默认是na.omit含缺失值的行在构建model.frame阶段就被删掉了模型对象里根本不记录它们原来的位置。解决把全局选项设成options(na.action na.exclude)这样residuals、fitted、predict返回的结果会自动用NA占位长度和原数据一致。更稳妥的做法是在建模前显式统计缺失比例缺失率高的列直接排除而不是让它在某个环节悄悄丢行。5.4 set.seed 在并行里失灵的踩坑记录现象单线程跑交叉验证set.seed(42)之后结果完全可复现改成并行后同样的种子每次跑出来的误差都不一样。原因并行 worker 是独立的进程或会话它们不共享主进程的随机数状态各自从系统熵源初始化。主进程那句set.seed对它们不生效。解决用RNGkind(LEcuyer-CMRG)切到支持多流的随机数发生器再用clusterSetRNGStream或者future.seed TRUE给每个 worker 分发独立的子流。另一个常见误用是把set.seed写在循环体内每折都重置成同一个种子结果各折的随机性完全相关交叉验证的方差估计被严重低估。5.5 as.numeric(factor) 拿到编码而不是数值的踩坑记录现象一列看起来是数字的变量as.numeric转完变成 1、2、3 这样的小整数画散点图全挤在一起。原因read.csv或data.frame在某些配置下把字符列自动转成了因子因子在 R 内部是整数编码加一张水平标签表as.numeric作用在因子上取的就是编码。解决as.numeric(as.character(x))先还原成字符再解析更彻底的办法是读入时显式指定stringsAsFactors FALSE和colClasses别让类型推断替你做决定。这个坑在 ID 列、邮编、年份上尤其阴险因为不报错只算错。6. 设计取舍公式接口与矩阵接口R 把哪一头放在了前面R 在这件事上的取舍非常明确公式是一等公民矩阵是二等公民。好处是统计语义可以被完整表达一个y ~ x * g就同时携带了变量、交互、截距和因子编码规则模型函数只需要约定“我吃设计矩阵”。代价是这条链路每一环都有开销model.frame要复制数据model.matrix生成稠密矩阵因子展开后列数可能翻几倍。维度公式接口矩阵接口表达力交互、变换、偏移量直接写需手工构造列内存稠密设计矩阵因子展开膨胀可用稀疏矩阵推断自动带terms、xlevels、对比编码全部自己管典型场景统计建模、小中规模数据高维、需要复用同一矩阵这也解释了为什么惩罚回归这类需要沿正则化路径反复求解的方法普遍提供矩阵接口而不是公式接口同一次分解要在几十个 lambda 上复用不能让model.matrix在每次迭代里重新展开因子。验证这套机制是否真的吃透了有个笨办法很有效随便挑一个模型把它的设计矩阵、梯度和收敛判据分别手算一遍再跟库函数的输出比。三处对得上说明你理解的是机制而不是调用方式。从那以后我每次换新算法都先用十几行样本数据把设计矩阵和梯度打印出来对一遍再上真实数据。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
