R机器学习公式接口与设计矩阵：一行代码背后的原理与避坑-开发者社区

简介：面向希望用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 2

factors这个矩阵是整套机制的关键：行是原始变量，列是模型项，数字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 NA

NA不是拟合失败，而是 QR 分解在选主元时判定这一列与前面对应列线性相关，把它踢出了独立列集合。判定阈值tol = 1e-7是相对于原始列范数的，不是绝对值，所以量纲变化不会轻易改变判定结果。

2.3 QR 分解：为什么 lm 不去求 (X'X) 的逆

教科书上的最小二乘解是(X'X)^{-1} X'y，但 R 的lm.fit走的是另一条路：底层 C 代码Cdqrls对设计矩阵做 Householder QR 分解，解R b = Q'y。原因是条件数会被平方。如果X的条件数是 10^8，X'X的条件数就是 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计算标准误时用的也是这个R：chol2inv(R)得到(X'X)^{-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 b

IRLS 的一轮可以写成：先算线性预测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 里，流程是：对每个预测变量，把样本按该变量排序，只在相邻的不同取值中点处试探分裂，然后计算左右子节点的加权不纯度。分类树默认用 Gini，1 - 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) # MeanDecreaseGini

mtry的默认值随任务类型变，分类是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("L'Ecuyer-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`、对比编码	全部自己管
典型场景	统计建模、小中规模数据	高维、需要复用同一矩阵