简介一份面向统计学、数据分析及Matlab建模初学者的多元线性回归完整实现文档。文档以研究生教材《数理统计》例4.4.1为背景在回归方程F检验基础上扩展了对各回归系数x1、x2、x3的t检验支持用户自定义显著性水平α提高计算精度。程序中包含数据读取、最小二乘参数估计、设计矩阵构建、F检验与t检验等关键步骤输出结果整洁易读。使用者只需更换Excel数据文件即可应用于其他维度数据集移植性强。资源共1个docx文件压缩包大小50KB文件包含程序说明、数据输入格式、完整Matlab源码及运行结果示例适合需要对照教材理解回归显著性检验原理、并希望获得可直接运行模板的读者。目前已有246人学习下载。1. 多元线性回归显著性检验到底查什么先看p值还是先看模型拿到“多元线性回归及显著性检验Matlab程序.docx”这个标题我第一反应是这不是一个求系数的任务而是一个做“归因判断”的任务。多元线性回归的Matlab程序很容易写调用一次 regress 就能拿到 b 和 stats但标题里特意带上“显著性检验”意味着要回答的不只是“y 和 x 是什么关系”而是“哪些 x 真正影响 y这种影响在统计上是否站得住脚”。在实际项目中这个区别很现实论文外审会问回归表和显著性标记业务评审会问某个系数为什么是 2.3以及它为什么能推广到新数据。这类程序适合两类读者一类是写课程作业、毕业论文或实证论文的人需要把回归结果按规范格式呈现另一类是做数据分析或算法验证的工程师手上有一组带标签的观测数据想快速判断特征是否有效。下面我按自己常用的做法把模型建立、整体显著性检验、系数显著性检验、诊断和封装串成一条完整路径。你会看到同一组数据用 regress 和 fitlm 两种方式跑出来结果一致但呈现思路不同而显著性检验才是决定模型能否交付的关键。2. 用Matlab跑通多元线性回归程序regress与fitlm的数据组织和参数表多元线性回归在 Matlab 里的入口很多最常见的是统计工具箱里的 regress 和 fitlm。前者面向矩阵运算输出维度直接适合写进自动化程序后者面向表格数据自带回归报告适合交互式分析和出图。在开始之前要做一件事确认当前环境已经安装了 Statistics and Machine Learning Toolbox。无论你是刚完成 Matlab 下载还是使用公司内网安装版本工具箱缺失时调用 regress 会报出 Undefined function这个问题在安装步骤里常被忽略值得先检查。2.1 多元回归的数据组织设计矩阵必须有一列截距项regress 的调用格式是[b,bint,r,rint,stats] regress(y,X,alpha)其中 y 是 n×1 列向量X 是 n×(p1) 设计矩阵。这里的“p1”是很多新手第一个踩坑点regress 不会自动帮你补截距列它把 X 的第一列当作 beta0 的系数列。也就是说如果你想拟合 y beta0 beta1x1 beta2x2X 必须是[ones(n,1), x1, x2]而不是直接把原始变量拼在一起。fitlm 的设计刚好相反它默认包含截距项从 table 类型的表格数据出发用公式y ~ x1 x2 x3描述模型。因此在数据组织上使用 fitlm 时只需要把变量放在一个 table 里函数会自动为截距项分配一行参数。两种方式没有好坏之分我通常的做法是需要在循环里批量跑模型、后续还要做残差诊断时用 regress需要快速查看模型报告、画预测图和残差图时用 fitlm。2.1.1 为什么 regress 要求自己加常数项列这和 regress 的底层实现有关。它直接对线性方程组 Xb y 做最小二乘求解X 的每一列都对应 b 的一个分量。如果你不给常数项列求出的回归直线会被强制穿过原点模型实际上是 y beta1x1 beta2*x2截距被丢弃。这在变量中心化且理论上确实过原点的问题里没问题但绝大多数业务数据并不满足这个前提。一个简单的检查方法是在调用 regress 前查看size(X,2) - size(unique(X,rows),2)如果设计矩阵没有全 1 列结果会很反常。2.1.2 一个可复现的示例数据集为了演示我构造一组已知真实关系的样本数据y 与三个变量 x1、x2、x3 满足线性关系并加入随机噪声。这样后续的显著性检验结果能和我们预设的系数对照方便判断程序是否正确。样本量设为 30足以支撑 3 个自变量的检验。2.2 regress 的返回值与 stats 参数含义regress 的五个输出参数中b 是回归系数向量bint 是每个系数在置信水平 1−alpha 下的置信区间r 是残差rint 是每个残差的置信区间stats 是模型统计量汇总。前四个参数对后续残差图很有用第五个参数则是显著性检验的核心数据。下面是一张参数速查表适合放在程序注释旁边。输出参数维度说明b(p1)×1回归系数b(1) 是截距项bint(p1)×2每个系数的置信区间第一列为下限第二列为上限rn×1残差向量等于 y 减去拟合值rintn×2残差置信区间用于 rcoplot 画诊断图stats1×4依次为 R²、F 统计量、F 对应的 p 值、误差方差估计要注意 stats 的四个分量顺序固定stats(1) 是确定系数 R²stats(2) 是回归方程显著性 F 统计量stats(3) 是 F 检验的 p 值stats(4) 是均方误差的估计值。alpha 参数只影响 bint 和 rint 的置信水平不会改变 b 和 stats 中的 F、p 值。所以在写程序时把 alpha 单独定义成一个变量方便后续换成 0.01 或 0.1 观察显著性变化。2.3 最小可运行程序示例% 多元线性回归及显著性检验最小程序 clear; clc; rng(0); % 固定随机种子保证可复现 % 生成 30 个样本3 个自变量 x1 rand(30,1) * 10; x2 rand(30,1) * 5; x3 rand(30,1) * 2; % 构造设计矩阵第一列为常数项 X [ones(30,1), x1, x2, x3]; % 真实关系为 y 3 2*x1 1.5*x2 - 1.2*x3 y 3 2*x1 1.5*x2 - 1.2*x3 randn(30,1) * 0.8; alpha 0.05; [b, bint, r, rint, stats] regress(y, X, alpha); disp(回归系数 b); disp(b); disp(每个系数的 95% 置信区间 bint); disp(bint); disp(模型统计量 [R^2, F, p, 误差方差估计]); disp(stats);这段代码先用固定随机种子生成数据再调用 regress。核心是X [ones(30,1), x1, x2, x3]它把截距项作为第 1 列。运行后可以看到 b 大约在 3、2、1.5、-1.2 附近bint 包围真实值stats 第二个数是很大的 F 值第三个数很小说明整体回归显著。这里每个变量的系数都以真实值为中心是因为样本噪声不大且样本量足够覆盖 3 个自由度。2.4 用 fitlm 输出回归表和 regress 相互验证如果觉得 regress 的结果不够直观可以用 fitlm 生成带文字描述的回归表。先把变量放入 table再指定公式。fitlm 会输出每个系数的估计值、标准误、t 统计量和 p 值底部还会给出 F 检验结果和 regress 的 stats 对应。% 用 fitlm 生成结构化回归报告 tbl table(x1, x2, x3, y, VariableNames, {x1, x2, x3, y}); mdl fitlm(tbl, y ~ x1 x2 x3); disp(mdl);从 fitlm 输出中可以同时看到回归系数和显著性标记。rows 对应 Intercept、x1、x2、x3Columns 包含 Estimate、SE、tStat、pValue。这个输出适合直接截图放进报告而 regress 的输出更适合后续程序化处理。两种方式用同一份数据时结果一致这可以作为一个验证手段先用 regress 计算再用 fitlm 复现如果系数差异超过浮点精度范围说明设计矩阵或数据组织存在问题。3. 回归方程显著性检验F统计量、决定系数R²和程序判定模型有没有解释力不是看 b 是否为 0而是看所有自变量的联合解释力是否显著。这一步对应的是回归方程显著性检验也就是通常说的 F 检验。它的作用域是“整个方程”而不是某一个具体变量。很多初学者把注意力放在 R² 上觉得 R² 达到 0.9 就把结果报上去这是不够的。R² 描述的是拟合优度F 检验描述的是统计显著性二者一个看数据解释比例一个看抽样波动下的置信度必须同时出现在程序输出里。3.1 F检验的原理从方差分解到统计量计算多元线性回归的 F 检验把总离差平方和 SST 分解成回归平方和 SSR 和残差平方和 SSE。原假设 H0 是所有自变量系数同时为 0备择假设是至少有一个系数不为 0。F 统计量定义为 SSR 除以回归自由度 p再除以 SSE 除以残差自由度 n−p−1。当 F 值较大超过 F 分布临界值时拒绝原假设认为回归方程整体显著。regress 已经把这个过程封装进 stats第二个元素就是 F 值第三个元素是对应的 p 值。用代码表示F 值也可以从基础量计算出来这样能加深理解。给定回归结果后n length(y)p size(X,2) - 1残差平方和由SSE sum(r.^2)得到那么% 从 regress 的结果中手动计算 F 统计量 n length(y); p size(X, 2) - 1; ybar mean(y); SST sum((y - ybar).^2); SSE sum(r.^2); SSR SST - SSE; F (SSR / p) / (SSE / (n - p - 1)); pF 1 - fcdf(F, p, n - p - 1);把这段结果和 stats(2)、stats(3) 对比会发现两者完全一致。fcdf 是 F 分布的累积分布函数1 - fcdf(F, p, n-p-1)得到右侧尾部概率也就是 p 值。F 值越大p 值越小这解释了为什么很多回归报告里只看 p 值就能判断整体显著性。注意自由度 n−p−1 是残差自由度p 是自变量个数不含截距项。3.2 R²与调整R²在显著性判断中的角色R² 等于 1 减去 SSE 除以 SST表示自变量解释了 y 中变异的比例。但它有一个固有缺陷只要增加自变量R² 就会上升哪怕新变量和 y 毫无关系。调整 R² 把自由度惩罚引入公式在 p 增大时抵消一部分 R² 的虚增。因此在做显著性检验时我通常同时输出 R² 和调整 R²并观察二者差距是否过大。如果差距明显说明模型里很可能混入了不显著的解释变量。指标计算公式使用建议R²1 − SSE/SST描述拟合优度不惩罚变量数量调整 R²1 − (1−R²)(n−1)/(n−p−1)比较不同自变量个数的模型时更可靠F 统计量(SSR/p)/(SSE/(n−p−1))判断整体显著性伴随 p 值误差方差估计SSE/(n−p−1)即 stats(4)用于计算系数标准误调整 R² 的计算在 Matlab 里一行就能完成adjR2 1 - (1-stats(1)) * (n-1) / (n-p-1)。把它放在报告输出中优先级高于单纯的 R²。当 R² 为 0.95、调整 R² 只有 0.60 时说明模型塞进了太多低解释力变量此时再看 F 检验和单个系数的 t 检验往往会发现大量不显著项。3.3 用程序判定整体显著性并输出结论实际使用中我习惯在回归程序后加一段自动判定代码把 stats 的三个关键量格式化输出。比如设置 alpha 0.05当stats(3) alpha时可以判断模型整体显著否则需要检查数据是否存在严重共线性、变量个数是否过多或者原始变量是否需要变换。下面这段代码可以直接嵌在上一章的程序后% 整体显著性判定 if stats(3) alpha fprintf(模型整体显著F%.3fp%.4g\n, stats(2), stats(3)); else fprintf(模型整体不显著F%.3fp%.4g\n, stats(2), stats(3)); end fprintf(R^2%.4f调整R^2%.4f\n, stats(1), 1-(1-stats(1))*(n-1)/(n-p-1));这里有个细节p 值在极小的时候disp会显示成 1.2345e-12用%.4g格式化后更友好。整体显著不等于每个自变量都显著这一结论直接引出下一章的单变量显著性检验。如果整体不显著后面的 t 检验单独看意义也不大应当先回到变量选择和数据质量问题上。因此程序里要把整体检验放在变量筛选之前。4. 回归系数显著性检验t检验、置信区间和变量筛选F 检验回答的是“模型整体有没有用”t 检验回答的是“每个变量是不是真的有用”。二者互补但不互相替代。一个典型的误读场景是F 检验显著业务方就认为所有自变量都该留下结果发现 x3 的 p 值是 0.35根本没有通过显著性检验。反过来也可能出现整体不显著但某个变量单独显著的情况这通常由多重共线性或样本量不足造成。多元线性回归的显著性检验最终要落在每个解释变量头上。4.1 为什么必须对每个回归系数做 t 检验t 检验的原假设是 H0: beta_j 0备择假设是 beta_j 不等于 0。检验统计量是参数估计值除以它的标准误t b_j / se(b_j)其中 se(b_j) 是 b_j 的抽样标准差由噪声方差和设计矩阵共同决定。回归方程里有 p 个自变量就有 p 个这样的 t 检验。regress 返回的 bint 就是基于 t 分布构造的置信区间如果区间跨过 0说明在给定置信水平下无法排除该系数为 0 的可能。这一环节必须结合自由度来解释。残差自由度是 n−p−1t 统计量服从自由度 n−p−1 的 t 分布。样本量越小临界值越大相同 t 值下显著性越难显现。比如 n10、p3 时自由度只有 695% 置信水平下的临界 t 值超过 2.4而 n100 时临界值接近 1.98。所以在写程序时不要用固定的 1.96 作为 t 临界值应当用 tinv 或 tcdf 按实际自由度计算。4.2 从 regress 结果手动计算 t 统计量和 p 值regress 没有直接返回 t 统计量但 stats(4) 给出了误差方差估计值 MSE。系数协方差矩阵的估计是 MSE 乘以 XX 的逆矩阵对角线开方就是每个系数对应的标准误。拿到标准误后t 统计量和 p 值可以用以下代码得到% 手算 t 统计量与 p 值和 bint 互相验证 [b, bint, ~, ~, stats] regress(y, X, alpha); n size(X, 1); k size(X, 2) - 1; df n - k - 1; MSE stats(4); C inv(X * X); se sqrt(diag(MSE * C)); % 每个系数的标准误 tstat b ./ se; pval 2 * (1 - tcdf(abs(tstat), df)); fprintf(变量\t系数\t标准误\tt值\tp值\n); for j 1:length(b) fprintf(x%d\t%.4f\t%.4f\t%.4f\t%.4g\n, j-1, b(j), se(j), tstat(j), pval(j)); end代码中用tcdf计算 t 分布累积概率双尾 p 值乘 2。注意这里 b(1) 对应截距项它的 p 值一般不作为变量筛选依据只有业务上关心截距是否为 0 时才需要解读。严格来说inv(X*X)在 X 存在严重复共线性时不稳定更稳的写法是使用pinv(X*X)或直接调用regress内部的结果对于常规数据这个公式便于理解。4.3 用置信区间和 t 检验结果快速筛选变量除了看 p 值bint 是一个更直观的判定工具。某个系数的 bint 区间如果落在 0 的同侧说明下限和上限同号此时区间不含 0可以拒绝原假设。如果下限为负、上限为正则意味着系数在正负之间波动不能确认它对 y 有真实作用。可以用符号判断写一段筛选代码% 用 bint 识别不显著的回归系数 coef_sign sign(bint(:, 1)) .* sign(bint(:, 2)); insig_idx find(coef_sign 0); % 区间包含 0不显著 if ~isempty(insig_idx) fprintf(以下系数未通过显著性检验); fprintf(%d , insig_idx - 1); fprintf(\n); else fprintf(所有系数均通过显著性检验\n); end对于上述模拟数据x1、x2、x3 通常都会通过检验。为了让筛选逻辑更通用下面用一张示意表说明不同情况下的判定结果。实际项目里x1 的系数可能是显著的x2 处于临界x3 不显著这时应该把 x3 从模型中去掉并重新拟合。变量bsetp判定x12.030.316.550.0002显著x21.470.364.080.0010显著x3-0.210.30-0.700.4950不显著x4-1.341.18-1.140.2680不显著需要特别提醒的是一次检验不显著不代表变量本身和 y 无关也可能是样本量不够或者该变量与其他变量高度相关。在删除变量之前先回到数据相关性矩阵和 VIF 检查这一步可以避免误删有效特征。如果只是做预测保留不显著变量虽然会增加参数复杂度但不一定降低预测精度只有要解释变量作用时才必须严格执行剔除规则。5. 显著性检验后的诊断残差图、多重共线性与stepwisefit显著性检验做完并不代表模型可以交付。回归模型的可靠性建立在几个假设上残差独立、方差齐性、残差近似正态以及自变量之间不存在严重多重共线性。如果这些假设被踩破刚才的 F 和 t 检验都会失真。最常见的例子是两个自变量相关系数达到 0.95模型整体 F 值很大但每个系数的 t 检验都不显著原因是系数估计的标准误被共线性放大。此时如果直接看 p 值删除变量会得出完全错误的结论。5.1 残差与异常点rcoplot和normplot怎么看regress 的 r 和 rint 是专门喂给残差图的。rcoplot(r, rint) 会画出每个样本的残差竖线和置信区间如果某条残差置信区间不跨过 0说明该样本对拟合结果影响较大可以视为异常点。normplot(r) 用于检查残差是否服从正态分布如果散点大致落在一条直线上说明正态假设合理如果两端发散严重后续 p 值的可信度就要打折。% 显著性检验后的残差诊断 figure; rcoplot(r, rint); % 残差及置信区间图 title(残差置信区间诊断); figure; normplot(r); % 正态概率图 title(残差正态性检验);在 rcoplot 中看到个别样本的区间明显偏离其他样本时不要急着删除。先检查该样本的原始数据是否存在录入错误再考虑是否属于业务上的特殊情况。一段稳健的做法是删除异常点后重新跑回归对比前后系数和显著性是否发生剧烈变化。如果变化不大说明模型稳定如果显著性结论反转说明原结果本身脆弱需要在报告里如实说明。5.2 多重共线性检验计算VIF替代直接删变量多重共线性的常用指标是方差膨胀因子 VIF。对第 j 个变量先把其他自变量作为输入对这种变量做回归得到决定系数 R_j²然后 VIF_j 1 / (1 − R_j²)。VIF 大于 10 通常被认为存在严重共线性5 到 10 之间需要警惕。计算 VIF 的代码可以这样写% 计算各自变量的 VIF Xall [x1, x2, x3]; [nobs, px] size(Xall); VIF zeros(1, px); for j 1:px Xm Xall; Xm(:, j) []; % 去掉当前变量 Xm [ones(nobs, 1), Xm]; % 补截距列 [~, ~, ~, ~, stat_j] regress(Xall(:, j), Xm); Rj2 stat_j(1); VIF(j) 1 / (1 - Rj2); end disp(各变量 VIF); disp(VIF);Xm 不需要截距之外的额外处理关键是用 regress 对当前变量做回归时它的设计矩阵要包含全 1 列。若某个变量的 VIF 超过 10我通常先在相关矩阵里找到与之相关系数最高的变量然后根据业务含义决定剔除哪一个。也可以保留该变量但改用岭回归不过岭回归的系数不再是无偏估计显著性检验的解释会变复杂。5.3 用stepwisefit自动筛选显著变量如果自变量数量较多手动检验每个 t 值很繁琐。Matlab 的 stepwisefit 可以基于显著性检验自动选择变量参数 penter 控制变量进入模型的 p 值阈值premove 控制变量被剔除的阈值。一般 penter 设 0.05premove 设 0.10使变量一旦进入就不要因为轻微波动立刻离开。% 逐步回归基于显著性检验自动筛变量 [inmodel, ~, ~] stepwisefit(Xall, y, penter, 0.05, premove, 0.10); disp(保留的变量索引); disp(find(inmodel));inmodel 是逻辑向量1 表示该变量保留在模型中。stepwisefit 返回的模型用全部样本重新估计得到的系数和 p 值可以直接作为最终报告。需要注意的是逐步回归的自动搜索会放大偶然显著性尤其在变量多、样本少的情况下最终入选的变量可能在一组新样本中不再显著。因此它只能作为变量初筛工具最终模型还是要回到业务判断和残差诊断上。6. 把多元回归程序封装成函数显著性标记与结果导出回顾前面的步骤从 regress 到显著性判断再到诊断零散的脚本在多个项目里复用时会很别扭。常见做法是封装一个函数输入 y、X 和显著性水平输出结构化的回归结果表把 t 值、p 值、置信区间和显著性标记一次性算好。这样无论是课程报告还是业务分析一个命令就能生成结论。function tbl mlr_report(y, X, alpha) % 多元线性回归及显著性检验结果表 if nargin 3 alpha 0.05; end [b, bint, ~, ~, stats] regress(y, X, alpha); n size(X, 1); k size(X, 2) - 1; df n - k - 1; MSE stats(4); se sqrt(diag(MSE * inv(X * X))); tval b ./ se; pval 2 * (1 - tcdf(abs(tval), df)); % 显著性星号标记 sig cell(size(b)); for j 1:numel(b) if pval(j) 0.01 sig{j} **; elseif pval(j) 0.05 sig{j} *; elseif pval(j) 0.1 sig{j} .; else sig{j} ; end end % 构造变量名 if all(X(:, 1) 1) varnames {Intercept}; start 2; else varnames {}; start 1; end for j start:size(X, 2) varnames{end 1} sprintf(x%d, j - start 1); end tbl table(b, bint(:, 1), bint(:, 2), se, tval, pval, sig, ... VariableNames, {beta, CI_low, CI_up, SE, t, p, sig}, ... RowNames, varnames); end调用时直接传入数据和设计矩阵即可tbl mlr_report(y, [ones(30,1), x1, x2, x3], 0.05); disp(tbl);。输出表里每一行是一个回归系数p 值小于 0.01 显示双星小于 0.05 显示单星小于 0.1 显示点没有星号说明不显著。接下来可以把这张表用 writetable 导出到 Excel 或 CSV方便写论文或做汇报writetable(tbl, regression_result.xlsx, WriteRowNames, true)。需要注意X*X在变量量纲差距极大时可能产生较大数值误差如果发现标准误异常可以把 inv 换成 pinv 并观察结果是否稳定。对于样本量小于 30 的小样本数据建议把 alpha 提高到 0.1或者用 fitlm 里的robust选项配合稳健标准误重复一遍如果两种方法给出的显著性标记一致这个模型才值得写进结论里。本文还有配套的精品资源点击获取