[学习笔记]Python for Data Analysis, 3E-12.Python建模库简介

在本书中，重点聚焦于提供在Python中进行数据分析的编程基础。由于数据分析时和科学家经常报告在数据整理和准备上话费不成比例的时间，因此本书的结构反应了掌握这些技术的重要性。
用于开发模型的库将取决于应用程序。许多统计问题可以通过更简单的技术来解决，如普通最小二乘回归，而其他问题可能需要更先进的机器学习方法。幸运的是，Python已经成为实现分析方法的首选语言之一，因此在完成本书之后，你可以探索许多工具。
在本章中，我们将回顾pandas的一些功能，当你在与pandas的数据整理、模型拟合和评分之间来回切换时候，这些功能可能会有所帮助。我将简要介绍两大流行的建模工具包，statsmodels和scikit-learn。

12.1模型代码和pandas之间的接口

模型开发的常见工作流程是使用pandas进行数据加载和清理，然后再切换到建模库来构建模型本身。模型开发过程的一个重要部分称为机器学习中的特征工程。这可以描述从原始数据集中提取在建模上下文中可能有用的信息的任何数据转换和分析。我们在本书中探讨的数据聚合和GroupBy工具经常在特征工程环境中使用。
虽然’好’的特征工程的细节超出了本书的范围，但我将展示一些方法，是pandas的数据操作和建模之间的切换尽可能轻松。
pandas和其他分析库之间的接触点通常是NumPy数组。若要将DataFrame转换为NumPy数组，请使用to_numpy方法：

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.]})
data.to_numpy() # 将DataFrame转化为NumPy数组(二维数组)
# 要转换回DataFrame，可以使用可选的列名传递二维数组
df2 = pd.DataFrame(data.to_numpy(), columns=['one', 'two', 'three'])

# to_numpy方法适用于数据是同构(例如，全是数值类型)的情况。
# 如果有异构数据，结果将是object组成的的ndarray
df3 = data.copy()
df3['strings'] = ['a', 'b', 'c', 'd', 'e']
df3.to_numpy()

# 对于某些模型，你可能希望仅使用列的子集。此时推荐将loc索引和to_numpy结合
model_cols = ['x0', 'x1']
data.loc[:, model_cols].to_numpy() # 先筛选列，再转化为数组

某些库具有对pandas的原生支持，并自动为你执行其中的一些工作：从DataFrame转换为NumPy，并将模型参数名称附加到输出表或者Series的列。在其他情况下，你必须手动执行此’元数据管理’。

在’章节7.5：分类数据’中，我们研究了pandas的Categorical类型和pandas.get_dummies函数。

# 假设我们的示例数据集中有一个非数字列
data['category'] = pd.Categorical(['a', 'b', 'a', 'a', 'b'], categories=['a', 'b'])

# 如果我们想用虚拟变量替换列，则我们创建虚拟变量，删除'category'列，然后连接结果：
dummies = pd.get_dummies(data.category, prefix='category') # 生成含'category_a'和'category_b'列的DataFrame
date_with_dummies = data.drop('category', axis=1).join(dummies) # 在列中删除'category'列，并拼接dummies

使用虚拟变量你和某些统计模型有一些细微差别。当你拥有的不仅仅是简单的数字列时，使用Patsy(下一节的主题)可能更简单且不容易出错。

12.2使用Patsy创建模型描述

Patsy是一个Python库，用于使用基于字符串的’公式语法’描述统计模型(尤其是线性模型)，该语法的灵感来自R和S统计编程语言使用的公式语法(但不完全相同)。当你安装stasmodels时候，Patsy库会自动安装。
Patsy在statsmodels中指定线性模型得到了很好的支持，因此我将重点介绍一些主要功能来帮你启动和运行。Patsy的公式是一种特殊的字符串语法，如下所示：

y ~ x0 + x1

a+b的语法并不意味着a加b，而是为模型创建的在设计矩阵中的项。patsy.dmatrices函数将公式字符串与数据集(可以是DataFrame或数组字典)一起获取，并为线性模型生成设计矩阵：

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.]})
import patsy
y, X = patsy.dmatrices('y ~ x0 + x1', data) # y为DesignMatrix with shape(5, 1)，项为'y'，X为DesignMatrix with shape(5, 3)，项为'截距，x0, x1'

# 这些Patsy DesignMatrix实例是具有附加元数据的NumPy数组：
np.asarray(y) # 转化为NumPy数组(元数据被去掉)
np.asarray(X)

你也许想知道Intercept项是从哪里来的。这是线性模型，如普通最小二乘(OLS)回归的约定。你可以通过向模型添加’+0’项来取消截距：

patsy.dmatrices('y ~ x0 + x1 + 0', data)[1]

# Patsy对象可以直接传递到类似numpy.linalg.lstsq的算法中，该算法执行普通最小二乘回归
coef, resid, _, _ = np.linalg.lstsq(X, y) # 得到系数等结果

# 模型的元数据保留在design_info属性中，因此你可以将模型列名称重新附加到拟合系数以获取Series，例如
coef = pd.Series(coef.squeeze(), index=X.design_info.column_names) # 构建索引为'Intercept', 'x0'和'x1'的Series

Patsy公式中的数据转换

你可以将Python代码混合到Patsy公式中；当计算公式时，库将尝试在封闭作用于中找到你使用的函数：

y, X = patsy.dmatrices('y ~ x0 + np.log(np.abs(x1) + 1)', data) # 设计矩阵'X'中的项为'x0'和'np.log(np.abs(x1)+1'

一些常用的变量变换包括标准化(转化为均值0和方差1)和中心化(减去均值)。Pasty具有用于此目的的内置函数：

y, X = patsy.dmatrices('y ~ standardize(x0) + center(x1)', data)

作为建模过程的一部分，你可以在一个数据集上拟合模型，然后基于另一个数据集评估模型。这可能是原数据的保留部分或者稍后观察到的新数据。应用中心化和标准化等转换时，使用模型基于新数据形成预测时应小心。这些成为有状态转换，因为在转换新数据集时必须使用原始数据集的平均值或标准差等统计信息。
patsy.build_design_matrices函数可以使用原始样本内数据集中保存的信息将转换应用于新的’样本外数据’：

new_data = pd.DataFrame({'x0': [6, 7, 8, 9], 'x1':[3.1, -0.5, 0, 2.3], 'y':[1, 2, 3, 4]}) # 创建新的DataFrame
new_X = patsy.build_design_matrices([X.design_info], new_data) # 将原始样本的转换应用到新的数据中，得到DesignMatrix对象

由于Patsy公式上下文的加号(+)并不意味着加法，因此当你要按照名称添加数据集中的列时，必须将它们包装在特定的I(大写的i)函数中：

y, X = patsy.dmtrices('y ~ I(x0 + x1)', data) # I(x0 + x1)是一个完整的项

Patsy在patsy.builtins模块中还有其他几个内置转换。有关详细信息可以参阅在线文档。
分类数据具有一类特殊的转换，将在下面介绍。

分类数据和Patsy

非数值数据可以通过很多种方式转化为模型设计矩阵。对这个主题的完整处理超出了本书的范围，最好与统计学课程一起研究。
在Patsy公式中使用非数值项时，默认情况下，它们将转化为虚拟变量。如果存在截距，则将省略其中一个水平以避免共线性：

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]})
y, X = patsy.dmatrices('v2 ~ key1', data) # 因为存在截距(没有+0，默认存在)，所以这里省略了'a'对应的水平key1[T.a]，只有key1[T.b]

# 如果从模型中省略截距，则每个类别值的列将包含在模型设计矩阵中
y, X = patsy.dmatrices('v2 ~ key1 + 0', data) # 由于没有截距，所以存在水平key1[a]和key1[b]

# 数值列可以使用C函数解释为分类变量
y, X = patsy.dmatrices('v2 ~ C(key2)', data) # 包含截距项和数值1对应的水平C(key2)[T.1]

当你在模型中使用多个分类项时，事情可能会更加复杂，因为你可以包含形式的交互作用项key1:key2，例如，可以在方差分析(ANOVA)模型中使用它

data['key2'] = data['key2'].map({0: 'zero', 1: 'one'}) # 将数值型变量转化为非数值型变量(分类变量)
y, X = patsy.dmatrices('v2 ~ key1 + key2', data) # DesignMatrix X中含三个项：Intercept，key1[T.b], key2[T.zero]
y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data) # DesignMatrix X中含四个项：截距，key1[T.b], key2[T.zero], key1[T.b]:key2[T.zero]

12.3统计模型简介

statsmodels是一个用于拟合各种统计模型，执行统计测试以及数据探索和可视化的Python库。statsmodels包含更多’经典’频率统计方法，而贝叶斯方法和机器学习模型可以在别的库找到。

在statsmodels找到的一些模型包括：

• 线性模型，生成线性模型和鲁棒线性模型
• 线性混合效应模型
• 方差分析方法(ANOVA)
• 时序过程和状态空间模型
• 广义矩方法

在接下来的页面中，我们将使用一些statsmodels中的基础工具并探索如何通过Patsy公式和DataFrame使用模型接口。

估计线性模型

在statsmodels中有一些线性回归模型，从比较基础(如，普通最小平方)到更复杂(如，重复再加权最小平方法)。statsmodels的线性模型有两个不同的接口：基于数组和基于公式。它们通过这些API模块导入访问：

import statsmodel.api as sm
import statsmodel.formula.api as smf

为了展示如何使用这些，我们从一些随机数据生成一个线性模型。在Jupyter单元中运行以下代码:

rng = np.random.default_rng(seed=12345)

def dnorm(mean, variance, size=1):
    if isinstance(size, int):
        size = size, # 将size转化为元组
    return mean + np.sqrt(variance) * rng.standard_normal(*size) # 返回均值和方差给定的正态分布的抽样数据

N = 100
X = np.c_[dnorm(0, 0.4, size=N),
          dnorm(0, 0.6, size=N),
          dnorm(0, 0.2, size=N)] # np.c_是按列拼接，np.r_是按行拼接，这里返回100*3的数组
eps = dnorm(0, 0.1, size=N)
beta = [0.1, 0.3, 0.5]

y = np.dot(X, beta) + eps # 构造数据

这里，'真实’模型的参数为beta。在本例中，dnorm是一个辅助函数，用于生成具有特定均值和方差的正态分布数据。
线性模型通常用截距项来拟合，就像我们之前看到的Patsy一样。sm.add_constant函数可以在现有矩阵中添加一个截距列:

X_model = sm.add_constant(X) # 在X的第0列添加截距列(全为1的列)

# sm.OLS类可以拟合一个普通的最小二乘线性回归
model = sm.OLS(y, X)

# 模型的fit方法返回回归结果对象，它包含估计模型的参数和其他诊断
results = model.fit()
results.params # 返回参数组成的数组

# results上的summary方法可以打印模型的细节诊断输出
print(results.summary())

这里的参数名被赋予了通用名称x1、x2等等。

# 假设所有的模型参数都在一个数据帧中
data = pd.DataFrame(X, columns=['col0', 'col1', 'col2'])
data['y'] = y # 构造包含所有数据的DataFrame

# 现在我们使用statsmodels公式接口和Patsy公式字符串：
results = smf.ols('y ~ col0 + col1 + col2', data=data)
results.params  # 返回各参数值
results.tvalues # 返回各参数的t-value

观察statmodels是如何返回带有DataFrame列名的Series结果的。在使用公式和pandas对象时，我们也不需要使用add_constant。

# 给定新的样本外数据，可以根据估计的模型参数计算预测值:
results.predict(data[:5])

在statsmodels中还有许多用于分析、诊断和可视化线性模型结果的其他工具，您可以探索这些工具。除了普通的最小二乘，还有其他类型的线性模型。

估计时间序列过程

statsmodels中的另一类模型用于时间序列分析。其中包括自回归过程、卡尔曼滤波等状态空间模型以及多元自回归模型。
让我们用自回归结构和噪声来模拟一些时间序列数据。

init_x = 4
values = [init_x, init_x]
N = 1000
b0 = 0.8
b1 = -0.4
noise = dnorm(0, 0.1, N)
for i in range(N):
    new_x = values[-1] * b0 + values[-2] * b1 + noise[i]
    values.append(new_x)

该数据具有AR(2)结构(两个滞后)，参数为0.8和-0.4。当你拟合一个AR模型时，你可能不知道要包含的滞后项的数量，所以你可以用一些更大的滞后数来拟合模型:

from statsmodels.tsa.ar_model import AutoReg
MAXLAGS = 5
model = AutoReg(values, MAXLAGS)
results = model.fit()

结果中的估计参数先有截距，然后是前两个滞后的估计:

results.params

关于这些模型的更深入的细节以及如何解释它们的结果超出了我在本书中所能介绍的范围，但是在statmodels文档中还有更多可以发现的内容。

12.4scikit-learn介绍

scikit-learn是使用最广泛、最值得信任的通用Python机器学习工具包之一。它包含了各种标准监督和无监督机器学习方法，以及用于模型选择和评估、数据转换、数据加载和模型持久性的工具。这些模型可以用于分类、聚类、预测和其他常见任务。
关于机器学习，以及如何应用scikit-learn这样的库来解决现实世界的问题，有很好的在线和印刷资源。在本节中，我将简要介绍一下scikit-learn API格式。
近年来，在scikit-learn中，pandas的集成已经有了显著的提高，当你读到这篇文章的时候，它可能已经提高了更多。我鼓励您查看最新的项目文档。作为本章的一个例子，我使用了一个来自Kaggle竞赛的关于1912年泰坦尼克号乘客存活率的经典数据集。

# 使用pandas加载训练和测试数据集:
train = pd.read_csv('datasets/titanic/train.csv')
test = pd.read_csv('datasets/titanic/test.csv')

像statmodels和scikit-learn这样的库通常不能被提供缺失的数据，所以我们查看列，看看是否有包含缺失数据的列:

train.isna().sum()
test.isna().sum()

在像这样的统计和机器学习例子中，一个典型的任务是根据数据中的特征预测乘客是否会幸存。模型在训练数据集上拟合，然后在样本外测试数据集上进行评估。
我想使用’Age’作为一个预测器，但它有缺失的数据。有很多方法可以进行缺失数据填充，但我将使用一个简单的方法，使用训练数据集的中位数来填补两个表中的空值:

impute_value = train['Age'].mmedian()
train['Age'] = train['Age'].fillna(impute_value)
test['Age'] = test['Age'].fillna(impute_value)

现在我们需要指定我们的模型。我添加了一个’IsFemale’列作为“Sex”列的编码版本:

train['Isfemale'] = (train['Sex'] == 'female').astype(int)
test['Isfemale'] = (train['Sex'] == 'female').astype(int)

然后我们决定一些模型变量并创建NumPy数组:

predictors = ['Pclass', 'IsFemale', 'Age'] # 用于预测的变量
X_train = train[predictors].to_numpy() # 训练集
X_test = test[predictors].to_numpy()   # 测试集
y_train = train['Survived'].to_numpy()

我不会说这是一个好模型或这些特征设计得很好。我们使用scikit-learn中的LogisticRegression模型并创建一个模型实例:

from sklearn.linear_model import LogisticRegression
model = LogisticRegression() # 创建LogisticRegression模型实例

# 我们可以使用模型的fit方法将该模型与训练数据拟合:
model.fit(X_train, y_train) # 返回LogisticRegression()对象

# 现在，我们可以使用model.predict对测试数据集进行预测:
y_predict = model.predict(X_test)

# 如果你有测试数据集的真实值，你可以计算一个准确率百分比或其他误差度量:
(y_true == y_predict).mean()

在实践中，在模型训练中通常有许多额外的复杂性层。许多模型都有可以调优的参数，还有一些技术，例如交叉验证，可以用于参数调优，以避免与训练数据过拟合。这通常可以对新数据产生更好的预测性能或鲁棒性。

交叉验证通过分割训练数据来模拟样本外预测。基于像均方误差这样的模型精度评分，您可以对模型参数执行网格搜索。有些模型，如逻辑回归，具有内置交叉验证的估计器类。例如，LogisticRegressionCV类可以与一个参数一起使用，该参数指示网格搜索在模型正则化参数C上执行的细粒度:

from sklearn.linear_model import LogisticRegressionCV
model_cv = LogisticRegressionCV(Cs=10)
model_cv.fit(X_train, y_train)

要手工进行交叉验证，可以使用cross_val_score辅助函数，它处理数据拆分过程。例如，要用训练数据的四个不重叠的分割交叉验证我们的模型，我们可以这样做:

from sklearn.model_selection import cross_val_score
model = LogisticRegression(C=10)
scores = cross_val_score(model, X_train, y_train, cv=4)

默认的评分指标是依赖模型的，但是可以选择显式的评分函数。交叉验证模型需要更长的训练时间，但通常可以产生更好的模型性能。