简介多元回归模型是数学建模与统计数据分析中的核心方法。这份docx文档围绕某市粮食年销售量与常住人口、人均收入以及肉、蛋、鱼销售量等变量的关系系统记录了完整实验报告适合学习统计建模、经济数据分析或准备数学建模竞赛的学生参考。资源仅含1个docx文件压缩包约205KB轻量易用。文档从因变量的经济意义分析、散点图考察各变量相关性到初始多元线性回归模型建立再借助Matlab逐步回归法剔除X3、X5、X6等不显著变量得到优化模型Y β0 β1X2 β2X4并完成R²、F统计量、P值与置信区间检验最后对特定年份进行预测及经济含义解读步骤清晰附有可复现的Matlab程序附录。当前已有159人学习。通过该文档可掌握多元回归模型的建立、估计、检验和解释变量筛选方法并可直接借鉴其实验报告结构用于课程作业或建模训练。1. 数学建模赛道上的多元回归模型先想清楚再动手数学建模赛题发下来很大一部分最后都会落在“找关系”上哪些因素在影响某个指标影响有多大能不能用一张系数表把结论讲清楚。多元回归模型几乎每年都出现在国赛、华为杯这类比赛里有些题表面写“分类预测”拆解之后仍然要先跑一个回归作为基线。所谓“完整版”的回归不是把数据丢进库函数就结束而是从 docx 里的原始赛题表格出发经过数据清洗、相关性检查、建模、变量筛选、残差诊断最后把结果写回报告文档形成一条可复现的流水线。这套流程既适合拿到赛题后想快速稳住阵脚的参赛者也适用于企业内部用回归做量化分析、还要交付可读文档的工程师。2. 多元回归模型的数据准备把 docx 赛题表格读成可计算的 data2.1 从赛题正文和附件把表格数据取出来赛题通常有两种数据形态第一种直接嵌在 docx 正文里第二种放在附件里作为 Excel 文件。前者要先把表提取成结构化数据后者数据相对干净但仍要先跑缺失值统计。第一步是用 python-docx 把 docx 内嵌表格抓出来我常用的读取代码是这样的import pandas as pd from docx import Document # python-docx 会把文档里的所有表按出现顺序放进 doc.tables doc Document(problem.docx) print(赛题文档内表格数量:, len(doc.tables)) tb doc.tables[0] # 第一张表通常是赛题主数据表 rows [] for row in tb.rows: # cell.text 可能带换行符和前后空格统一 strip 干净 rows.append([cell.text.strip() for cell in row.cells]) df pd.DataFrame(rows[1:], columnsrows[0]) # 首行当列名 print(df.head())这段代码的关键在于doc.tables返回文档中的全部表格对象即使赛题正文里“表1”“表2”之间隔了大段文字也能按序号拿到对应表。rows[1:]是把表头行单独切走如果表格里有合并单元格row.cells长度会不齐常见做法是做前向填充把合并后空白的单元格填成上一格的值再用drop_duplicates清理重复行。第二种情况是数据放在附件里读取通常简单很多但缺失值和脏数据还是要处理# 附件数据统一用 read_excel 读入顺便检查缺测分布 df pd.read_excel(data.xlsx, sheet_nameSheet1, header0) miss df.isna().sum() / len(df) # 按列计算缺失比例 print(miss[miss 0.05]) # 缺失超 5% 的列单独摘出来 # 如果某列的值明显是“文本数字”混排先转 numeric 再处理 df[wind] pd.to_numeric(df[wind], errorscoerce)pd.isna()拿到的是每个位置的缺失标记除以总行数得到列级缺失率。只有缺失率大于 5% 的列才需要专门决策低于这个阈值直接丢行效果通常可以接受。errorscoerce会把无法转成数字的文本变成 NaN这能在进入模型之前显式暴露脏数据而不是让模型报一份看不懂的类型错误。表格数据进回归之前会遇到缺失值怎么补的问题我一般按这张表决策缺失情况处理方案适用前提随机缺失且占比小于 5%直接删除对应行总体样本量不低于 100删除后仍保留多数样本按时间或序号顺序缺失前向填充或线性插值数据本身有顺序插值不破坏时间结构单列缺失超过 20%删除该变量该列不是赛题要讨论的核心解释变量核心解释变量缺失多重插补或用同类变量均值论文中必须写明插补方式和最终样本量2.2 相关性分析热力图先给模型一个预期多元回归有一条隐含前提入选变量至少和因变量有一定关系同时自变量之间不要过分相关。建模前先看相关性矩阵等于给自己一张“预期表”挂在这里后面跑出来的符号和显著性如果和热力图对不上就要回头找原因import seaborn as sns import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] # 让热力图正常显示中文 plt.rcParams[axes.unicode_minus] False # 避免负号变成方块 num_cols [temp, humidity, wind, price, power] corr df[num_cols].corr(methodpearson) # 连续变量用 pearson sns.heatmap(corr, annotTrue, fmt.2f, cmapRdBu_r, linewidths0.5) plt.show()methodpearson适合数值型连续变量如果存在排序型的量比如“等级 1/2/3”改成spearman更合适。看热力图时我只关注两件事解释变量和因变量之间是否存在较高相关性以及解释变量之间是否出现 |r| 0.8 的组合。这里做一层过滤并不是要严格剔除变量而是给后续逐步回归提供一个候选池参考——相关性过高的两个变量同时进模型结论往往不可靠其中一个系数很可能变得不显著甚至变号。3. 多元回归模型的拟合与参数解读statsmodels 一行出完整结果3.1 最小二乘拟合与回归摘要的四个先读项数据整理完之后就可以直接拟合。比赛和业务里我更推荐statsmodels而不是 sklearnstatsmodels 输出的是带检验统计量和 P 值的完整摘要写论文时可以直接引用sklearn 的LinearRegression侧重预测系数倒是能拿但标准误和显著性要自己算。import statsmodels.api as sm # 解释变量和因变量分开便于后续调整变量集合 X df[[temp, humidity, wind]].copy() y df[power] # add_constant 是在 X 前面补一列 1对应回归方程里的截距项 # 不补这一列模型会被强制过原点多数情况下并不合理 X sm.add_constant(X) model sm.OLS(y, X).fit() # 普通最小二乘OLS 即 Ordinary Least Squares print(model.summary())sm.add_constant这一步非常关键漏掉之后截距项被吸收成 0系数会整体偏离真实值。model.summary()的输出里有四个位置值得在论文里展开R-squared表示模型对因变量整体方差的解释比例每个系数的P|t|决定显著性F-statistic检验所有解释变量的系数是否同时为 0Cond. No.是条件数接近 30 以上就要怀疑共线性这是下一章要专门处理的问题。回归摘要和论文写作的对应关系我一般按这张表组织摘要里的指标论文里的措辞判断标准R-squared模型的拟合优度0.6 以上说明解释力可接受但仍要看样本量Pt单个变量对因变量的影响是否显著F-statistic 对应 P 值整体回归方程是否成立小于 0.05说明至少存在一个显著的解释变量coef 的符号影响方向正向还是负向与业务常识不符时优先检查数据是否放错列写论文时很多人直接贴summary()原文截图这并不专业。规范体例是把系数、标准误、P 值和显著性标记做成三线表summary()只作为中间检查工具。把关键结果抽成 DataFrame 是顺手的事# 抽取系数、标准误和 P 值格式化出三线表雏形 coef_df pd.DataFrame({ 变量: model.params.index, 系数: model.params.values, 标准误: model.bse.values, P值: model.pvalues.values }) print(coef_df)model.params返回每个变量的点估计model.bse是参数的标准误model.pvalues是显著性检验的 P 值。这三组数据在 OLS 结果里是成对出现的导出时保持顺序一致就不会写串。3.2 类别变量和交互项把赛题里的定性条件写进公式赛题经常给出“工作日/节假日”或“旺季/淡季”这样的定性变量直接把它们当数字塞进模型没有意义。statsmodels 的公式接口支持哑变量转换C()包装一个变量会自动生成 k-1 个虚拟变量import statsmodels.formula.api as smf # C() 表示把 season 转为哑变量season 有 4 个水平就会生成 3 列 # C(season):humidity 是交互项表示湿度对因变量的影响随季节变化 model2 smf.ols(power ~ temp C(season) C(season):humidity, datadf).fit() print(model2.summary())smf.ols使用 patsy 公式语法写法更贴近回归方程本身也省去手动创建哑变量的过程。交互项的解读要小心交互项显著只能说明“斜率在不同组之间存在差异”不能直接解读成“某变量在该组内单独起作用”。分组的基准组是C()生成的第一个水平论文里解释系数时要先说明以哪一组作为基准否则读者会误解系数的参照点。4. 多元回归模型的变量筛选逐步回归帮你留下该留的变量4.1 全模型的问题在哪把能拿到的变量一次性全放入模型结果不一定更好。变量过多意味着参数空间变大在赛题这种样本量有限的情况下模型很容易把噪声吸收进系数里同时高度相关的变量会让系数估计值的方差变大某个变量单看可能不显著删掉它之后另一个变量又变得显著。这就是我们通常不用全模型的原因而是用变量筛选先框定一个“解释得动、也写得清”的变量集合。筛选的评判指标常用 AIC赤池信息准则。AIC 在残差平方和的基础上加入一项参数个数的惩罚变量加得越多、惩罚越大所以筛选会倾向停在某个平衡点上。AIC 越小模型越优BIC 则是惩罚更严格的版本样本量不大时直接用 BIC 会把变量删得过多建模竞赛里用 AIC 更多。4.2 前向逐步回归的循环实现前向逐步回归的流程是从空模型开始逐个尝试把候选变量加入模型计算每次的 AIC挑出能让 AIC 下降最多的变量真正加入重复这个过程直到加入任何剩余变量都不能让 AIC 继续下降。用循环实现并不复杂import statsmodels.api as sm def forward_aic(data, target, candidate): 前向逐步回归按 AIC 最小原则逐个选入变量 data: 数据集target: 因变量列名candidate: 候选解释变量列表 selected [] old_aic None # 上一次模型的 AICNone 代表空模型 while candidate: best_col, best_aic None, None for col in candidate: trial selected [col] # 假设把 col 加进当前模型 trial_model sm.OLS(data[target], sm.add_constant(data[trial])).fit() if best_aic is None or trial_model.aic best_aic: best_aic, best_col trial_model.aic, col # 如果加入最佳候选之后 AIC 反而上升就停止筛选 if old_aic is None or best_aic old_aic: selected.append(best_col) candidate.remove(best_col) old_aic best_aic else: print(AIC 不再下降停止) break return selected final_vars forward_aic(df, power, [temp, humidity, wind, price]) print(最终入选变量:, final_vars)这段代码里有两个容易忽略的细节。第一每次尝试候选变量时都会重新拟合一次模型总拟合次数等于“候选变量数 × 迭代轮数”所以只适合候选变量在 20 个以内的情况候选变量上百个时改用 Lasso 这类带惩罚的回归更现实。第二old_aic记录上一轮已选模型的 AIC比较时必须在同一样本上进行中途不能丢行或换样本否则前后 AIC 不具备可比性。自动筛选只是机器视角最终模型还要过一遍三问检查检查问题具体做法符号是否符合赛题背景对比回归系数的正负号与业务常识、文献结论核心变量有没有被筛掉赛题明确要求讨论的变量必须保留不管 AIC 怎么说加入交互项后符号是否突变在主效应模型和交互模型之间各跑一遍确认不是共线性所致5. 多元回归模型诊断残差、VIF 与稳健标准误5.1 残差四件套先过一遍回归模型的推断建立在残差满足若干假设的基础上。论文里不用把每个检验都写全但评委复查时残差图是最常用的快速判断手段。拿到初始模型后我习惯一次画 Q-Q 图和残差对拟合值散点图再配合 Shapiro 检验做数值判断import scipy.stats as stats import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(11, 4)) # Q-Q 图散点越贴在对角线上残差越接近正态分布 stats.probplot(model.resid, distnorm, plotaxes[0]) axes[0].set_title(残差 Q-Q 图) # 残差对拟合值散点观察是否存在单调的喇叭口或弯曲 axes[1].scatter(model.fittedvalues, model.resid, alpha0.6) axes[1].axhline(0, colorgray, linestyle--) axes[1].set_xlabel(拟合值) axes[1].set_ylabel(残差) plt.tight_layout() plt.show() w, p_value stats.shapiro(model.resid) print(残差 Shapiro P 值 , p_value) # P 0.05 认为残差近似正态Q-Q 图不需要复杂的统计知识就能看散点整体贴着 45 度线说明残差没有明显偏离正态散点两端出现“S”形说明尾部比正态更厚。残差对拟合值散点如果呈现以 0 为中心的均匀条带说明线性关系和同方差假设基本成立出现从左到右扩散的喇叭口就说明存在异方差对应 5.3 节要处理的场景。stats.shapiro对样本量有限制n 超过 5000 时检验功效变得很强任何小偏差都会被判定显著这时改看残差直方图更可靠。建模比赛里评委通常不会只看 P 值而是先看散点形态再决定是否接受正态性结论所以图要有、文字结论也要有。5.2 VIF 共线性检查判断系数能不能站得住相关性热力图只能暴露两两相关多个变量之间的线性关系要靠方差膨胀因子 VIF 来判断。VIF 的做法是把每个自变量当作因变量、用其余自变量对它做一次回归看这个回归的 R² 有多高R² 越高说明该变量越能被其他变量解释它的系数方差就会被放大得越厉害from statsmodels.stats.outliers_influence import variance_inflation_factor X_use df[final_vars] # final_vars 是逐步回归选出来的变量 X_const sm.add_constant(X_use) # VIF 计算同样需要常数项 vif_df pd.DataFrame({ 变量: X_const.columns, VIF: [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])] }) print(vif_df)variance_inflation_factor(X, i)的第二个参数是列下标传 0 时对应常数项它的 VIF 没有解释意义只看后面自变量即可。VIF 的阈值习惯用 10 作线5 到 10 之间则要结合变量重要性权衡VIF 区间处理建议小于 5正常维持原模型5 到 10系数方向没有明显跳动则保留否则考虑删一个大于等于 10优先删除 VIF 最高的变量再重新拟合观察效果5.3 异方差的补救换稳健标准误散点图出现喇叭口时系数的点估计不会变但标准误被低估结果是 P 值偏小、显著性被夸大。这种情况不需要马上放弃 OLS更务实的做法是换异方差稳健标准误# HC1 是赛题里最常用的稳健标准误对中小样本比较温和 model_robust model.get_robustcov_results(cov_typeHC1) print(model_robust.summary())cov_typeHC1对应 Stata 回归里的默认稳健标准误在比赛和论文里是主流选择如果数据集存在明显的离群点换成HC3对小样本更敏感。替换成稳健标准误之后系数保持相同P 值和置信区间会变宽这才是更真实的显著性水平。写入报告时应注明“标准误为异方差稳健标准误”评委会注意到这一步。6. 多元回归模型结果进 docx系数表导出与数值一致性验证6.1 一键生成规范的三线表模型定稿后把系数表手动敲进 Word 不但慢而且容易抄错。直接用 python-docx 在报告末尾追加表格可以让格式统一、数字和模型输出保持一致from docx import Document def write_coef_table(doc_path, coef_df): 把系数表追加到 docx 文档末尾写成三线表雏形 coef_df 至少包含 变量、系数、标准误、P值 四列 report Document(doc_path) table report.add_table(rowscoef_df.shape[0] 1, cols4) table.style Table Grid headers [变量, 系数, 标准误, 显著性] for j, h in enumerate(headers): table.rows[0].cells[j].text h for i, row in coef_df.iterrows(): # 按 P 值生成论文通用的星号标记 if row[P值] 0.01: sig ** elif row[P值] 0.05: sig * else: sig values [str(row[变量]), f{row[系数]:.4f}, f{row[标准误]:.4f}, sig] for j, v in enumerate(values): table.rows[i 1].cells[j].text v report.save(doc_path) return table显著性星号在论文里要统一写在表注中** p0.01* p0.05。系数保留四位小数即可保留太多会让表显得拥挤也没必要。表格追加完成后在正文里补一句“估计结果见表 1”这类引导并给出回归方程形式例如power β0 β1·temp β2·humidity β3·wind ε解释时把每个 β 对应到表里的行。6.2 提交前把 docx 里的数字和模型回读比对最后一步是验证。模型可能在交稿前又调过报告里的数字却还停留在上一版这问题在比赛里很常见。做法是把 docx 里的系数抽出来和当前模型输出比对偏差超过设定阈值就报警from docx import Document def verify_report(docx_path, expected_dict): 从 docx 表格中抽取系数与最新模型输出比对 expected_dict: {变量名: 最新系数值} doc Document(docx_path) coef_dict {} for table in doc.tables: for row in table.rows[1:]: cells [c.text.strip() for c in row.cells] if len(cells) 2 and cells[0] in expected_dict: try: coef_dict[cells[0]] float(cells[1]) except ValueError: continue # 跳过非数字单元格 return {k: coef_dict.get(k) - expected_dict[k] for k in expected_dict} # 把最新 coef_df 转成字典回读 report.docx 做差值检查 expected dict(zip(coef_df[变量], coef_df[系数])) diff verify_report(report.docx, expected) print(docx 与当前模型的系数差值:, diff)这个回读函数能识别变量名并跳过表头和非数字行所以把报告里任何一张系数表都过一遍偏差一目了然。比赛前一天把新模型重跑一遍让报告完整再生一次再执行这个校验基本能堵住“模型更新了、文档忘了同步”的低级错误这也是把回归做完整之后最值得花时间的环节。本文还有配套的精品资源点击获取