多元线性回归（Multiple Linear Regression, MLR）是一种统计模型，被广泛认为是计量经济学的核心基础。多元线性回归为经济研究者提供了一种有效的方法来建模和分析多个自变量与因变量之间的线性关系。
在计量经济学中，研究者常常面临复杂的经济现象，这些现象往往受多种因素影响。通过建立多元线性回归模型，研究者能够定量分析各个自变量对因变量的影响程度，从而揭示经济关系的本质。例如，在分析工资水平时，研究者可能会考虑教育、工作经验、行业和地区等多个因素，通过多元线性回归，能够有效识别这些因素的边际贡献。
多元线性回归不仅可以用于模型构建，还可以用于假设检验，检验经济理论的有效性。计量经济工作者可以使用回归分析来验证某一理论预测是否成立，如需求弹性或生产函数的形式。通过统计检验（如t检验和F检验），研究者可以判断自变量对因变量的影响是否显著，进而为经济理论提供实证支持。
多元线性回归具有良好的统计性质，如无偏性和一致性，使得其在经济数据分析中成为一种稳健的工具。随着数据分析技术的发展，多元线性回归也可以与现代机器学习方法结合，扩展其应用范围。通过对模型进行适当的调整和扩展，研究者能够应对非线性关系、多重共线性、异方差性等问题，从而提高模型的准确性和适用性。
一、多元线性回归模型
在面对复杂的非线性系统时，将其转化为线性系统进行处理的原因和考量主要集中在以下几个方面。
简化计算与模型求解：线性模型的计算相对简单，参数估计通常使用最小二乘法，可以直接求解。而对于非线性模型，通常需要迭代算法（如牛顿法、梯度下降法）来估计参数，这些方法计算复杂且容易陷入局部最优解。将非线性问题通过适当的变换处理成线性问题，可以利用现成的线性回归方法，更快速、稳定地得到解。
稳健性与理论基础：线性回归有着良好的统计学性质，如无偏性、有效性和一致性。这些性质使得线性模型在估计和预测中表现得更加稳健。相比之下，非线性模型因为更复杂，容易受到过拟合、参数估计不准确等问题的影响。通过线性化处理，研究者可以保留这些稳健性特征，减少模型不稳定的风险。
易于解释和理解：线性模型的结果易于解释，回归系数表示每个自变量对因变量的边际贡献，直观而且容易被理解。而非线性模型通常难以直接解释自变量与因变量之间的关系。将非线性问题转化为线性模型，可以让研究者和读者更清楚地理解变量之间的相互作用。
1.1 模型假设
多元线性回归模型的基本形式为：
其中：
- $ Y$ 是因变量，表示我们希望解释或预测的结果。
- $ X_1, X_2, \cdots, X_k$ 是自变量，表示影响因变量的因素。
- $ \beta_0$ 是截距项，$ \beta_1, \beta_2, \cdots, \beta_k$ 是各个自变量的回归系数。
- $ \varepsilon$ 是随机误差项，表示无法通过自变量解释的部分。
为了保证回归模型能够给出有效的估计结果，必须满足以下几个经典假设：
- 线性关系：因变量与自变量之间的关系是线性的，即因变量是自变量的线性组合。
- 独立性：误差项彼此之间是相互独立的，不存在相关性。引出自相关等
- 同方差性：误差项的方差是恒定的，不随自变量的取值变化而变化。引出异方差性
- 正态性：误差项服从正态分布，尤其在样本量较小时这一假设尤为重要。引出非正态Logit回归等
- 无共线性：自变量之间没有完全的线性相关性，即不存在多重共线性问题。引出多重共线性
1.2 参数估计
多元线性回归的参数估计方法通常采用普通最小二乘法（Ordinary Least Squares, OLS），这是最经典的估计方法之一，其目标是最小化回归模型中的残差平方和，即预测值与实际观测值之间差距的平方和最小化。OLS 在经典线性回归模型的框架下表现出优良的统计性质，如无偏性、最小方差性和一致性，使得它成为经济学、社会科学等领域分析变量关系的重要工具。
普通最小二乘法的原理
在多元线性回归模型中，我们假设因变量 $ Y $ 与一组自变量 $ X $ 之间的关系可以用以下线性方程表示：
其中：
- $ Y $ 是因变量的向量；
- $ X $ 是 $ n \times k $ 维的设计矩阵，包含 $ n $ 个样本和 $ k $ 个自变量；
- $ \beta $ 是回归系数向量；
- $ \varepsilon $ 是误差项向量。
普通最小二乘法的核心思想是通过最小化模型预测值 $ X \beta $ 与实际观测值 $ Y $ 之间的差距，即最小化残差平方和来估计回归系数 $ \beta $。残差定义为 $ \varepsilon = Y - X \beta $，因此残差平方和 $ S(\beta) $ 为：
最小化这个残差平方和的过程可以通过求解偏导数得到最优的 $ \beta $。该优化问题的解是：
其中，$ \hat{\beta} $ 是估计的回归系数向量，$ X^T X $ 是设计矩阵 $ X $ 的转置矩阵，$ (X^T X)^{-1} $ 是矩阵 $ X^T X $ 的逆矩阵。这个方程是 OLS 的闭式解，用于计算最佳的回归系数。
普通最小二乘法的性质
在经典线性回归模型的假设条件下，OLS 估计具有以下三个主要性质：
- 无偏性：OLS 估计是无偏的，即在重复抽样的条件下，OLS 的估计值 $ \hat{\beta} $ 的期望等于真实值 $ \beta $。这意味着，尽管单个样本的估计可能有误差，但在大量重复样本中，估计结果平均来说是准确的。
- 最小方差性：根据高斯-马尔科夫定理，在所有线性无偏估计量中，OLS 具有最小的方差。这意味着在满足模型假设的条件下，OLS 是最有效的线性估计量，它能够提供最小的估计误差。
- 一致性：随着样本量的增加，OLS 估计量趋近于真实参数值，即估计是“收敛”的。无论样本量如何大小，OLS 估计量在大样本情况下会更接近真实的参数。
这些性质使 OLS 成为统计和计量经济学中广泛应用的估计方法。
1.3 常见模型问题及其解决方法
异方差性
异方差性是指模型中误差项的方差不是恒定的，即残差的波动随着自变量的变化而变化。这违反了经典线性回归模型的同方差性假设，导致回归系数的标准误估计值失效，进而影响t检验和F检验的准确性。
异方差性的检验
- 图形法：通过绘制残差与拟合值的散点图，观察是否存在系统性的模式。如果残差图中出现漏斗形或其他形态，则可能存在异方差性。
- 布鲁斯-帕根检验（Breusch-Pagan Test）：该检验通过检验残差平方是否与自变量存在关联，来判断是否存在异方差性。
- 怀特检验（White Test）：通过构建自变量及其平方项来检验模型中的异方差性。
消除异方差性的方法
- 加权最小二乘法（WLS）：通过对不同观测值赋予不同的权重，以解决误差项方差不等的问题。加权较大的值代表方差较小的观测，降低其对整体估计的影响。
- 稳健标准误（Robust Standard Errors）：若难以消除异方差性，稳健标准误可以在异方差性存在的情况下给出更为准确的标准误估计，保证显著性检验的有效性。
- 广义最小二乘法（GLS）：GLS 是一种扩展的OLS方法，能够处理误差项的异方差性问题。
多重共线性
多重共线性是指模型中的自变量之间存在较高的线性相关性，使得回归系数无法准确估计。共线性问题不会影响模型的拟合优度，但会导致回归系数的标准误增大，影响自变量的显著性检验，甚至会使某些重要变量被错误地认为是不显著的。
多重共线性的检验
①方差膨胀因子
方差膨胀因子是检测多重共线性的最常用方法之一。它衡量某个自变量 $ X_i $ 与其他所有自变量之间的线性相关性程度，具体公式为：
其中，$ R_i^2 $ 是将自变量 $ X_i $ 对其他自变量进行回归后的决定系数。如果 $ R_i^2 $ 较高，说明 $ X_i $ 与其他自变量之间存在较强的线性关系，VIF 值也会随之增大。
判断依据：
- 当 $ \text{VIF}_i \approx 1 $ 时，表明没有多重共线性；
- 当 $ \text{VIF}_i > 5 $ 时，通常认为存在中等程度的多重共线性；
- 当 $ \text{VIF}_i > 10 $ 时，认为存在严重的多重共线性。
VIF 提供了一个简单易用的工具来判断某个特定变量是否与其他自变量存在共线性问题。如果所有自变量的 VIF 值都较高，说明回归模型整体上可能存在共线性问题。
②特征值分解法（Eigenvalue Decomposition）和条件数（Condition Number）
特征值分解法是通过分析设计矩阵 $ X $ 的特征值来判断多重共线性的方法。特征值表示的是矩阵中的线性独立信息的多少，特征值越小，说明矩阵中的变量越不独立，存在多重共线性。在特征值分解的基础上，条件数（Condition Number）常用于进一步判定共线性。条件数定义为设计矩阵 $ X $ 的最大特征值 $ \lambda_{\text{max}} $ 与最小特征值 $ \lambda_{\text{min}} $ 的比值：
判断依据：
- 当 $ \kappa(X) < 10 $ 时，说明没有显著的多重共线性；
- 当 $ \kappa(X) > 30 $ 时，通常认为存在严重的多重共线性。
特征值分解法是一种更加数学化的方法，适合用于深入分析多重共线性问题，特别是在有多个自变量的情况下。
③相关矩阵法（Correlation Matrix）
相关矩阵法是通过检查自变量之间的相关系数来识别共线性。计算每对自变量 $ X_i $ 和 $ X_j $ 之间的皮尔逊相关系数 $ \rho_{ij}$：
其中，$ \text{Cov}(X_i, X_j) $ 是 $ X_i $ 和 $ X_j $ 之间的协方差，$ \sigma_{X_i} $ 和 $ \sigma_{X_j} $ 分别是自变量的标准差。
判断依据：
- 当 $ |\rho_{ij}| > 0.8 $ 时，通常认为这对自变量之间存在较强的共线性；
- 当 $ |\rho_{ij}| > 0.9 $ 时，说明存在严重的共线性。
相关矩阵法是识别共线性的初步手段，可以帮助判断变量之间的线性关系，但它无法解决变量与其他多个变量之间的线性相关性问题，因此需要与其他方法配合使用。
④条件指数（Condition Index）
条件指数是通过特征值分析进一步量化共线性程度的指标。它基于特征值分解，将特征值的平方根用于衡量：
其中，$ \lambda_i $ 是第 $ i $ 个特征值。条件指数反映了设计矩阵的不同维度之间的线性相关性。
判断依据：
- 当条件指数大于 10 时，表明存在一定程度的共线性；
- 当条件指数大于 30 时，认为存在严重的多重共线性。
条件指数与条件数相似，但它通过对各个特征值进行逐步分解，能够进一步帮助识别出具体的变量之间共线性的程度。
⑤逐步回归法（Stepwise Regression）
在多重共线性情况下，某些自变量的显著性检验可能会受到影响。因此，逐步回归法（Stepwise Regression）可以通过引入或剔除变量来检测和处理共线性问题。逐步回归法逐步将自变量引入回归模型或将不显著的变量移除，并在每一步都检测回归模型的显著性和自变量之间的相关性。虽然逐步回归法不能直接量化共线性，但它可以通过调整模型结构来减弱共线性的影响。此外，逐步回归法可以与 VIF 等方法结合使用，以提高模型的稳定性。
消除多重共线性的方法
删除高度相关的自变量：如果某些自变量之间高度相关，可以通过删除其中一个变量来减少共线性。
主成分回归（PCR）：通过主成分分析将自变量降维，将相关性较强的变量转化为若干无相关的主成分，再用这些主成分进行回归分析。
岭回归（Ridge Regression）：在回归模型中引入正则化惩罚项，通过控制回归系数的大小，减少共线性对模型的影响。
二、案例分析1
考虑一家快递公司送货：X1： 运输里程 X2： 运输次数 Y：总运输时间，数据见下表所示。
| ID | X1: Miles Traveled | X2: Number of Deliveries | Y: Travel Time (Hours) |
|---|---|---|---|
| 1 | 100 | 4 | 9.3 |
| 2 | 50 | 3 | 4.8 |
| 3 | 100 | 4 | 8.9 |
| 4 | 100 | 2 | 6.5 |
| 5 | 50 | 2 | 4.2 |
| 6 | 80 | 2 | 6.2 |
| 7 | 75 | 3 | 7.4 |
| 8 | 65 | 4 | 6.0 |
| 9 | 90 | 3 | 7.6 |
| 10 | 90 | 2 | 6.1 |
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor# 数据输入
data = {'X1': [100, 50, 100, 100, 50, 80, 75, 65, 90, 90], # Miles traveled'X2': [4, 3, 4, 2, 2, 2, 3, 4, 3, 2], # Number of deliveries'Y': [9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6.0, 7.6, 6.1] # Travel time
}# 创建DataFrame
df = pd.DataFrame(data)# 自变量（X1 和 X2）
X = df[['X1', 'X2']]
X = sm.add_constant(X) # 添加常数项# 因变量（Y）
Y = df['Y']# OLS回归
model = sm.OLS(Y, X).fit()# 输出回归结果
summary_table = model.summary2()
print(summary_table)# Breusch-Pagan 异方差检验
bp_test = het_breuschpagan(model.resid, X)
bp_test_results = {'Lagrange multiplier statistic': bp_test[0],'p-value': bp_test[1],'f-value': bp_test[2],'f p-value': bp_test[3]
}
print("\nBreusch-Pagan Test Results:")
for key, value in bp_test_results.items():print(f"{key}: {value:.4f}")# 方差膨胀因子（VIF） 共线性检验
vif = pd.DataFrame()
vif["Variable"] = X.columns
vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]print("\nVariance Inflation Factor (VIF):")
print(vif)
OLS 回归结果：
Dependent Variable: Y
Model: OLS
R-squared: 0.904
Adj. R-squared: 0.876
No. Observations: 10
F-statistic: 32.88
Prob (F-statistic): 0.000276
Log-Likelihood: -6.8398
AIC: 19.6796
BIC: 20.5873
| Variable | Coefficient | Std. Err. | t-statistic | P-value | 95% Confidence Interval |
|---|---|---|---|---|---|
| const | -0.8687 | 0.9515 | -0.9129 | 0.3916 | [-3.1188, 1.3814] |
| X1 | 0.0611 | 0.0099 | 6.1824 | 0.0005 | [0.0378, 0.0845] |
| X2 | 0.9234 | 0.2211 | 4.1763 | 0.0042 | [0.4006, 1.4463] |
Durbin-Watson: 2.515
Condition No.: 435
Breusch-Pagan 异方差检验结果：
- Lagrange Multiplier Statistic: 0.369
- p-value: 0.8315 （无异方差问题）
方差膨胀因子（VIF）：
| Variable | VIF |
|---|---|
| const | 27.5636 |
| X1 | 1.027 |
| X2 | 1.027 |
该模型表明没有显著的异方差和多重共线性问题。
三、案例分析2
这里采用波士顿房价数据集进行多元回归分析，各个特征变量(自变量)含义如下表所示。
| Attribute | Description |
|---|---|
| CRIM | 该区域内犯罪率（每人均犯罪次数） |
| ZN | 该区域内住宅地的占地面积比例 |
| INDUS | 该区域内非零售商业用地比例 |
| CHAS | 是否靠河，1为靠河，0为不靠河 |
| NOX | 一氧化氮浓度（每千万分之一） |
| RM | 该区域内房屋的平均房间数量 |
| AGE | 该区域内房屋的平均房龄 |
| DIS | 该区域内房屋到市中心的加权距离 |
| RAD | 该区域内房屋到径向高速公路的距离指数 |
| TAX | 该区域内每一万美元的财产税率 |
| PTRATIO | 该区域内教师与学生数之比 |
| B | 该区域内黑人所占比例 |
| LSTAT | 该区域内人口中有多少比例属于低收入阶层 |
| PRICE | 自主房屋房价中位数(被解释变量-因变量) |
3.1 描述性分析
# 导入必要的库
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt# 读取波士顿房价数据
data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)# 数据处理
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]]) # 特征
target = raw_df.values[1::2, 2] # 房价（目标变量）# 创建DataFrame
columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']
df = pd.DataFrame(data, columns=columns)
df['PRICE'] = target # 添加房价列# 计算相关矩阵
corr_matrix = df.corr()# 绘制相关矩阵的热力图
plt.figure(figsize=(12, 8))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Heatmap of Boston Housing Dataset')
plt.show()
可以看出RM，lSTAT，INDUS，PTRATIO都是与PRICE的相关性比较大的特征。
3.2 多元回归分析
波士顿房价数据集是20世纪70年代波士顿郊区房价的中位数，统计了当时教区部分的犯罪率、房产税等共计13个输入变量，统计出房价（1个输出变量）。哪些指标对房价的影响较大？通过对波士顿房价数据集的作图和定性分析，可以得出以下结论：
- 房屋的平均房间数 (RM)：对房价的影响较大，系数为正数，说明房间越多，房价越高。
- 距离五个波士顿就业中心的加权距离 (DIS)：对房价的影响较大，系数为负数，说明距离就业中心越远，房价越低。
- 教育水平较高的人口所占比例 (LSTAT)：对房价的影响较大，系数为负数，说明该比例越高，房价越低。
- 住宅的年龄 (AGE)：对房价的影响较大，系数为负数，说明住宅越老，房价越低。
其他自变量对房价的影响较小，甚至为负数，说明这些因素与房价呈反比例关系。为简化计算仅对上面四个变量做回归，需要注意的是，这样的结论仅代表波士顿房价数据集的情况，不一定适用于其他地区。
# 导入必要的库
import pandas as pd
import numpy as np
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt# 读取波士顿房价数据
data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)# 数据处理
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]]) # 特征
target = raw_df.values[1::2, 2] # 房价（目标变量）# 创建DataFrame
columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']
df = pd.DataFrame(data, columns=columns)
df['PRICE'] = target # 添加房价列# 选择变量：RM, DIS, LSTAT, AGE 和 PRICE
X = df[['RM', 'DIS', 'LSTAT', 'AGE']] # 自变量
y = df['PRICE'] # 因变量# 添加常数项（截距项）
X = sm.add_constant(X)# 使用OLS进行回归分析
model = sm.OLS(y, X).fit()# 输出回归结果
print(model.summary())# 绘制残差图
plt.figure(figsize=(10,6))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True, line_kws={'color': 'red'})
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted values')
plt.show()# 绘制预测值与实际值对比图
plt.figure(figsize=(10,6))
plt.scatter(model.fittedvalues, y, alpha=0.5, color='b')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2)
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.title('Predicted vs Actual')
plt.show()
OLS Regression Results
==============================================================================
Dep. Variable: PRICE R-squared: 0.649
Model: OLS Adj. R-squared: 0.646
Method: Least Squares F-statistic: 231.6
Date: Mon, 14 Oct 2024 Prob (F-statistic): 2.09e-112
Time: 14:41:23 Log-Likelihood: -1575.4
No. Observations: 506 AIC: 3161.
Df Residuals: 501 BIC: 3182.
Df Model: 4
Covariance Type: nonrobust
==============================================================================coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.0799 3.437 1.187 0.236 -2.672 10.832
RM 4.9904 0.449 11.126 0.000 4.109 5.872
DIS -0.6588 0.175 -3.767 0.000 -1.002 -0.315
LSTAT -0.6848 0.054 -12.721 0.000 -0.791 -0.579
AGE -0.0255 0.014 -1.771 0.077 -0.054 0.003
==============================================================================
Omnibus: 124.062 Durbin-Watson: 0.869
Prob(Omnibus): 0.000 Jarque-Bera (JB): 355.890
Skew: 1.171 Prob(JB): 5.24e-78
Kurtosis: 6.376 Cond. No. 1.08e+03
==============================================================================
结果解释
- R-squared (R²):
- 值为 0.649，表示模型解释了约 64.9% 的房价变动。这是一个中等的拟合度，表明还有 35.1% 的变动未被模型解释。
- Adj. R-squared (调整后的 R²):
- 值为 0.646，相较于 R²，这个指标考虑了自变量的数量。调整后的 R² 用于评估模型的拟合度，当增加自变量时，调整后的 R² 可能会降低。
- F-statistic:
- 值为 231.6，表明模型整体是显著的（P值为 2.09e-112）。这意味着至少有一个自变量与因变量 PRICE 之间存在显著关系。
- 各个自变量的系数 (coef):
- RM (房间数): 系数为 4.9904，表明房间数每增加一个单位，房价将增加约 4.99 千美元，且 p 值为 0.000，说明显著。
- DIS (与就业中心的距离): 系数为 -0.6588，说明距离每增加一个单位，房价将降低约 0.66 千美元，且 p 值也显著。
- LSTAT (低收入人群比例): 系数为 -0.6848，表明比例每增加一个单位，房价将降低约 0.68 千美元，且 p 值显著。
- AGE (房屋年龄): 系数为 -0.0255，虽然 p 值为 0.077（接近 0.05），但未达到显著水平，表明房屋年龄对房价的影响不明显。
- Omnibus 和 Jarque-Bera 检验:
- 这些检验用于检查残差的正态性。Omnibus 的 p 值为 0.000，表明残差不服从正态分布，Jarque-Bera 检验的 p 值也非常小，这进一步支持了这一结论。
- Durbin-Watson:
- 值为 0.869，接近于 0，表明可能存在自相关。理想情况下，Durbin-Watson 值应在 1.5 至 2.5 之间。
- Condition Number (条件数):
- 值为 1.08e+03，条件数大于 30 时，表明模型中可能存在多重共线性问题。
3.3 检查结果
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_breusch_godfrey
from sklearn.datasets import fetch_california_housing# 加载数据集：使用加州房价数据作为替代
housing = fetch_california_housing(as_frame=True)
df = housing.frame# 选择自变量和因变量
X = df[['MedInc', 'HouseAge', 'AveRooms', 'AveOccup']]
y = df['MedHouseVal']# 添加常数项
X = sm.add_constant(X)# 拟合 OLS 回归模型
model = sm.OLS(y, X).fit()# 打印回归结果
print(model.summary())# 1. 共线性检测（VIF）
def calculate_vif(X):vif = pd.DataFrame()vif["Variables"] = X.columnsvif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]return vifvif = calculate_vif(X)
print("\n共线性检测结果 (VIF):")
print(vif)# 2. 异方差性检测 (Breusch-Pagan 检验)
bp_test = het_breuschpagan(model.resid, model.model.exog)
bp_test_results = pd.Series(bp_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n异方差性检测结果 (Breusch-Pagan Test):")
print(bp_test_results)# 3. 自相关性检测 (Breusch-Godfrey 检验)
bg_test = acorr_breusch_godfrey(model, nlags=1)
bg_test_results = pd.Series(bg_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n自相关性检测结果 (Breusch-Godfrey Test):")
print(bg_test_results)
检查共线性:VIF 值超过 10 通常表明存在严重的多重共线性。
检查异方差性:Breusch-Pagan 检验的 p 值如果小于 0.05，表明存在异方差性。
检查自相关:Breusch-Godfrey 检验的 p 值如果小于 0.05，说明存在自相关。
3.4 剔除变量HouseAge后的结果
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_breusch_godfrey
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing# 加载加州房价数据集
housing = fetch_california_housing(as_frame=True)
df = housing.frame# 选择自变量，保留两个变量：MedInc, AveRooms
X = df[['MedInc', 'AveRooms']] # 保留所需变量
y = df['MedHouseVal']# 添加常数项
X = sm.add_constant(X)# 计算共线性（VIF）并输出
def calculate_vif(X):vif = pd.DataFrame()vif["Variables"] = X.columnsvif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]return vifvif = calculate_vif(X)
print("\n共线性检测结果 (VIF):")
print(vif)# 如果VIF检测通过，进行异方差性检测和回归分析
model = sm.OLS(y, X).fit(cov_type='HC3') # 使用robust标准误差解决异方差性问题# 打印回归结果, 仅输出变量的回归系数部分，去掉'const'
coef_df = pd.DataFrame({"Variable": X.columns[1:], # 去掉'const'"Coefficient": model.params[1:], # 仅显示变量的系数"P>|t|": model.pvalues[1:], # 仅显示P值"Std. Err.": model.bse[1:] # 显示标准误差
})print("\n最终回归结果（保留两个变量后，带robust标准误差）：")
print(coef_df)# 自相关性检测（Breusch-Godfrey）
bg_test = acorr_breusch_godfrey(model, nlags=1)
bg_test_results = pd.Series(bg_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n自相关性检测结果 (Breusch-Godfrey Test):")
print(bg_test_results)# 异方差性检测 (Breusch-Pagan)
bp_test = het_breuschpagan(model.resid, model.model.exog)
bp_test_results = pd.Series(bp_test, index=['LM Stat', 'LM p-value', 'F-stat', 'F p-value'])
print("\n异方差性检测结果 (Breusch-Pagan Test):")
print(bp_test_results)# 残差图绘制
residuals = model.resid # 获取残差
fitted_values = model.fittedvalues # 获取预测值plt.figure(figsize=(10, 6))
plt.scatter(fitted_values, residuals, edgecolors='k', facecolors='none')
plt.axhline(y=0, color='r', linestyle='--')
plt.title("Residual Plot")
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.show()# 预测部分 - 使用训练好的模型进行预测
# 在这里使用数据集的前5行作为预测示例
new_data = X.iloc[:5, :]
predictions = model.predict(new_data)
print("\n预测结果：")
print(predictions)
从回归结果来看，虽然系数显著且共线性问题不明显，但自相关性问题可能影响了模型的稳健性。序列自相关意味着残差项并非相互独立，这会导致模型的系数估计不再是最优的，并影响模型预测的可靠性。共线性问题：基本不存在；系数显著性：两个自变量均对房价有显著影响，但 AveRooms 的负系数与预期不符，可能需要更复杂的模型来解释；自相关问题：存在显著自相关，影响了回归模型的有效性和预测能力，残差图如下图。
总结
多元线性回归的核心在于通过多个自变量的线性组合来预测因变量的变化，但模型的有效性依赖于假设的满足。当假设未被满足时，估计的回归系数可能不再准确，并且显著性检验结果会失效。异方差性和多重共线性是多元线性回归中常见的两个问题，通过适当的检验和修正方法，可以提高模型的解释能力和预测效果。此外，研究者应结合实际问题，选择适合的模型和方法来处理回归分析中的复杂性。多元线性回归虽然是一种经典的统计方法，但其应用必须谨慎。通过合理的假设检验和对潜在问题的处理，可以确保回归结果的可靠性，并为决策和政策分析提供有力的依据。
参考资料
- 7.3 多元回归分析(multiple regression)
- 关于波士顿房价数据集的回归预测
