Python的统计建模：利用`Statsmodels`进行回归分析和假设检验。 - 智猿学院-前后端，数据库，人工智能，云计算等领域前沿技术讲座

Python统计建模：利用`Statsmodels`进行回归分析和假设检验

大家好，今天我们来探讨Python中的统计建模，重点关注Statsmodels库在回归分析和假设检验中的应用。Statsmodels是一个强大的Python库，提供了丰富的统计模型、推断和评估工具，使得我们能够轻松地进行回归分析、时间序列分析等。

1. Statsmodels简介

Statsmodels是一个开源的Python库，旨在为用户提供统计建模和计量经济学工具。它构建在NumPy和SciPy之上，并与Pandas紧密集成，能够处理各种类型的数据。Statsmodels提供了一系列模型，包括线性回归、广义线性模型、时间序列模型等，并提供了模型诊断、假设检验和结果可视化等功能。

2. 线性回归分析

线性回归是统计建模中最基本也是最重要的模型之一。它试图建立自变量和因变量之间的线性关系。Statsmodels提供了多种线性回归模型，其中最常用的是普通最小二乘法（OLS）。

2.1 OLS回归模型

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 生成模拟数据
np.random.seed(0)
n = 100
X = np.random.rand(n, 2) # 两个自变量
y = 2*X[:, 0] + 3*X[:, 1] + np.random.randn(n) * 0.1 # 带有噪声的因变量

# 将数据转换为DataFrame
data = pd.DataFrame({'y': y, 'x1': X[:, 0], 'x2': X[:, 1]})

# 使用OLS进行回归分析
model = smf.ols('y ~ x1 + x2', data=data)
results = model.fit()

# 打印回归结果
print(results.summary())

上述代码首先生成了一组模拟数据，其中y是因变量，x1和x2是自变量。然后，使用smf.ols函数创建了一个OLS回归模型，其中'y ~ x1 + x2'是模型公式，表示y是x1和x2的函数。data=data指定了数据源。model.fit()方法用于拟合模型，并将结果存储在results对象中。最后，results.summary()方法打印了回归结果的详细信息，包括回归系数、标准误差、t值、p值、置信区间和R方等。

结果解读:

results.summary()输出的表格包含了以下信息：

Dependent Variable: 因变量的名称。
Model: 使用的模型类型 (OLS)。
Method: 估计方法 (Least Squares)。
Date: 运行日期。
Time: 运行时间。
No. Observations: 观测值的数量。
Df Residuals: 残差的自由度 (观测数量减去参数数量)。
Df Model: 模型自由度 (自变量的数量)。
Covariance Type: 协方差类型 (nonrobust表示未使用任何稳健标准误差)。
R-squared: 决定系数，表示模型解释因变量变异的程度。
Adj. R-squared: 调整后的决定系数，考虑了模型中自变量的数量。
F-statistic: F统计量，用于检验模型整体的显著性。
Prob (F-statistic): F统计量的p值，表示模型整体不显著的可能性。
Log-Likelihood: 对数似然函数值。
AIC: Akaike信息准则，用于比较不同模型的优劣。
BIC: Bayesian信息准则，用于比较不同模型的优劣。
coef: 回归系数的估计值。
std err: 回归系数的标准误差。
t: t统计量，用于检验单个回归系数的显著性。
P>|t|: t统计量的p值，表示回归系数不显著的可能性。
[0.025 0.975]: 回归系数的95%置信区间。
Omnibus: Omnibus检验，用于检验残差是否服从正态分布。
Prob(Omnibus): Omnibus检验的p值。
Skew: 偏度，用于描述残差分布的对称性。
Kurtosis: 峰度，用于描述残差分布的尖峭程度。
Durbin-Watson: Durbin-Watson统计量，用于检验残差是否存在自相关。
Jarque-Bera (JB): Jarque-Bera检验，用于检验残差是否服从正态分布。
Prob(JB): Jarque-Bera检验的p值。
Cond. No.: 条件数，用于评估模型是否存在多重共线性。

2.2 多重共线性

多重共线性是指自变量之间存在高度相关性。多重共线性会导致回归系数的估计值不稳定，并可能影响模型的预测能力。可以使用方差膨胀因子（VIF）来检测多重共线性。

from statsmodels.stats.outliers_influence import variance_inflation_factor

# 计算VIF
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(X, i) for i in range(X.shape[1])]
vif["features"] = ['x1', 'x2']
print(vif)

上述代码使用variance_inflation_factor函数计算了每个自变量的VIF。一般来说，VIF大于5或10表示存在严重的多重共线性。

2.3 异方差性

异方差性是指残差的方差不是常数。异方差性会影响回归系数的标准误差，并可能导致错误的假设检验结果。可以使用Breusch-Pagan检验或White检验来检测异方差性。

import statsmodels.stats.diagnostic as diag

# Breusch-Pagan检验
bp_test = diag.het_breuschpagan(results.resid, results.model.exog)
print("Breusch-Pagan test statistic:", bp_test[0])
print("Breusch-Pagan test p-value:", bp_test[1])

# White检验
white_test = diag.het_white(results.resid, results.model.exog)
print("White test statistic:", white_test[0])
print("White test p-value:", white_test[1])

上述代码分别使用了het_breuschpagan和het_white函数进行了Breusch-Pagan检验和White检验。如果p值小于显著性水平（例如0.05），则拒绝原假设，认为存在异方差性。

如果存在异方差性，可以使用加权最小二乘法（WLS）或稳健标准误差来解决。

2.4 自相关性

自相关性是指残差之间存在相关性。自相关性会影响回归系数的标准误差，并可能导致错误的假设检验结果。可以使用Durbin-Watson检验来检测自相关性。

# Durbin-Watson检验
dw_test = sm.stats.durbin_watson(results.resid)
print("Durbin-Watson statistic:", dw_test)

Durbin-Watson统计量的值在0到4之间。如果统计量接近2，则表示不存在自相关性。如果统计量接近0，则表示存在正自相关性。如果统计量接近4，则表示存在负自相关性。

如果存在自相关性，可以使用广义最小二乘法（GLS）或Newey-West标准误差来解决。

3. 广义线性模型（GLM）

广义线性模型（GLM）是线性回归的推广，允许因变量服从非正态分布。Statsmodels提供了多种GLM模型，包括Logistic回归、泊松回归等。

3.1 Logistic回归

Logistic回归用于预测二元因变量。

# 生成模拟数据
np.random.seed(0)
n = 100
X = np.random.rand(n, 2)
p = 1 / (1 + np.exp(-(2*X[:, 0] + 3*X[:, 1]))) # 计算概率
y = np.random.binomial(1, p, size=n) # 生成二元因变量

# 将数据转换为DataFrame
data = pd.DataFrame({'y': y, 'x1': X[:, 0], 'x2': X[:, 1]})

# 使用GLM进行Logistic回归
model = smf.glm('y ~ x1 + x2', data=data, family=sm.families.Binomial())
results = model.fit()

# 打印回归结果
print(results.summary())

上述代码使用smf.glm函数创建了一个Logistic回归模型，其中family=sm.families.Binomial()指定了因变量服从二项分布。

3.2 泊松回归

泊松回归用于预测计数型因变量。

# 生成模拟数据
np.random.seed(0)
n = 100
X = np.random.rand(n, 2)
rate = np.exp(2*X[:, 0] + 3*X[:, 1]) # 计算发生率
y = np.random.poisson(rate, size=n) # 生成计数型因变量

# 将数据转换为DataFrame
data = pd.DataFrame({'y': y, 'x1': X[:, 0], 'x2': X[:, 1]})

# 使用GLM进行泊松回归
model = smf.glm('y ~ x1 + x2', data=data, family=sm.families.Poisson())
results = model.fit()

# 打印回归结果
print(results.summary())

上述代码使用smf.glm函数创建了一个泊松回归模型，其中family=sm.families.Poisson()指定了因变量服从泊松分布。

4. 假设检验

Statsmodels提供了多种假设检验方法，用于检验关于模型参数的假设。

4.1 t检验

t检验用于检验单个回归系数的显著性。

# t检验
t_test = results.t_test('x1 = 0')
print(t_test)

上述代码使用results.t_test方法对x1的回归系数进行了t检验，检验假设是x1的回归系数等于0。

4.2 F检验

F检验用于检验多个回归系数的联合显著性。

# F检验
f_test = results.f_test('x1 = 0, x2 = 0')
print(f_test)

上述代码使用results.f_test方法对x1和x2的回归系数进行了F检验，检验假设是x1和x2的回归系数都等于0。

4.3 Wald检验

Wald检验是一种通用的假设检验方法，可以用于检验关于模型参数的任何线性约束。

# Wald检验
R = np.array([[1, 0, 0], [0, 1, 0]]) # 约束矩阵
q = np.array([0, 0]) # 约束值
wald_test = results.wald_test((R, q))
print(wald_test)

上述代码使用results.wald_test方法对x1和x2的回归系数进行了Wald检验，检验假设是x1和x2的回归系数都等于0。

5. 模型诊断

模型诊断是评估模型是否满足假设的重要步骤。Statsmodels提供了多种模型诊断工具。

5.1 残差分析

残差分析用于评估残差是否满足正态性、同方差性和独立性等假设。

import matplotlib.pyplot as plt

# 残差图
plt.scatter(results.fittedvalues, results.resid)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()

# QQ图
sm.qqplot(results.resid, line='s')
plt.show()

上述代码绘制了残差图和QQ图，用于评估残差的分布情况。

5.2 Cook’s distance

Cook’s distance用于识别对模型有较大影响的观测值。

# Cook's distance
influence = results.get_influence()
cooks_d = influence.cooks_distance
plt.stem(np.arange(len(cooks_d[0])), cooks_d[0], markerfmt=",")
plt.show()

上述代码计算了每个观测值的Cook’s distance，并绘制了Cook’s distance图。一般来说，Cook’s distance大于4/n的观测值被认为是具有较大影响力的观测值。

6. 时间序列分析

Statsmodels也提供了强大的时间序列分析功能，包括ARIMA模型、季节性分解等。由于时间关系，这里只简单介绍ARIMA模型。

# 生成模拟时间序列数据
np.random.seed(0)
n = 100
time_series = np.cumsum(np.random.randn(n))

# 使用ARIMA模型
from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(time_series, order=(5, 1, 0)) # AR=5, I=1, MA=0
results = model.fit()

# 预测未来值
predictions = results.predict(start=n, end=n+10)
print(predictions)

这段代码展示了如何使用Statsmodels拟合一个ARIMA模型并进行预测。其中order=(5, 1, 0)指定了AR、I和MA的阶数。

7. 代码调试建议

数据预处理： 确保数据类型正确，处理缺失值和异常值。
模型公式： 仔细检查模型公式，确保自变量和因变量之间的关系正确。
参数设置： 了解模型的参数含义，并根据数据特点进行合理设置。
错误信息： 仔细阅读错误信息，并根据提示进行调试。
文档查阅： 查阅Statsmodels官方文档，了解模型和函数的详细用法。
逐步调试： 将代码分解成小块，逐步调试，更容易发现问题。
版本兼容： 确保使用的Statsmodels版本与Python和其他库的版本兼容。

8. 表格总结：常用 Statsmodels 函数和类

函数/类	描述
`smf.ols()`	创建普通最小二乘法（OLS）回归模型。
`smf.glm()`	创建广义线性模型（GLM）。
`results.summary()`	打印模型结果的详细信息。
`variance_inflation_factor()`	计算方差膨胀因子（VIF），用于检测多重共线性。
`diag.het_breuschpagan()`	进行 Breusch-Pagan 检验，用于检测异方差性。
`diag.het_white()`	进行 White 检验，用于检测异方差性。
`sm.stats.durbin_watson()`	进行 Durbin-Watson 检验，用于检测自相关性。
`results.t_test()`	进行 t 检验，用于检验单个回归系数的显著性。
`results.f_test()`	进行 F 检验，用于检验多个回归系数的联合显著性。
`results.wald_test()`	进行 Wald 检验，用于检验关于模型参数的任何线性约束。
`results.resid`	获取模型的残差。
`results.fittedvalues`	获取模型的拟合值。
`results.get_influence()`	获取模型的影响力度量，例如 Cook’s distance。
`ARIMA()`	创建 ARIMA 时间序列模型。

关键信息回顾

今天我们学习了如何使用Statsmodels进行回归分析和假设检验。Statsmodels提供了强大的统计建模工具，能够帮助我们理解数据、建立模型和进行推断。我们学习了线性回归、广义线性模型、模型诊断和假设检验等内容。掌握这些知识能够帮助我们更好地进行数据分析和决策。

Python统计建模：利用Statsmodels进行回归分析和假设检验

发表回复 取消回复

Python统计建模：利用`Statsmodels`进行回归分析和假设检验

发表回复取消回复