Statsmodels 统计建模：构建复杂统计模型与假设检验

好的，没问题！咱们今天就来聊聊 Statsmodels 这个神器，一起玩转复杂的统计模型和假设检验。别担心，我会尽量用大白话，让你听得懂，学得会，还能笑出声。

开场白：Statsmodels，统计界的瑞士军刀

各位观众，晚上好！（或者早上好，取决于你啥时候看这篇文章）。今天咱们的主题是 Statsmodels。这玩意儿啊，就像统计界的瑞士军刀，啥都能干一点。你想搞回归分析？没问题！你想做时间序列？小菜一碟！你想做假设检验？安排！

总之，Statsmodels 是 Python 中一个强大的统计建模库，它提供了各种统计模型、假设检验、以及数据探索的功能。如果你想用 Python 做一些严肃的统计分析，那 Statsmodels 绝对是你的好帮手。

第一部分：Statsmodels 的基础操作

首先，咱们得先把 Statsmodels 安装好。如果你用的是 Anaconda，那就直接：

conda install statsmodels

如果你用的是 pip，那就：

pip install statsmodels

安装完了，咱们就可以开始玩了。

1. 导入 Statsmodels

import statsmodels.api as sm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

这里，sm 是 Statsmodels 的常用别名，numpy 用来处理数值数据，pandas 用来处理表格数据，matplotlib 用来画图。

2. 数据准备

咱们先来搞点数据。这里咱们用一个简单的例子，假设我们想研究广告投入和销售额之间的关系。

# 模拟数据
np.random.seed(0) # 设置随机种子，保证结果可重复
X = np.random.rand(100, 1) * 10  # 广告投入 (0-10)
y = 2 * X + np.random.randn(100, 1) * 2  # 销售额 (带噪声)

# 将数据转换为 Pandas DataFrame
data = pd.DataFrame({'广告投入': X.flatten(), '销售额': y.flatten()})
print(data.head())

这段代码模拟了一些广告投入和销售额的数据，并且加入了随机噪声，让数据更真实一些。flatten()函数用于将二维数组转换为一维数组，方便存储到DataFrame中。

3. 线性回归模型

咱们用 Statsmodels 来做一个简单的线性回归。

# 添加截距项 (constant)
X = sm.add_constant(X)

# 创建 OLS (Ordinary Least Squares) 模型
model = sm.OLS(y, X)

# 拟合模型
results = model.fit()

# 打印模型结果
print(results.summary())

这段代码做了几件事：

• sm.add_constant(X)：给自变量 X 添加一个常数项（截距），这是线性回归中必不可少的一步。
• sm.OLS(y, X)：创建一个 OLS 模型，其中 y 是因变量（销售额），X 是自变量（广告投入）。
• model.fit()：用最小二乘法拟合模型。
• results.summary()：打印模型的结果，包括回归系数、标准误差、t 值、p 值、R 方等等。

结果会告诉你广告投入对销售额的影响有多大，以及这个影响是否显著。results.summary() 的输出非常详细，包含了模型的所有重要信息。

第二部分：高级模型：广义线性模型 (GLM)

线性回归很棒，但它也有局限性。比如，当因变量不是连续的，而是二元的（比如“是”或“否”）或者计数型的（比如“有多少个”），线性回归就不太适用了。这时候，咱们就需要广义线性模型 (GLM) 出马了。

GLM 允许咱们使用不同的分布来拟合数据，比如 Logistic 分布（用于二元数据）或者 Poisson 分布（用于计数数据）。

1. Logistic 回归

假设咱们想预测用户是否会点击广告，点击（1）或者不点击（0）。

# 模拟数据
np.random.seed(0)
X = np.random.rand(100, 1) * 10  # 用户特征 (0-10)
p = 1 / (1 + np.exp(-(0.5 * X - 3)))  # 点击概率
y = np.random.binomial(1, p.flatten(), 100)  # 模拟点击行为 (0 或 1)

# 将数据转换为 Pandas DataFrame
data = pd.DataFrame({'用户特征': X.flatten(), '是否点击': y})
print(data.head())

# 添加截距项
X = sm.add_constant(X)

# 创建 Logistic 模型
model = sm.GLM(y, X, family=sm.families.Binomial())

# 拟合模型
results = model.fit()

# 打印模型结果
print(results.summary())

这段代码的关键在于 sm.GLM(y, X, family=sm.families.Binomial())。这里，family=sm.families.Binomial() 指定了使用 Logistic 分布，也就是二项分布。

2. Poisson 回归

假设咱们想预测每天有多少顾客光顾店铺。

# 模拟数据
np.random.seed(0)
X = np.random.rand(100, 1) * 10  # 营销投入 (0-10)
rate = np.exp(0.2 * X)  # 顾客光顾率
y = np.random.poisson(rate.flatten(), 100)  # 模拟顾客数量

# 将数据转换为 Pandas DataFrame
data = pd.DataFrame({'营销投入': X.flatten(), '顾客数量': y})
print(data.head())

# 添加截距项
X = sm.add_constant(X)

# 创建 Poisson 模型
model = sm.GLM(y, X, family=sm.families.Poisson())

# 拟合模型
results = model.fit()

# 打印模型结果
print(results.summary())

这段代码的关键在于 sm.GLM(y, X, family=sm.families.Poisson())。这里，family=sm.families.Poisson() 指定了使用 Poisson 分布。

第三部分：假设检验：检验你的猜想

统计分析不仅仅是拟合模型，更重要的是检验你的猜想。假设检验就是用来干这个的。

1. t 检验

t 检验用于检验样本均值是否与某个值相等。比如，咱们想检验某个班级的平均成绩是否高于 80 分。

from scipy import stats

# 模拟数据
np.random.seed(0)
scores = np.random.normal(75, 10, 100)  # 模拟成绩 (平均 75, 标准差 10)

# 进行单样本 t 检验
t_statistic, p_value = stats.ttest_1samp(scores, 80)

print("t 统计量:", t_statistic)
print("p 值:", p_value)

# 判断是否拒绝原假设 (原假设: 平均成绩等于 80)
alpha = 0.05  # 显著性水平
if p_value < alpha:
    print("拒绝原假设，平均成绩显著不等于 80")
else:
    print("不能拒绝原假设，平均成绩可能等于 80")

这段代码做了以下几件事：

• stats.ttest_1samp(scores, 80)：进行单样本 t 检验，其中 scores 是样本数据，80 是要检验的均值。
• t_statistic 是 t 统计量，p_value 是 p 值。
• alpha = 0.05：设置显著性水平为 0.05，也就是容忍 5% 的犯错概率。
• 如果 p_value < alpha，则拒绝原假设，认为样本均值显著不等于 80。否则，不能拒绝原假设。

2. 方差分析 (ANOVA)

方差分析用于比较多个组的均值是否有显著差异。比如，咱们想比较三种不同的广告策略对销售额的影响。

# 模拟数据
np.random.seed(0)
strategy_A = np.random.normal(100, 15, 50)  # 策略 A 的销售额
strategy_B = np.random.normal(110, 15, 50)  # 策略 B 的销售额
strategy_C = np.random.normal(120, 15, 50)  # 策略 C 的销售额

# 进行单因素方差分析
f_statistic, p_value = stats.f_oneway(strategy_A, strategy_B, strategy_C)

print("F 统计量:", f_statistic)
print("p 值:", p_value)

# 判断是否拒绝原假设 (原假设: 三种策略的平均销售额相等)
alpha = 0.05  # 显著性水平
if p_value < alpha:
    print("拒绝原假设，至少有一种策略的平均销售额与其他策略不同")
else:
    print("不能拒绝原假设，三种策略的平均销售额可能相等")

这段代码做了以下几件事：

• stats.f_oneway(strategy_A, strategy_B, strategy_C)：进行单因素方差分析，其中 strategy_A, strategy_B, strategy_C 是三组样本数据。
• f_statistic 是 F 统计量，p_value 是 p 值。
• 如果 p_value < alpha，则拒绝原假设，认为至少有一种策略的平均销售额与其他策略不同。否则，不能拒绝原假设。

3. 卡方检验

卡方检验用于检验分类变量之间是否存在关联。比如，咱们想检验性别和是否喜欢某个产品之间是否存在关联。

from scipy.stats import chi2_contingency

# 创建列联表
observed = np.array([[30, 20], [20, 30]]) # 男生喜欢/不喜欢，女生喜欢/不喜欢

# 进行卡方检验
chi2, p, dof, expected = chi2_contingency(observed)

print("卡方统计量:", chi2)
print("p 值:", p)
print("自由度:", dof)
print("期望值:", expected)

# 判断是否拒绝原假设 (原假设: 性别和是否喜欢该产品没有关联)
alpha = 0.05  # 显著性水平
if p < alpha:
    print("拒绝原假设，性别和是否喜欢该产品之间存在关联")
else:
    print("不能拒绝原假设，性别和是否喜欢该产品之间可能没有关联")

这段代码做了以下几件事：

• observed = np.array([[30, 20], [20, 30]])：创建一个列联表，表示观察到的数据。
• chi2_contingency(observed)：进行卡方检验。
• chi2 是卡方统计量，p 是 p 值，dof 是自由度，expected 是期望值。
• 如果 p < alpha，则拒绝原假设，认为性别和是否喜欢该产品之间存在关联。否则，不能拒绝原假设。

第四部分：时间序列分析：预测未来

Statsmodels 也是时间序列分析的一把好手。时间序列分析用于分析随时间变化的数据，比如股票价格、气温、销售额等等。

1. ARIMA 模型

ARIMA (Autoregressive Integrated Moving Average) 模型是时间序列分析中最常用的模型之一。它结合了自回归 (AR)、差分 (I) 和移动平均 (MA) 三种成分。

import statsmodels.tsa.arima.model as arima
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 模拟时间序列数据
np.random.seed(0)
n = 200
ar_params = np.array([0.7])
ma_params = np.array([0.3])
errors = np.random.normal(0, 1, n)
data = arima.arma_generate_sample(ar_params, ma_params, nsample=n, scale=1, distrvs=np.random.normal)

# 绘制时间序列图
plt.figure(figsize=(12, 4))
plt.plot(data)
plt.title("模拟时间序列数据")
plt.show()

# 进行 ADF 单位根检验 (检验时间序列是否平稳)
result = adfuller(data)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('t%s: %.3f' % (key, value))

# 绘制自相关函数 (ACF) 和偏自相关函数 (PACF) 图
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(data, ax=axes[0], lags=20)
plot_pacf(data, ax=axes[1], lags=20)
plt.show()

# 拟合 ARIMA 模型
# 根据 ACF 和 PACF 图选择合适的 p, d, q 值
model = arima.ARIMA(data, order=(1, 0, 1)) # 这里p=1, d=0, q=1
results = model.fit()

# 打印模型结果
print(results.summary())

# 进行预测
predictions = results.predict(start=n, end=n+10) # 预测未来10个时间点的值
print("预测结果:", predictions)

这段代码比较复杂，咱们一步一步来：

• 模拟时间序列数据：arima.arma_generate_sample 函数用于生成 ARMA 模型的数据。
• 绘制时间序列图：用于观察数据的趋势和季节性。
• ADF 单位根检验：用于检验时间序列是否平稳。平稳的时间序列更容易进行预测。如果 p 值小于显著性水平，则拒绝原假设，认为时间序列是平稳的。
• 绘制自相关函数 (ACF) 和偏自相关函数 (PACF) 图：用于选择合适的 ARIMA 模型的 p, d, q 值。
• 拟合 ARIMA 模型：arima.ARIMA(data, order=(p, d, q)) 创建一个 ARIMA 模型，其中 p 是自回归阶数，d 是差分阶数，q 是移动平均阶数。
• 进行预测：results.predict(start=n, end=n+10) 预测未来 10 个时间点的值。

第五部分：模型评估与选择

模型建好了，咱们还得评估一下模型的好坏，以及选择哪个模型更好。

1. R 方 (R-squared)

R 方是衡量模型拟合程度的指标，取值范围是 0 到 1。R 方越大，模型拟合得越好。

# 以线性回归为例
# ... (前面的代码) ...

# 打印 R 方
print("R 方:", results.rsquared)

2. AIC 和 BIC

AIC (Akaike Information Criterion) 和 BIC (Bayesian Information Criterion) 是用于模型选择的指标。AIC 和 BIC 越小，模型越好。

# 以线性回归为例
# ... (前面的代码) ...

# 打印 AIC 和 BIC
print("AIC:", results.aic)
print("BIC:", results.bic)

3. 交叉验证

交叉验证是一种评估模型泛化能力的常用方法。它将数据分成训练集和测试集，用训练集训练模型，然后用测试集评估模型。

from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

# 以线性回归为例
# ... (前面的代码) ...

# 创建 K 折交叉验证
kf = KFold(n_splits=5, shuffle=True, random_state=0) # 5折交叉验证

# 存储每次迭代的 MSE
mse_scores = []

# 循环进行交叉验证
for train_index, test_index in kf.split(X):
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]

    # 创建和拟合模型
    model = sm.OLS(y_train, X_train)
    results = model.fit()

    # 预测测试集
    y_pred = results.predict(X_test)

    # 计算 MSE
    mse = mean_squared_error(y_test, y_pred)
    mse_scores.append(mse)

# 打印平均 MSE
print("平均 MSE:", np.mean(mse_scores))

这段代码做了以下几件事：

• KFold(n_splits=5, shuffle=True, random_state=0)：创建一个 5 折交叉验证对象。
• 循环进行交叉验证：每次迭代都将数据分成训练集和测试集，用训练集训练模型，然后用测试集评估模型。
• mean_squared_error(y_test, y_pred)：计算均方误差 (MSE)，MSE 越小，模型越好。
• 打印平均 MSE：用于评估模型的整体泛化能力。

总结：Statsmodels，你的统计分析好帮手

好了，各位，今天的 Statsmodels 之旅就到这里了。咱们从基础操作开始，学习了线性回归、广义线性模型、假设检验、时间序列分析，以及模型评估与选择。希望你能通过这篇文章，对 Statsmodels 有一个更深入的了解，并且能够用它来解决实际问题。

记住，Statsmodels 只是一个工具，更重要的是你的统计思维。多思考，多实践，你一定能成为统计分析的高手！