【翻译】《利用Python进行数据分析·第2版》第13章（中） Python建模库介绍

作者：SeanCheney   Python爱好者社区专栏作者

简书专栏：https://www.jianshu.com/u/130f76596b02

前文传送门：

13.2 用Patsy创建模型描述

Patsy是Python的一个库，使用简短的字符串“公式语法”描述统计模型（尤其是线性模型），可能是受到了R和S统计编程语言的公式语法的启发。

Patsy适合描述statsmodels的线性模型，因此我会关注于它的主要特点，让你尽快掌握。Patsy的公式是一个特殊的字符串语法，如下所示：

y ~ x0 + x1

a+b不是将a与b相加的意思，而是为模型创建的设计矩阵。patsy.dmatrices函数接收一个公式字符串和一个数据集（可以是DataFrame或数组的字典），为线性模型创建设计矩阵：

In [29]: data = pd.DataFrame({   ....:     'x0': [1, 2, 3, 4, 5],   ....:     'x1': [0.01, -0.01, 0.25, -4.1, 0.],   ....:     'y': [-1.5, 0., 3.6, 1.3, -2.]}) In [30]: data Out[30]:   x0    x1    y 0   1  0.01 -1.5 1   2 -0.01  0.0 2   3  0.25  3.6 3   4 -4.10  1.3 4   5  0.00 -2.0 In [31]: import patsy In [32]: y, X = patsy.dmatrices('y ~ x0 + x1', data)

现在有：

In [33]: y Out[33]: DesignMatrix with shape (5, 1)     y  -1.5   0.0   3.6   1.3  -2.0  Terms:    'y' (column 0) In [34]: X Out[34]: DesignMatrix with shape (5, 3)  Intercept  x0     x1          1   1   0.01          1   2  -0.01          1   3   0.25          1   4  -4.10          1   5   0.00  Terms:    'Intercept' (column 0)    'x0' (column 1)    'x1' (column 2)

这些Patsy的DesignMatrix实例是NumPy的ndarray，带有附加元数据：

In [35]: np.asarray(y) Out[35]: array([[-1.5],       [ 0. ],       [ 3.6],       [ 1.3],       [-2. ]]) In [36]: np.asarray(X) Out[36]: array([[ 1.  ,  1.  ,  0.01],       [ 1.  ,  2.  , -0.01],       [ 1.  ,  3.  ,  0.25],       [ 1.  ,  4.  , -4.1 ],       [ 1.  ,  5.  ,  0.  ]])

你可能想Intercept是哪里来的。这是线性模型（比如普通最小二乘胡桂）的惯例用法。添加 +0 到模型可以不显示intercept：

In [37]: patsy.dmatrices('y ~ x0 + x1 + 0', data)[1] Out[37]: DesignMatrix with shape (5, 2)  x0     x1   1   0.01   2  -0.01   3   0.25   4  -4.10   5   0.00  Terms:    'x0' (column 0)    'x1' (column 1)

Patsy对象可以直接传递到算法（比如numpy.linalg.lstsq）中，它执行普通最小二乘回归：

In [38]: coef, resid, _, _ = np.linalg.lstsq(X, y)

模型的元数据保留在design_info属性中，因此你可以重新附加列名到拟合系数，以获得一个Series，例如：

In [39]: coef Out[39]: array([[ 0.3129],       [-0.0791],       [-0.2655]]) In [40]: coef = pd.Series(coef.squeeze(), index=X.design_info.column_names) In [41]: coef Out[41]: Intercept    0.312910 x0          -0.079106 x1          -0.265464 dtype: float64

用Patsy公式进行数据转换

你可以将Python代码与patsy公式结合。在评估公式时，库将尝试查找在封闭作用域内使用的函数：

In [42]: y, X = patsy.dmatrices('y ~ x0 + np.log(np.abs(x1) + 1)', data) In [43]: X Out[43]: DesignMatrix with shape (5, 3)  Intercept  x0  np.log(np.abs(x1) + 1)          1   1                 0.00995          1   2                 0.00995          1   3                 0.22314          1   4                 1.62924          1   5                 0.00000  Terms:    'Intercept' (column 0)    'x0' (column 1)    'np.log(np.abs(x1) + 1)' (column 2)

常见的变量转换包括标准化（平均值为0，方差为1）和中心化（减去平均值）。Patsy有内置的函数进行这样的工作：

In [44]: y, X = patsy.dmatrices('y ~ standardize(x0) + center(x1)', data) In [45]: X Out[45]: DesignMatrix with shape (5, 3)  Intercept  standardize(x0)  center(x1)          1         -1.41421        0.78          1         -0.70711        0.76          1          0.00000        1.02          1          0.70711       -3.33          1          1.41421        0.77  Terms:    'Intercept' (column 0)    'standardize(x0)' (column 1)    'center(x1)' (column 2)

作为建模的一步，你可能拟合模型到一个数据集，然后用另一个数据集评估模型。另一个数据集可能是剩余的部分或是新数据。当执行中心化和标准化转变，用新数据进行预测要格外小心。因为你必须使用平均值或标准差转换新数据集，这也称作状态转换。

patsy.build_design_matrices函数可以应用于转换新数据，使用原始样本数据集的保存信息：

In [46]: new_data = pd.DataFrame({   ....:     'x0': [6, 7, 8, 9],   ....:     'x1': [3.1, -0.5, 0, 2.3],   ....:     'y': [1, 2, 3, 4]}) In [47]: new_X = patsy.build_design_matrices([X.design_info], new_data) In [48]: new_X Out[48]: [DesignMatrix with shape (4, 3)   Intercept  standardize(x0)  center(x1)           1          2.12132        3.87           1          2.82843        0.27           1          3.53553        0.77           1          4.24264        3.07   Terms:     'Intercept' (column 0)     'standardize(x0)' (column 1)     'center(x1)' (column 2)]

因为Patsy中的加号不是加法的意义，当你按照名称将数据集的列相加时，你必须用特殊I函数将它们封装起来：

In [49]: y, X = patsy.dmatrices('y ~ I(x0 + x1)', data) In [50]: X Out[50]: DesignMatrix with shape (5, 2)  Intercept  I(x0 + x1)          1        1.01          1        1.99          1        3.25          1       -0.10          1        5.00  Terms:    'Intercept' (column 0)    'I(x0 + x1)' (column 1)

Patsy的patsy.builtins模块还有一些其它的内置转换。请查看线上文档。

分类数据有一个特殊的转换类，下面进行讲解。

分类数据和Patsy

非数值数据可以用多种方式转换为模型设计矩阵。完整的讲解超出了本书范围，最好和统计课一起学习。

当你在Patsy公式中使用非数值数据，它们会默认转换为虚变量。如果有截距，会去掉一个，避免共线性：

In [51]: data = pd.DataFrame({   ....:     'key1': ['a', 'a', 'b', 'b', 'a', 'b', 'a', 'b'],   ....:     'key2': [0, 1, 0, 1, 0, 1, 0, 0],   ....:     'v1': [1, 2, 3, 4, 5, 6, 7, 8],   ....:     'v2': [-1, 0, 2.5, -0.5, 4.0, -1.2, 0.2, -1.7]   ....: }) In [52]: y, X = patsy.dmatrices('v2 ~ key1', data) In [53]: X Out[53]: DesignMatrix with shape (8, 2)  Intercept  key1[T.b]          1          0          1          0          1          1          1          1          1          0          1          1          1          0          1          1  Terms:    'Intercept' (column 0)    'key1' (column 1)

如果你从模型中忽略截距，每个分类值得列都会包括在设计矩阵的模型中：

In [54]: y, X = patsy.dmatrices('v2 ~ key1 + 0', data) In [55]: X Out[55]: DesignMatrix with shape (8, 2)  key1[a]  key1[b]        1        0        1        0        0        1        0        1        1        0        0        1        1        0        0        1  Terms:    'key1' (columns 0:2)

使用C函数，数值列可以截取为分类量：

In [56]: y, X = patsy.dmatrices('v2 ~ C(key2)', data) In [57]: X Out[57]: DesignMatrix with shape (8, 2)  Intercept  C(key2)[T.1]          1             0          1             1          1             0          1             1          1             0          1             1          1             0          1             0  Terms:    'Intercept' (column 0)    'C(key2)' (column 1)

当你在模型中使用多个分类名，事情就会变复杂，因为会包括key1:key2形式的相交部分，它可以用在方差（ANOVA）模型分析中：

In [58]: data['key2'] = data['key2'].map({0: 'zero', 1: 'one'}) In [59]: data Out[59]:  key1  key2  v1   v2 0    a  zero   1 -1.0 1    a   one   2  0.0 2    b  zero   3  2.5 3    b   one   4 -0.5 4    a  zero   5  4.0 5    b   one   6 -1.2 6    a  zero   7  0.2 7    b  zero   8 -1.7 In [60]: y, X = patsy.dmatrices('v2 ~ key1 + key2', data) In [61]: X Out[61]: DesignMatrix with shape (8, 3)  Intercept  key1[T.b]  key2[T.zero]          1          0             1          1          0             0          1          1             1          1          1             0          1          0             1          1          1             0          1          0             1          1          1             1  Terms:    'Intercept' (column 0)    'key1' (column 1)    'key2' (column 2) In [62]: y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data) In [63]: X Out[63]: DesignMatrix with shape (8, 4)  Intercept  key1[T.b]  key2[T.zero] key1[T.b]:key2[T.zero]          1          0             1                       0          1          0             0                       0          1          1             1                       1          1          1             0                       0          1          0             1                       0          1          1             0                       0          1          0             1                       0          1          1             1                       1  Terms:    'Intercept' (column 0)    'key1' (column 1)    'key2' (column 2)    'key1:key2' (column 3)

Patsy提供转换分类数据的其它方法，包括以特定顺序转换。请参阅线上文档。

赞赏作者

好课推荐，名额有限，下图扫码即可免费学

Python爱好者社区历史文章大合集：

福利：文末扫码立刻关注公众号，“Python爱好者社区”，开始学习Python课程：

关注后在公众号内回复“课程”即可获取：

小编的Python入门视频课程！！！

崔老师爬虫实战案例免费学习视频。

丘老师数据科学入门指导免费学习视频。

陈老师数据分析报告制作免费学习视频。

玩转大数据分析！Spark2.X+Python 精华实战课程免费学习视频。

丘老师Python网络爬虫实战免费学习视频。