简介：多元回归模型是数学建模与统计数据分析中的核心方法。这份docx文档围绕某市粮食年销售量与常住人口、人均收入以及肉、蛋、鱼销售量等变量的关系，系统记录了完整实验报告，适合学习统计建模、经济数据分析或准备数学建模竞赛的学生参考。资源仅含1个docx文件，压缩包约205KB，轻量易用。文档从因变量的经济意义分析、散点图考察各变量相关性，到初始多元线性回归模型建立，再借助Matlab逐步回归法剔除X3、X5、X6等不显著变量，得到优化模型Y = β0 + β1X2 + β2X4，并完成R²、F统计量、P值与置信区间检验，最后对特定年份进行预测及经济含义解读，步骤清晰，附有可复现的Matlab程序附录。当前已有159人学习。通过该文档可掌握多元回归模型的建立、估计、检验和解释变量筛选方法，并可直接借鉴其实验报告结构用于课程作业或建模训练。
1. 数学建模赛道上的多元回归模型：先想清楚再动手
数学建模赛题发下来，很大一部分最后都会落在“找关系”上：哪些因素在影响某个指标，影响有多大，能不能用一张系数表把结论讲清楚。多元回归模型几乎每年都出现在国赛、华为杯这类比赛里，有些题表面写“分类预测”，拆解之后仍然要先跑一个回归作为基线。所谓“完整版”的回归，不是把数据丢进库函数就结束，而是从 docx 里的原始赛题表格出发，经过数据清洗、相关性检查、建模、变量筛选、残差诊断，最后把结果写回报告文档，形成一条可复现的流水线。这套流程既适合拿到赛题后想快速稳住阵脚的参赛者，也适用于企业内部用回归做量化分析、还要交付可读文档的工程师。
2. 多元回归模型的数据准备：把 docx 赛题表格读成可计算的 data
2.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:], columns=rows[0]) # 首行当列名 print(df.head())这段代码的关键在于doc.tables返回文档中的全部表格对象，即使赛题正文里“表1”“表2”之间隔了大段文字，也能按序号拿到对应表。rows[1:]是把表头行单独切走；如果表格里有合并单元格，row.cells长度会不齐，常见做法是做前向填充：把合并后空白的单元格填成上一格的值，再用drop_duplicates清理重复行。
第二种情况是数据放在附件里，读取通常简单很多，但缺失值和脏数据还是要处理：
# 附件数据统一用 read_excel 读入，顺便检查缺测分布 df = pd.read_excel("data.xlsx", sheet_name="Sheet1", header=0) miss = df.isna().sum() / len(df) # 按列计算缺失比例 print(miss[miss > 0.05]) # 缺失超 5% 的列单独摘出来 # 如果某列的值明显是“文本+数字”混排，先转 numeric 再处理 df["wind"] = pd.to_numeric(df["wind"], errors="coerce")pd.isna()拿到的是每个位置的缺失标记，除以总行数得到列级缺失率。只有缺失率大于 5% 的列才需要专门决策，低于这个阈值直接丢行，效果通常可以接受。errors="coerce"会把无法转成数字的文本变成 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(method="pearson") # 连续变量用 pearson sns.heatmap(corr, annot=True, fmt=".2f", cmap="RdBu_r", linewidths=0.5) plt.show()method="pearson"适合数值型连续变量；如果存在排序型的量，比如“等级 1/2/3”，改成spearman更合适。看热力图时我只关注两件事：解释变量和因变量之间是否存在较高相关性，以及解释变量之间是否出现 |r| > 0.8 的组合。这里做一层过滤并不是要严格剔除变量，而是给后续逐步回归提供一个候选池参考——相关性过高的两个变量同时进模型，结论往往不可靠，其中一个系数很可能变得不显著甚至变号。
3. 多元回归模型的拟合与参数解读：statsmodels 一行出完整结果
3.1 最小二乘拟合与回归摘要的四个先读项
数据整理完之后就可以直接拟合。比赛和业务里我更推荐statsmodels而不是 sklearn：statsmodels 输出的是带检验统计量和 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检验所有解释变量的系数是否同时为 0；Cond. No.是条件数，接近 30 以上就要怀疑共线性，这是下一章要专门处理的问题。
回归摘要和论文写作的对应关系，我一般按这张表组织：
| 摘要里的指标 | 论文里的措辞 | 判断标准 |
|---|---|---|
| R-squared | 模型的拟合优度 | 0.6 以上说明解释力可接受，但仍要看样本量 |
| P> | t | 单个变量对因变量的影响是否显著 |
| 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", data=df).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 # 上一次模型的 AIC，None 代表空模型 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, dist="norm", plot=axes[0]) axes[0].set_title("残差 Q-Q 图") # 残差对拟合值散点：观察是否存在单调的喇叭口或弯曲 axes[1].scatter(model.fittedvalues, model.resid, alpha=0.6) axes[1].axhline(0, color="gray", 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_type="HC1") print(model_robust.summary())cov_type="HC1"对应 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(rows=coef_df.shape[0] + 1, cols=4) 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显著性星号在论文里要统一写在表注中：** p<0.01，* p<0.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)这个回读函数能识别变量名并跳过表头和非数字行，所以把报告里任何一张系数表都过一遍，偏差一目了然。比赛前一天把新模型重跑一遍，让报告完整再生一次，再执行这个校验，基本能堵住“模型更新了、文档忘了同步”的低级错误，这也是把回归做完整之后最值得花时间的环节。
本文还有配套的精品资源，点击获取
