9.4 Pandas 与 Statsmodels

9.4 Pandas 与 Statsmodels

9.4 Pandas 与 Statsmodels 的集成

9.4.1 数据准备

首先，我们需要准备数据，并将其转换为 Pandas DataFrame。Statsmodels 可以直接使用 DataFrame 作为输入，这简化了数据处理步骤。

import pandas as pd
import numpy as np
# 创建示例数据
np.random.seed(0)
data = {
    'X1': np.random.rand(100),
    'X2': np.random.rand(100),
    'Y': 2 * np.random.rand(100) + 3 * np.random.rand(100) + np.random.randn(100)
}
df = pd.DataFrame(data)
print(df.head())
​

这段代码创建了一个包含三个变量（X1, X2, Y）的 DataFrame。X1 和 X2 是自变量，Y 是因变量。

9.4.2 使用 Statsmodels 进行线性回归

Statsmodels 提供了多种统计模型，其中线性回归是最常用的之一。我们可以使用 statsmodels.formula.api 模块中的 ols 函数来拟合线性回归模型。

import statsmodels.formula.api as smf
# 使用公式接口拟合线性回归模型
model = smf.ols('Y ~ X1 + X2', data=df)
results = model.fit()
# 打印模型摘要
print(results.summary())
​

代码解释:

1. import statsmodels.formula.api as smf: 导入 Statsmodels 的公式接口模块。这个模块允许我们使用类似于 R 的公式语法来指定模型。

2. model = smf.ols('Y ~ X1 + X2', data=df): 使用 ols 函数创建线性回归模型。

   • 'Y ~ X1 + X2' 是公式，表示 Y 是因变量，X1 和 X2 是自变量。~ 符号表示“由…解释”。

   • data=df 指定使用 df DataFrame 作为数据源。

3. results = model.fit(): 使用最小二乘法拟合模型。fit() 方法返回一个包含模型结果的对象。

4. print(results.summary()): 打印模型摘要。摘要包含了模型的各种统计信息，例如 R-squared、系数、标准误差、t 值和 p 值。

模型摘要解读:

模型摘要提供了关于模型拟合效果的重要信息。以下是一些关键指标：

• R-squared: 决定系数，表示模型解释因变量变异的比例。值越高，模型拟合得越好。

• Adj. R-squared: 调整后的决定系数，考虑了模型中自变量的数量。

• coef: 系数，表示自变量对因变量的影响程度。

• std err: 标准误差，表示系数估计的准确性。

• t: t 值，用于检验系数是否显著不为零。

• P>|t|: p 值，表示在零假设（系数为零）下观察到当前或更极端结果的概率。通常，p 值小于 0.05 表示系数是显著的。

• Confidence Interval: 置信区间，表示系数的可能取值范围。

9.4.3 使用 Statsmodels 进行广义线性模型 (GLM)

广义线性模型 (GLM) 是一种更通用的模型，可以处理非正态分布的因变量。例如，我们可以使用 GLM 来拟合泊松回归模型，用于计数数据。

import statsmodels.api as sm
import numpy as np
# 创建泊松分布的因变量
np.random.seed(0)
df['Y_poisson'] = np.random.poisson(lam=5, size=100)
# 拟合泊松回归模型
model_poisson = sm.GLM(df['Y_poisson'], df[['X1', 'X2']], family=sm.families.Poisson())
results_poisson = model_poisson.fit()
# 打印模型摘要
print(results_poisson.summary())
​

代码解释:

1. df['Y_poisson'] = np.random.poisson(lam=5, size=100): 创建一个服从泊松分布的因变量 Y_poisson。lam=5 表示泊松分布的均值为 5。

2. model_poisson = sm.GLM(df['Y_poisson'], df[['X1', 'X2']], family=sm.families.Poisson()): 使用 GLM 函数创建泊松回归模型。

   • df['Y_poisson'] 是因变量。

   • df[['X1', 'X2']] 是自变量。

   • family=sm.families.Poisson() 指定使用泊松分布作为误差分布。

3. results_poisson = model_poisson.fit(): 拟合模型。

4. print(results_poisson.summary()): 打印模型摘要。

9.4.4 使用 DataFrame 的列名直接进行建模

Statsmodels 可以直接使用 DataFrame 的列名进行建模，无需手动创建设计矩阵。这简化了建模过程。

import statsmodels.formula.api as smf
# 使用列名直接建模
model = smf.ols('Y ~ X1 + X2', data=df)
results = model.fit()
# 打印模型摘要
print(results.summary())
​

这与之前的例子相同，但强调了使用 DataFrame 列名的便利性。

9.4.5 进行预测

拟合模型后，我们可以使用 predict() 方法进行预测。

# 创建新的数据用于预测
new_data = pd.DataFrame({'X1': [0.5, 0.75], 'X2': [0.2, 0.9]})
# 进行预测
predictions = results.predict(new_data)
print(predictions)
​

代码解释:

1. new_data = pd.DataFrame({'X1': [0.5, 0.75], 'X2': [0.2, 0.9]}): 创建一个新的 DataFrame，包含用于预测的自变量值。

2. predictions = results.predict(new_data): 使用 predict() 方法进行预测。该方法接受一个新的 DataFrame 作为输入，并返回预测值。

3. print(predictions): 打印预测值。

9.4.6 假设检验

Statsmodels 提供了各种假设检验工具，可以用于检验模型的假设。例如，我们可以使用 statsmodels.stats.diagnostic.het_breuschpagan 函数检验异方差性。

import statsmodels.stats.diagnostic as diag
# 检验异方差性
bp_test = diag.het_breuschpagan(results.resid, results.model.exog)
print("Breusch-Pagan test:")
print("LM statistic:", bp_test[0])
print("LM p-value:", bp_test[1])
print("F statistic:", bp_test[2])
print("F p-value:", bp_test[3])
​

代码解释:

1. bp_test = diag.het_breuschpagan(results.resid, results.model.exog): 使用 het_breuschpagan 函数进行 Breusch-Pagan 检验。

   • results.resid 是模型的残差。

   • results.model.exog 是模型的自变量。

2. print("Breusch-Pagan test:"): 打印检验结果的标题。

3. print("LM statistic:", bp_test[0]): 打印 LM 统计量。

4. print("LM p-value:", bp_test[1]): 打印 LM p 值。如果 p 值小于显著性水平（例如 0.05），则拒绝原假设，认为存在异方差性。

5. print("F statistic:", bp_test[2]): 打印 F 统计量。

6. print("F p-value:", bp_test[3]): 打印 F p 值。

9.4.7 时间序列分析

Statsmodels 也提供了强大的时间序列分析工具。我们可以使用 statsmodels.tsa.arima.model.ARIMA 类来拟合 ARIMA 模型。

import statsmodels.tsa.arima.model as arima
import pandas as pd
import numpy as np
# 创建时间序列数据
np.random.seed(0)
dates = pd.date_range('2023-01-01', periods=100, freq='D')
data = np.random.randn(100)
ts = pd.Series(data, index=dates)
# 拟合 ARIMA 模型
model_arima = arima.ARIMA(ts, order=(5, 1, 0)) # 示例参数：AR(5), I(1), MA(0)
results_arima = model_arima.fit()
# 打印模型摘要
print(results_arima.summary())
​

代码解释:

1. dates = pd.date_range('2023-01-01', periods=100, freq='D'): 创建一个日期索引，从 2023-01-01 开始，持续 100 天，频率为每天。

2. data = np.random.randn(100): 创建一个包含 100 个随机数的时间序列数据。

3. ts = pd.Series(data, index=dates): 创建一个 Pandas Series，将数据与日期索引关联起来。

4. model_arima = arima.ARIMA(ts, order=(5, 1, 0)): 使用 ARIMA 类创建 ARIMA 模型。

   • ts 是时间序列数据。

   • order=(5, 1, 0) 指定 ARIMA 模型的阶数。p=5 表示自回归 (AR) 阶数为 5，d=1 表示差分 (I) 阶数为 1，q=0 表示移动平均 (MA) 阶数为 0。

5. results_arima = model_arima.fit(): 拟合模型。

6. print(results_arima.summary()): 打印模型摘要。

9.4.8 使用 Mermaid 绘制流程图

以下是一个使用 Mermaid 绘制的流程图，展示了 Pandas 和 Statsmodels 集成的典型流程：

流程图解释:

1. A[Pandas DataFrame]: 流程从 Pandas DataFrame 开始，这是数据的存储和操作中心。

2. B(数据准备): 对数据进行清洗、转换和预处理，以满足模型的要求。

3. C{选择模型}: 根据数据的特点和研究目标选择合适的统计模型。

4. D[smf.ols]: 如果选择线性回归模型，则使用 statsmodels.formula.api.ols 函数创建模型。

5. E[sm.GLM]: 如果选择广义线性模型，则使用 statsmodels.api.GLM 函数创建模型。

6. F(model.fit()): 使用 fit() 方法拟合模型。

7. G{模型诊断}: 对模型进行诊断，例如检验残差的正态性、异方差性等。

8. H[results.summary()]: 使用 summary() 方法查看模型摘要，包含各种统计信息。

9. I[假设检验]: 进行假设检验，例如检验系数是否显著不为零。

10. J(模型评估): 根据模型诊断和假设检验的结果评估模型的性能。

11. K{预测}: 使用模型进行预测。

12. L[results.predict()]: 使用 predict() 方法进行预测。

13. M(结果分析): 分析预测结果，并得出结论。

9.4.9 总结

Pandas 和 Statsmodels 的集成提供了一个强大的统计建模平台。Pandas 提供了灵活的数据结构和操作工具，Statsmodels 提供了各种统计模型和假设检验工具。通过将两者结合使用，我们可以简化统计建模流程，提高效率，并进行更深入的数据分析。本节介绍了如何使用 Pandas 和 Statsmodels 进行线性回归、广义线性模型和时间序列分析，以及如何进行预测和假设检验。希望这些代码实践和详细解释能够帮助你更好地理解和使用 Pandas 和 Statsmodels。