• 13.2 用Patsy创建模型描述
    • 用Patsy公式进行数据转换
    • 分类数据和Patsy

    13.2 用Patsy创建模型描述

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

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

    1. y ~ x0 + x1

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

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

    现在有:

    1. In [33]: y
    2. Out[33]:
    3. DesignMatrix with shape (5, 1)
    4. y
    5. -1.5
    6. 0.0
    7. 3.6
    8. 1.3
    9. -2.0
    10. Terms:
    11. 'y' (column 0)
    12. In [34]: X
    13. Out[34]:
    14. DesignMatrix with shape (5, 3)
    15. Intercept x0 x1
    16. 1 1 0.01
    17. 1 2 -0.01
    18. 1 3 0.25
    19. 1 4 -4.10
    20. 1 5 0.00
    21. Terms:
    22. 'Intercept' (column 0)
    23. 'x0' (column 1)
    24. 'x1' (column 2)

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

    1. In [35]: np.asarray(y)
    2. Out[35]:
    3. array([[-1.5],
    4. [ 0. ],
    5. [ 3.6],
    6. [ 1.3],
    7. [-2. ]])
    8. In [36]: np.asarray(X)
    9. Out[36]:
    10. array([[ 1. , 1. , 0.01],
    11. [ 1. , 2. , -0.01],
    12. [ 1. , 3. , 0.25],
    13. [ 1. , 4. , -4.1 ],
    14. [ 1. , 5. , 0. ]])

    你可能想Intercept是哪里来的。这是线性模型(比如普通最小二乘回归)的惯例用法。添加 +0 到模型可以不显示intercept:

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

    Patsy对象可以直接传递到算法(比如numpy.linalg.lstsq)中,它执行普通最小二乘回归:

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

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

    1. In [39]: coef
    2. Out[39]:
    3. array([[ 0.3129],
    4. [-0.0791],
    5. [-0.2655]])
    6. In [40]: coef = pd.Series(coef.squeeze(), index=X.design_info.column_names)
    7. In [41]: coef
    8. Out[41]:
    9. Intercept 0.312910
    10. x0 -0.079106
    11. x1 -0.265464
    12. dtype: float64

    用Patsy公式进行数据转换

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

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

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

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

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

    patsy.build_design_matrices函数可以使用原始样本数据集的保存信息,来转换新数据,:

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

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

    1. In [49]: y, X = patsy.dmatrices('y ~ I(x0 + x1)', data)
    2. In [50]: X
    3. Out[50]:
    4. DesignMatrix with shape (5, 2)
    5. Intercept I(x0 + x1)
    6. 1 1.01
    7. 1 1.99
    8. 1 3.25
    9. 1 -0.10
    10. 1 5.00
    11. Terms:
    12. 'Intercept' (column 0)
    13. 'I(x0 + x1)' (column 1)

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

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

    分类数据和Patsy

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

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

    1. In [51]: data = pd.DataFrame({
    2. ....: 'key1': ['a', 'a', 'b', 'b', 'a', 'b', 'a', 'b'],
    3. ....: 'key2': [0, 1, 0, 1, 0, 1, 0, 0],
    4. ....: 'v1': [1, 2, 3, 4, 5, 6, 7, 8],
    5. ....: 'v2': [-1, 0, 2.5, -0.5, 4.0, -1.2, 0.2, -1.7]
    6. ....: })
    7. In [52]: y, X = patsy.dmatrices('v2 ~ key1', data)
    8. In [53]: X
    9. Out[53]:
    10. DesignMatrix with shape (8, 2)
    11. Intercept key1[T.b]
    12. 1 0
    13. 1 0
    14. 1 1
    15. 1 1
    16. 1 0
    17. 1 1
    18. 1 0
    19. 1 1
    20. Terms:
    21. 'Intercept' (column 0)
    22. 'key1' (column 1)

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

    1. In [54]: y, X = patsy.dmatrices('v2 ~ key1 + 0', data)
    2. In [55]: X
    3. Out[55]:
    4. DesignMatrix with shape (8, 2)
    5. key1[a] key1[b]
    6. 1 0
    7. 1 0
    8. 0 1
    9. 0 1
    10. 1 0
    11. 0 1
    12. 1 0
    13. 0 1
    14. Terms:
    15. 'key1' (columns 0:2)

    使用C函数,数值列可以截取为分类量:

    1. In [56]: y, X = patsy.dmatrices('v2 ~ C(key2)', data)
    2. In [57]: X
    3. Out[57]:
    4. DesignMatrix with shape (8, 2)
    5. Intercept C(key2)[T.1]
    6. 1 0
    7. 1 1
    8. 1 0
    9. 1 1
    10. 1 0
    11. 1 1
    12. 1 0
    13. 1 0
    14. Terms:
    15. 'Intercept' (column 0)
    16. 'C(key2)' (column 1)

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

    1. In [58]: data['key2'] = data['key2'].map({0: 'zero', 1: 'one'})
    2. In [59]: data
    3. Out[59]:
    4. key1 key2 v1 v2
    5. 0 a zero 1 -1.0
    6. 1 a one 2 0.0
    7. 2 b zero 3 2.5
    8. 3 b one 4 -0.5
    9. 4 a zero 5 4.0
    10. 5 b one 6 -1.2
    11. 6 a zero 7 0.2
    12. 7 b zero 8 -1.7
    13. In [60]: y, X = patsy.dmatrices('v2 ~ key1 + key2', data)
    14. In [61]: X
    15. Out[61]:
    16. DesignMatrix with shape (8, 3)
    17. Intercept key1[T.b] key2[T.zero]
    18. 1 0 1
    19. 1 0 0
    20. 1 1 1
    21. 1 1 0
    22. 1 0 1
    23. 1 1 0
    24. 1 0 1
    25. 1 1 1
    26. Terms:
    27. 'Intercept' (column 0)
    28. 'key1' (column 1)
    29. 'key2' (column 2)
    30. In [62]: y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data)
    31. In [63]: X
    32. Out[63]:
    33. DesignMatrix with shape (8, 4)
    34. Intercept key1[T.b] key2[T.zero]
    35. key1[T.b]:key2[T.zero]
    36. 1 0 1 0
    37. 1 0 0 0
    38. 1 1 1 1
    39. 1 1 0 0
    40. 1 0 1 0
    41. 1 1 0 0
    42. 1 0 1 0
    43. 1 1 1 1
    44. Terms:
    45. 'Intercept' (column 0)
    46. 'key1' (column 1)
    47. 'key2' (column 2)
    48. 'key1:key2' (column 3)

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