入门¶
这个非常简单的案例研究旨在帮助您快速上手使用 statsmodels
。从原始数据开始,我们将展示估计统计模型和绘制诊断图所需的步骤。我们只使用 statsmodels
或其 pandas
和 patsy
依赖项提供的函数。
加载模块和函数¶
在安装 statsmodels 及其依赖项之后,我们加载一些模块和函数
In [1]: import statsmodels.api as sm
In [2]: import pandas
In [3]: from patsy import dmatrices
pandas 基于 numpy
数组提供丰富的 数据结构和数据分析工具。 pandas.DataFrame
函数提供标记的 (可能异构的) 数据数组,类似于 R
的“data.frame”。pandas.read_csv
函数可用于将逗号分隔值文件转换为 DataFrame
对象。
patsy 是一个 Python 库,用于描述统计模型和使用类似于 R
的公式构建设计矩阵。
注意
本示例使用 API 接口。有关导入 API 接口 (statsmodels.api
和 statsmodels.tsa.api
) 与直接从定义模型的模块导入之间的区别,请参见导入路径和结构。
数据¶
我们下载了Guerry 数据集,这是一个用于支持 Andre-Michel Guerry 1833 年的《法国道德统计学论文》的历史数据集合。该数据集以逗号分隔值格式 (CSV) 形式托管在Rdatasets 存储库中。我们可以将文件下载到本地,然后使用 read_csv
加载它,但是 pandas
会自动为我们完成所有这些操作
In [4]: df = sm.datasets.get_rdataset("Guerry", "HistData").data
输入/输出文档页面 显示了如何从其他各种格式导入。
我们选择感兴趣的变量并查看最后 5 行
In [5]: vars = ['Department', 'Lottery', 'Literacy', 'Wealth', 'Region']
In [6]: df = df[vars]
In [7]: df[-5:]
Out[7]:
Department Lottery Literacy Wealth Region
81 Vienne 40 25 68 W
82 Haute-Vienne 55 13 67 C
83 Vosges 14 62 82 E
84 Yonne 51 47 30 C
85 Corse 83 49 37 NaN
请注意,Region 列中有一个缺失值。我们使用 pandas
提供的 DataFrame
方法将其消除
In [8]: df = df.dropna()
In [9]: df[-5:]
Out[9]:
Department Lottery Literacy Wealth Region
80 Vendee 68 28 56 W
81 Vienne 40 25 68 W
82 Haute-Vienne 55 13 67 C
83 Vosges 14 62 82 E
84 Yonne 51 47 30 C
实质动机和模型¶
我们想知道 1820 年代法国 86 个省份的识字率是否与人均皇家彩票赌注有关。我们需要控制每个省份的财富水平,并且我们还想在回归方程的右侧包含一系列虚拟变量,以控制由于区域效应导致的未观察到的异质性。该模型使用普通最小二乘回归 (OLS) 进行估计。
设计矩阵 (endog & exog)¶
要拟合 statsmodels
涵盖的大多数模型,您将需要创建两个设计矩阵。第一个是内生变量矩阵 (即因变量、响应变量、回归变量等)。第二个是外生变量矩阵 (即自变量、预测变量、回归变量等)。OLS 系数估计值按预期计算
其中 \(y\) 是 \(N \times 1\) 人均彩票赌注数据列 (Lottery)。\(X\) 是 \(N \times 7\),包含截距、Literacy 和 Wealth 变量以及 4 个区域二元变量。
patsy
模块提供了一个方便的函数,可以使用类似于 R
的公式来准备设计矩阵。您可以在此处找到更多信息。
我们使用 patsy
的 dmatrices
函数来创建设计矩阵
In [10]: y, X = dmatrices('Lottery ~ Literacy + Wealth + Region', data=df, return_type='dataframe')
生成的矩阵/数据框如下所示
In [11]: y[:3]
Out[11]:
Lottery
0 41.0
1 38.0
2 66.0
In [12]: X[:3]
Out[12]:
Intercept Region[T.E] Region[T.N] ... Region[T.W] Literacy Wealth
0 1.0 1.0 0.0 ... 0.0 37.0 73.0
1 1.0 0.0 1.0 ... 0.0 51.0 22.0
2 1.0 0.0 0.0 ... 0.0 13.0 61.0
[3 rows x 7 columns]
请注意,dmatrices
已
将分类变量 Region 分割成一组指示变量。
在外部回归变量矩阵中添加了一个常数。
返回了
pandas
DataFrame 而不是简单的 numpy 数组。这很有用,因为 DataFrame 允许statsmodels
在报告结果时保留元数据 (例如变量名称)。
当然,上述行为可以更改。请参见patsy 文档页面。
模型拟合和摘要¶
在 statsmodels
中拟合模型通常涉及 3 个简单的步骤
使用模型类描述模型
使用类方法拟合模型
使用摘要方法检查结果
对于 OLS,这是通过以下方式实现的
In [13]: mod = sm.OLS(y, X) # Describe model
In [14]: res = mod.fit() # Fit model
In [15]: print(res.summary()) # Summarize model
OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.338
Model: OLS Adj. R-squared: 0.287
Method: Least Squares F-statistic: 6.636
Date: Thu, 03 Oct 2024 Prob (F-statistic): 1.07e-05
Time: 16:15:27 Log-Likelihood: -375.30
No. Observations: 85 AIC: 764.6
Df Residuals: 78 BIC: 781.7
Df Model: 6
Covariance Type: nonrobust
===============================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 38.6517 9.456 4.087 0.000 19.826 57.478
Region[T.E] -15.4278 9.727 -1.586 0.117 -34.793 3.938
Region[T.N] -10.0170 9.260 -1.082 0.283 -28.453 8.419
Region[T.S] -4.5483 7.279 -0.625 0.534 -19.039 9.943
Region[T.W] -10.0913 7.196 -1.402 0.165 -24.418 4.235
Literacy -0.1858 0.210 -0.886 0.378 -0.603 0.232
Wealth 0.4515 0.103 4.390 0.000 0.247 0.656
==============================================================================
Omnibus: 3.049 Durbin-Watson: 1.785
Prob(Omnibus): 0.218 Jarque-Bera (JB): 2.694
Skew: -0.340 Prob(JB): 0.260
Kurtosis: 2.454 Cond. No. 371.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
res
对象具有许多有用的属性。例如,我们可以通过键入以下内容来提取参数估计值和 R 平方
In [16]: res.params
Out[16]:
Intercept 38.651655
Region[T.E] -15.427785
Region[T.N] -10.016961
Region[T.S] -4.548257
Region[T.W] -10.091276
Literacy -0.185819
Wealth 0.451475
dtype: float64
In [17]: res.rsquared
Out[17]: np.float64(0.337950869192882)
键入 dir(res)
以获取属性的完整列表。
有关更多信息和示例,请参见回归文档页面
诊断和规格测试¶
statsmodels
允许您进行一系列有用的回归诊断和规格测试。例如,应用彩虹检验以检查线性 (原假设是关系被正确地建模为线性)
In [18]: sm.stats.linear_rainbow(res)
Out[18]: (np.float64(0.847233997615691), np.float64(0.6997965543621644))
诚然,上面生成的输出不是特别详细,但我们从阅读docstring (还有 print(sm.stats.linear_rainbow.__doc__)
) 中知道,第一个数字是 F 统计量,第二个数字是 p 值。
statsmodels
还提供图形函数。例如,我们可以通过以下方式为一组回归变量绘制偏回归图
In [19]: sm.graphics.plot_partregress('Lottery', 'Wealth', ['Region', 'Literacy'],
....: data=df, obs_labels=False)
....:
Out[19]: <Figure size 640x480 with 1 Axes>
文档¶
可以使用 webdoc
从 IPython 会话访问文档。
打开浏览器并显示在线文档 |
更多¶
恭喜!您已准备好继续学习目录中的其他主题