线性回归实战：从长江水质预测案例掌握建模全流程与Python实现

1. 从“长江水质”到“线性回归”：一个经典赛题的启示

2005年的全国大学生数学建模竞赛A题，关于长江水质的评价与预测，至今仍是许多建模新手入门的“第一课”。这道题之所以经典，不仅仅是因为它贴近现实，更因为它完美地串联起了数据处理、模型建立、结果分析这一整套建模流程。而其中，线性回归模型作为预测水质指标（如高锰酸盐指数、氨氮含量）随时间或上游污染源变化的工具，扮演了至关重要的角色。很多同学初次接触时，往往直接套用公式，却对模型背后的假设、适用条件以及如何从一堆数据中“说服”模型为自己工作感到迷茫。结果就是，论文里的回归方程看起来像模像样，但一深究“为什么用这个变量？”“模型可靠吗？”，就露了怯。

这篇文章，我们就以这个经典案例为背景，抛开那些枯燥的数学推导教科书，从一个实际建模者的视角，带你重新走一遍用线性回归解决预测问题的完整路径。你会发现，线性回归远不止是y = kx + b那么简单。它涉及到你如何看待数据、如何与数据“对话”、以及如何向评委证明你的模型不是“瞎猜”。我们将重点放在模型的思想、实现的关键步骤、以及那些容易踩坑的细节上，并附上可运行的Python代码（基于scikit-learn和statsmodels）。无论你是正在备赛的数学建模新手，还是希望巩固数据分析基础的同学，这篇文章都能让你对线性回归有一个“既知其然，更知其所以然”的透彻理解。

2. 问题重述与数据理解：建模的第一步不是写公式

在2005年的赛题中，我们需要处理长江流域多个观测站近年的水质数据。假设我们现在聚焦于一个具体问题：预测某个断面未来一段时间的高锰酸盐指数（CODMn）。这是一个典型的时间序列预测问题，但我们可以将其转化为回归问题来处理。

2.1 数据特征工程：从原始表格到模型“食材”

原始数据通常是一张Excel表格，包含“时间”、“观测站”、“CODMn”、“氨氮”、“pH值”等列。直接把这些扔进模型是行不通的。我们需要进行特征工程，把原始数据加工成模型能“消化”的形式。

首先，明确预测目标（因变量y）：我们预测的是“CODMn浓度”。这是一个连续数值，符合线性回归的基本要求。

其次，构造或选择特征（自变量X）：这是最关键的一步，决定了模型的洞察力。我们可以从多个角度构造特征：

1. 时间趋势特征：这是最直接的。将“年月”转化为数值型变量，例如“距离起始月份的月数”。这可以捕捉水质随时间的长期缓慢变化趋势（如治理效果显现）。
2. 周期性特征：水质受季节影响显著。我们可以创建“月份”的循环编码（sin/cos），或者直接使用月份（1-12）作为类别特征（需独热编码）。例如，month_sin = np.sin(2 * np.pi * month/12)，month_cos = np.cos(2 * np.pi * month/12)。这样模型就能学到夏季丰水期稀释作用与冬季枯水期污染物浓度升高的规律。
3. 滞后特征：过去的水质会影响现在。我们可以加入前1个月、前2个月甚至前12个月的CODMn值作为特征。这实际上引入了自回归的思想，对于时间序列预测非常有效。
4. 上游站点特征：如果数据包含多个站点的信息，那么上游站点的水质指标（滞后一期）可以作为下游站点预测的强特征。这体现了污染物的输送效应。
5. 交互特征：例如“时间趋势 * 季节”，可以捕捉季节效应随时间的变化（比如治理措施可能削弱了季节影响的强度）。

注意：特征不是越多越好。过多的特征，尤其是高度相关的特征，会导致“多重共线性”问题，使模型系数估计不稳定，难以解释。我们后续的步骤会处理这个问题。

2.2 数据清洗与探索性分析（EDA）：看清数据的“脾气”

在建模前，必须花时间了解你的数据。这步偷懒，后面模型出问题都找不到北。

1. 处理缺失值：水质数据常有缺失。对于时间序列，常用前向填充（用前一个时间点的值填充）或线性插值。对于非时间序列特征，可以考虑均值、中位数填充，或直接删除缺失过多的样本。关键是要记录你的处理方法，并在论文中说明。
2. 异常值检测：由于监测误差或特殊事件（如突发污染），数据中可能存在异常值。可以用箱线图或3σ原则（三倍标准差）初步识别。对于异常值，需要结合业务判断：是录入错误（可修正或删除）还是真实发生的极端情况（需保留，但模型可能难以拟合）？
3. 可视化：画出CODMn随时间变化的折线图，观察趋势、周期性和异常点。画出特征与目标值的散点图，初步判断线性关系是否成立。计算特征间的相关系数矩阵，观察是否存在高度相关的特征对。

import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns # 假设 df 是已经加载的包含‘time’, ‘station’, ‘CODMn’等列的数据框 df['time'] = pd.to_datetime(df['time']) # 转换时间格式 df = df.sort_values('time') # 按时间排序 # 1. 创建基础特征 df['months_from_start'] = (df['time'].dt.year - df['time'].dt.year.min()) * 12 + (df['time'].dt.month - df['time'].dt.month.min()) df['month'] = df['time'].dt.month df['month_sin'] = np.sin(2 * np.pi * df['month']/12) df['month_cos'] = np.cos(2 * np.pi * df['month']/12) # 2. 创建滞后特征 (以滞后1期为例) df['CODMn_lag1'] = df['CODMn'].shift(1) # 3. 简单处理缺失值（滞后特征首行会产生NaN） df = df.dropna(subset=['CODMn_lag1']).copy() # 简单删除，也可用插值 # 4. 探索性可视化 fig, axes = plt.subplots(2, 2, figsize=(12, 8)) # 时间序列图 axes[0, 0].plot(df['time'], df['CODMn'], marker='o', markersize=3) axes[0, 0].set_title('CODMn Time Series') axes[0, 0].set_xlabel('Time') axes[0, 0].set_ylabel('CODMn') axes[0, 0].grid(True, linestyle='--', alpha=0.7) # 与时间趋势的散点图 axes[0, 1].scatter(df['months_from_start'], df['CODMn'], alpha=0.6) axes[0, 1].set_title('CODMn vs. Time Trend') axes[0, 1].set_xlabel('Months from Start') axes[0, 1].set_ylabel('CODMn') # 可以尝试画一条趋势线 z = np.polyfit(df['months_from_start'], df['CODMn'], 1) p = np.poly1d(z) axes[0, 1].plot(df['months_from_start'], p(df['months_from_start']), "r--", label=f'Trend: y={z[0]:.4f}x+{z[1]:.2f}') axes[0, 1].legend() # 与滞后值的散点图 axes[1, 0].scatter(df['CODMn_lag1'], df['CODMn'], alpha=0.6) axes[1, 0].set_title('CODMn vs. Its Lag-1 Value') axes[1, 0].set_xlabel('CODMn (t-1)') axes[1, 0].set_ylabel('CODMn (t)') axes[1, 0].plot([df['CODMn_lag1'].min(), df['CODMn_lag1'].max()], [df['CODMn_lag1'].min(), df['CODMn_lag1'].max()], "r--", label='y=x line') axes[1, 0].legend() # 月份箱线图 month_order = range(1, 13) boxplot_data = [df[df['month']==m]['CODMn'].dropna() for m in month_order] axes[1, 1].boxplot(boxplot_data, labels=month_order) axes[1, 1].set_title('CODMn Distribution by Month') axes[1, 1].set_xlabel('Month') axes[1, 1].set_ylabel('CODMn') plt.tight_layout() plt.show() # 5. 相关系数矩阵 features_for_corr = ['CODMn', 'months_from_start', 'month_sin', 'month_cos', 'CODMn_lag1'] corr_matrix = df[features_for_corr].corr() plt.figure(figsize=(8,6)) sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0, square=True) plt.title('Feature Correlation Matrix') plt.show()

通过这组图，你能直观看到数据是否存在明显的下降（治理有效）或上升趋势，季节性波动是否明显，当前值与前期值关系是否紧密，以及不同月份水质的中位数和离散程度。相关系数图则能警告你，month_sin和month_cos这类构造特征之间是否独立，以及CODMn_lag1与目标值CODMn的相关性有多强（通常很强，是好特征）。

3. 模型建立、训练与评估：不只是调用fit()

数据准备好了，终于可以建模了。但请记住，建立模型是一个循环迭代的过程，而不是一蹴而就。

3.1 模型选择与假设检验：理解你手中的“武器”

我们选择普通最小二乘（OLS）线性回归。它的核心是找到一组系数，使得预测值与真实值之差的平方和最小。但在使用前，必须心里清楚它的四大基本假设：

1. 线性关系：自变量与因变量之间存在线性关系。
2. 独立性：观测值之间相互独立（对于时间序列数据，这是一个强假设，通常不成立，残差可能存在自相关）。
3. 同方差性：残差的方差在所有观测点上应保持恒定。
4. 正态性：残差应近似服从正态分布。

我们的EDA部分已经初步检验了假设1（通过散点图）。假设2对于时间序列通常是挑战。假设3和4需要在模型拟合后，通过分析残差来检验。

3.2 训练-测试集划分：避免“自欺欺人”

绝对不能用全部数据来训练和评估模型！对于时间序列数据，不能随机划分，必须按时间顺序划分，以模拟真实的预测场景。

from sklearn.model_selection import TimeSeriesSplit from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score # 准备特征矩阵X和目标向量y feature_cols = ['months_from_start', 'month_sin', 'month_cos', 'CODMn_lag1'] X = df[feature_cols].values y = df['CODMn'].values # 按时间顺序划分：前80%训练，后20%测试 split_idx = int(len(X) * 0.8) X_train, X_test = X[:split_idx], X[split_idx:] y_train, y_test = y[:split_idx], y[split_idx:] time_train, time_test = df['time'].iloc[:split_idx], df['time'].iloc[split_idx:] print(f"训练集样本数: {X_train.shape[0]}, 测试集样本数: {X_test.shape[0]}") # 创建并训练模型 model = LinearRegression() model.fit(X_train, y_train) # 在训练集和测试集上进行预测 y_train_pred = model.predict(X_train) y_test_pred = model.predict(X_test) # 评估指标 def print_metrics(y_true, y_pred, set_name): mse = mean_squared_error(y_true, y_pred) rmse = np.sqrt(mse) mae = mean_absolute_error(y_true, y_pred) r2 = r2_score(y_true, y_pred) print(f"{set_name}集评估:") print(f" MSE: {mse:.4f}") print(f" RMSE: {rmse:.4f}") # 与目标值同量纲，更易解释 print(f" MAE: {mae:.4f}") print(f" R²: {r2:.4f}") return rmse, mae, r2 rmse_train, mae_train, r2_train = print_metrics(y_train, y_train_pred, "训练") rmse_test, mae_test, r2_test = print_metrics(y_test, y_test_pred, "测试")

这里的关键是观察训练集和测试集性能的对比。如果训练集R²很高（如>0.9），但测试集R²很低（如<0.5），甚至为负，说明模型过拟合了，它在死记硬背训练数据的噪声，而没有学到普适规律。如果两者都低，则可能是欠拟合，模型太简单（特征不够或关系非线性的）。

3.3 模型诊断：深入残差，发现隐藏问题

得到预测值和评估指标只是第一步。一个负责任的建模者必须诊断模型残差，检验之前提到的基本假设。

# 计算训练集残差 residuals_train = y_train - y_train_pred fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # 1. 残差 vs. 拟合值图 (检验同方差性) axes[0, 0].scatter(y_train_pred, residuals_train, alpha=0.6) axes[0, 0].axhline(y=0, color='r', linestyle='--') axes[0, 0].set_xlabel('Fitted Values (Predicted CODMn)') axes[0, 0].set_ylabel('Residuals') axes[0, 0].set_title('Residuals vs. Fitted Values') axes[0, 0].grid(True, linestyle='--', alpha=0.7) # 理想情况：残差随机均匀分布在0线上下，无明显模式（如漏斗形、弧形）。 # 2. 残差的正态Q-Q图 (检验正态性) import scipy.stats as stats stats.probplot(residuals_train, dist="norm", plot=axes[0, 1]) axes[0, 1].set_title('Normal Q-Q Plot of Residuals') axes[0, 1].grid(True, linestyle='--', alpha=0.7) # 理想情况：点大致落在对角线上。 # 3. 残差 vs. 时间序列图 (检验独立性与时间相关结构) axes[1, 0].plot(time_train, residuals_train, marker='o', markersize=3) axes[1, 0].axhline(y=0, color='r', linestyle='--') axes[1, 0].set_xlabel('Time') axes[1, 0].set_ylabel('Residuals') axes[1, 0].set_title('Residuals vs. Time Order') axes[1, 0].grid(True, linestyle='--', alpha=0.7) # 理想情况：残差在0线上下随机波动，无明显的趋势或周期性。 # 4. 残差自相关图 (ACF) - 严格检验时间序列独立性 from statsmodels.graphics.tsaplots import plot_acf plot_acf(residuals_train, lags=20, ax=axes[1, 1], title='Autocorrelation of Residuals') axes[1, 1].set_ylim(-0.5, 1.1) # 理想情况：除了0阶自相关为1，其他阶的自相关系数均落在置信区间内（蓝色阴影区域），说明残差是白噪声。 plt.tight_layout() plt.show()

如何解读这些诊断图？

• 残差vs拟合值：如果残差随拟合值增大而扩散（漏斗形），说明存在异方差性，可能需要对因变量做变换（如取对数）。
• Q-Q图：如果点严重偏离对角线，尤其是两端偏离，说明残差非正态。对于大样本量，线性回归对正态性假设有一定鲁棒性，但严重偏离可能影响显著性检验。
• 残差vs时间和自相关图：这是时间序列建模的重灾区！如果残差图显示出明显的趋势或周期性，或者自相关图在滞后1、2期等位置超出置信区间，则强烈表明残差存在自相关。这意味着模型没有捕捉到数据中的所有时间依赖结构，违反了独立性假设。这是05年水质预测问题中非常常见的情况，因为水质变化具有连续性和记忆性。

3.4 处理自相关：引入ARIMA思想或更复杂的模型

当诊断出残差自相关时，我们不能视而不见。有几种处理思路：

1. 增加滞后特征：我们已经做了（CODMn_lag1），可以尝试增加更多期滞后（lag2,lag3...），这相当于在回归模型中引入了自回归项。
2. 对残差建模：拟合一个ARIMA模型来捕捉线性回归残差中的时间序列结构，形成回归-ARIMA组合模型。这在statsmodels中可以实现，但复杂度较高。
3. 使用专门的时间序列模型：如ARIMA、SARIMA（带季节性的）直接对原始序列建模。这通常更纯粹，但可能丢失一些我们构造的外部特征（如月份循环编码）的解释能力。
4. 使用树模型：如随机森林、梯度提升树（如XGBoost, LightGBM）。这些模型对特征间的复杂关系和交互作用捕捉能力更强，且对残差的自相关不那么敏感，在时间序列预测中表现往往优于简单线性回归。在当今的数学建模竞赛中，这已经是主流且高效的选择。

这里，我们演示第一种思路的增强版，并引入特征选择。

from sklearn.feature_selection import RFE # 递归特征消除 from sklearn.preprocessing import StandardScaler from statsmodels.stats.outliers_influence import variance_inflation_factor # 创建更多滞后特征 lags = [1, 2, 3, 12] # 滞后1,2,3个月和1年 for lag in lags: df[f'CODMn_lag{lag}'] = df['CODMn'].shift(lag) # 删除因创建滞后特征产生的NaN行 df_expanded = df.dropna().copy() # 重新定义特征和目标 expanded_feature_cols = ['months_from_start', 'month_sin', 'month_cos'] + [f'CODMn_lag{lag}' for lag in lags] X_exp = df_expanded[expanded_feature_cols] y_exp = df_expanded['CODMn'] # 划分训练测试集 (注意索引对齐) split_idx_exp = int(len(X_exp) * 0.8) X_train_exp, X_test_exp = X_exp.iloc[:split_idx_exp], X_exp.iloc[split_idx_exp:] y_train_exp, y_test_exp = y_exp.iloc[:split_idx_exp], y_exp.iloc[split_idx_exp:] # 标准化特征（对于某些特征选择和正则化方法很重要） scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train_exp) X_test_scaled = scaler.transform(X_test_exp) # 方法1：使用RFE进行特征选择 estimator = LinearRegression() selector = RFE(estimator, n_features_to_select=5, step=1) # 选择5个最重要的特征 selector = selector.fit(X_train_scaled, y_train_exp) selected_features_mask = selector.support_ selected_features = X_train_exp.columns[selected_features_mask] print("RFE选出的特征:", list(selected_features)) # 用选出的特征重新训练模型 model_rfe = LinearRegression() model_rfe.fit(X_train_scaled[:, selected_features_mask], y_train_exp) y_train_pred_rfe = model_rfe.predict(X_train_scaled[:, selected_features_mask]) y_test_pred_rfe = model_rfe.predict(X_test_scaled[:, selected_features_mask]) print_metrics(y_train_exp, y_train_pred_rfe, "RFE-训练") print_metrics(y_test_exp, y_test_pred_rfe, "RFE-测试") # 方法2：检查多重共线性 - 方差膨胀因子(VIF) # 注意：VIF计算需要常数项，且对原始数据（未标准化）计算更有意义。 from statsmodels.tools.tools import add_constant X_for_vif = add_constant(X_train_exp) # 添加常数项 vif_data = pd.DataFrame() vif_data["feature"] = X_for_vif.columns vif_data["VIF"] = [variance_inflation_factor(X_for_vif.values, i) for i in range(X_for_vif.shape[1])] print("\n方差膨胀因子(VIF):") print(vif_data) # 通常，VIF > 10 表示存在严重的多重共线性，需要考虑删除或合并特征。

通过RFE，我们可以让模型自动筛选出对预测目标贡献最大的特征，避免冗余。通过VIF，我们可以量化特征间的共线性。如果months_from_start和某个滞后项的VIF很高，说明它们信息重叠，可能需要只保留一个。

4. 模型解释、结果可视化与报告撰写

模型建好了，指标也看了，最后一步是如何将你的工作清晰、有说服力地呈现出来。这是数学建模论文拿高分的关键。

4.1 解释模型系数：赋予数字以物理意义

线性回归的一大优势是可解释性。每个特征前的系数，代表了在其他特征不变的情况下，该特征每增加一个单位，预测的CODMn浓度平均变化多少。

# 获取最终模型的系数和截距 final_model = model_rfe # 假设我们使用RFE筛选后的模型 final_features = selected_features final_coef = final_model.coef_ final_intercept = final_model.intercept_ print("最终线性回归模型方程：") equation = f"预测CODMn = {final_intercept:.4f}" for feat, coef in zip(final_features, final_coef): equation += f" + ({coef:.4f}) * {feat}" print(equation + "\n") print("特征重要性（基于系数绝对值，需注意特征量纲已标准化）：") importance_df = pd.DataFrame({ 'Feature': final_features, 'Coefficient': final_coef, 'Abs_Coefficient': np.abs(final_coef) }).sort_values('Abs_Coefficient', ascending=False) print(importance_df)

解读示例：

• CODMn_lag1的系数为0.65（假设）：这意味着，在其他条件不变的情况下，上一个月的CODMn每升高1个单位，本月的CODMn平均会升高0.65个单位。这符合水质具有连续性的认知。
• month_sin的系数为负：这可能意味着在正弦波对应的某个相位（如夏季），CODMn浓度较低，反映了季节性稀释效应。
• months_from_start的系数为负且很小：这可能暗示了一个长期缓慢的下降趋势，可能与多年的治理投入有关。

注意：由于我们之前对特征进行了标准化，此时的系数大小可以直接比较，用于判断特征重要性。但解释“每增加一个单位”时，要明白这是指标准化后的一个标准差单位。如果想得到原始单位的解释，需要用未标准化的数据训练模型，或者进行系数转换。

4.2 可视化预测结果：让评委一目了然

一张好的结果图，胜过千言万语。

# 将训练集和测试集的预测结果合并回时间序列 df_expanded['prediction'] = np.nan train_idx = df_expanded.index[:split_idx_exp] test_idx = df_expanded.index[split_idx_exp:] df_expanded.loc[train_idx, 'prediction'] = y_train_pred_rfe df_expanded.loc[test_idx, 'prediction'] = y_test_pred_rfe plt.figure(figsize=(14, 7)) # 绘制真实值 plt.plot(df_expanded['time'], df_expanded['CODMn'], 'b-', label='Observed CODMn', linewidth=1.5, alpha=0.8) # 绘制预测值 plt.plot(df_expanded['time'], df_expanded['prediction'], 'r--', label='Predicted CODMn', linewidth=2) # 标记训练集和测试集分界线 split_time = df_expanded['time'].iloc[split_idx_exp] plt.axvline(x=split_time, color='g', linestyle=':', linewidth=2, label='Train/Test Split') plt.fill_betweenx(y=[df_expanded['CODMn'].min(), df_expanded['CODMn'].max()], x1=df_expanded['time'].min(), x2=split_time, color='gray', alpha=0.1, label='Training Period') plt.fill_betweenx(y=[df_expanded['CODMn'].min(), df_expanded['CODMn'].max()], x1=split_time, x2=df_expanded['time'].max(), color='yellow', alpha=0.1, label='Testing Period') plt.xlabel('Time') plt.ylabel('CODMn Concentration') plt.title('Linear Regression Model: Observed vs Predicted CODMn') plt.legend(loc='best') plt.grid(True, linestyle='--', alpha=0.5) plt.tight_layout() plt.show() # 单独绘制测试集的预测效果 plt.figure(figsize=(10, 6)) plt.scatter(y_test_exp, y_test_pred_rfe, alpha=0.6, edgecolors='k') # 绘制理想预测线 y=x min_val = min(y_test_exp.min(), y_test_pred_rfe.min()) max_val = max(y_test_exp.max(), y_test_pred_rfe.max()) plt.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Prediction (y=x)') plt.xlabel('True CODMn (Test Set)') plt.ylabel('Predicted CODMn (Test Set)') plt.title('Test Set: True vs Predicted Values') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) # 在图上标注R²和RMSE plt.text(0.05*max_val, 0.9*max_val, f'$R^2$ = {r2_test:.3f}\nRMSE = {rmse_test:.3f}', bbox=dict(boxstyle='round,pad=0.5', facecolor='wheat', alpha=0.8)) plt.tight_layout() plt.show()

第一张时序对比图，能清晰展示模型在整个时间轴上的拟合效果，以及在测试集（未来时段）上的预测能力。评委最关注的就是模型在“未见过的数据”上表现如何。第二张散点图则能直观看出预测值与真实值的偏离程度，点越靠近红色对角线，预测越准。

4.3 撰写建模报告的核心要点

在论文中，关于线性回归模型部分，你需要清晰地阐述以下内容，这体现了你的建模素养：

1. 变量选择依据：你为什么选择这些特征？是基于物理意义（如滞后项代表持续性）、统计检验（如相关性分析）还是领域知识（如季节性）？
2. 模型建立过程：你使用了什么方法处理多重共线性？（如VIF分析、特征选择）。你如何划分训练集和测试集？为什么这样划分？
3. 模型检验：不要只给出R²和RMSE。必须展示残差分析图（特别是残差vs时间和自相关图），并讨论是否满足线性回归假设。如果存在自相关，你采取了什么措施？（如增加滞后特征，或说明这是本模型的局限性，建议后续采用时间序列模型）。
4. 模型结果与解释：给出最终的回归方程，并解释关键系数的实际意义。例如：“滞后一期系数为0.82，表明长江水质具有强烈的自相关性，上月污染状况对本月有显著影响。”
5. 预测与评估：展示在测试集上的预测效果图，并给出具体的评估指标。客观分析模型的优缺点。例如：“模型能较好地捕捉水质的季节趋势和短期依赖，但在突变点（如突发污染事件）预测能力不足，因为模型未引入降雨量、排污量等外部冲击变量。”

5. 进阶思考与常见陷阱

走通了整个流程，我们再来探讨几个更深层次的问题和比赛中容易踩的坑。

5.1 线性回归的“非线性”扩展

水质变化与某些因素可能并非简单的线性关系。例如，污染物浓度与流量可能呈负指数关系。此时，我们可以通过变量变换将非线性关系线性化。

• 多项式特征：引入months_from_start的平方项、立方项，可以拟合更复杂的时间趋势。
• 交互项：引入month_sin * months_from_start，可以检验季节效应是否随时间变化。
• 对数变换：如果因变量CODMn呈现指数增长或衰减趋势，可以尝试对CODMn取对数，建立log(CODMn)与自变量的线性模型。注意：这要求CODMn值全为正，且解释系数时需要反变换。

在scikit-learn中，可以使用PolynomialFeatures来自动生成多项式特征和交互项，但要警惕特征维度爆炸和过拟合。

5.2 过拟合与正则化

当我们引入大量特征（如多个滞后项、多项式项）时，模型很容易过拟合。除了使用RFE等特征选择方法，正则化是控制过拟合的利器。

• 岭回归（Ridge）：在损失函数中加入L2正则项，惩罚过大的系数，使模型更稳定。
• Lasso回归（Lasso）：在损失函数中加入L1正则项，它可以将不重要的特征的系数直接压缩为0，从而实现自动特征选择。

from sklearn.linear_model import Ridge, LassoCV from sklearn.preprocessing import StandardScaler # 使用LassoCV自动交叉验证选择最佳的正则化强度alpha lasso_cv = LassoCV(cv=5, random_state=42, max_iter=10000).fit(X_train_scaled, y_train_exp) print(f"Lasso选出的最佳alpha: {lasso_cv.alpha_:.6f}") print(f"Lasso模型非零系数个数: {np.sum(lasso_cv.coef_ != 0)}") # 比较模型性能 models = { 'OLS': LinearRegression(), 'Ridge (alpha=1.0)': Ridge(alpha=1.0), 'Lasso (CV)': lasso_cv } for name, model in models.items(): if name != 'Lasso (CV)': model.fit(X_train_scaled, y_train_exp) y_pred_test = model.predict(X_test_scaled) r2 = r2_score(y_test_exp, y_pred_test) print(f"{name:15} - 测试集 R²: {r2:.4f}")

在论文中，如果你使用了正则化，需要说明原因（防止过拟合）和选择正则化参数的方法（如交叉验证）。

5.3 数学建模中的“软技巧”

1. 结果稳健性检验：改变训练测试集划分比例（如70/30, 85/15），看模型性能是否发生剧烈变化。使用时间序列交叉验证（TimeSeriesSplit）来获得更稳健的性能估计。
2. 多模型对比：不要只提交一个线性回归模型。可以尝试建立多个模型（如线性回归、带滞后项的线性回归、ARIMA、随机森林），在测试集上比较它们的RMSE或MAE。在论文中展示一个简单的模型对比表格，并说明你最终选择某个模型的理由（可能是精度最高，也可能是可解释性最好）。
3. 敏感性分析：探讨关键参数（如滞后阶数、正则化强度）对模型结果的影响。这能体现你对模型的理解深度。
4. 承认局限性：没有完美的模型。明确指出你的模型假设（如线性、无自相关）在哪些情况下可能不成立，以及模型未考虑哪些重要因素（如政策突变、极端气候）。提出改进方向，这往往是论文的加分项。

线性回归是数学建模中最基础、最常用的工具之一。通过长江水质这个具体案例，我们从数据预处理、特征工程、模型建立、诊断检验到结果解释，完整地走了一遍实战流程。希望这篇文章能让你明白，建立一个可靠的回归模型，其核心不在于复杂的数学，而在于严谨的数据思维、系统的检验流程和清晰的逻辑表述。下次当你再看到“预测”、“关联分析”这类关键词时，希望你能自信地拿起线性回归这个工具，并且知道如何用它做出一个经得起推敲的答案。