Question

7.2.1　描述性统计

下面先来分析葡萄酒质量数据集。首先，计算出每列的总体描述性统计量、质量列中的唯一值以及和这个唯一值对应的观测数量。要完成这个任务，需要创建一个新脚本 wine_quality.py，然后输入以下代码：

 1 #!/usr/bin/env python3
 2 import numpy as np
 3 import pandas as pd
 4 import seaborn as sns
 5 import matplotlib.pyplot as plt
 6 import statsmodels.api as sm
 7 import statsmodels.formula.api as smf
 8 from statsmodels.formula.api import ols, glm
 9 # 将数据集读入到pandas数据框中
10 wine = pd.read_csv('winequality-both.csv', sep=',', header=0)
11 wine.columns = wine.columns.str.replace(' ', '_')
12 print(wine.head())
13 # 显示所有变量的描述性统计量
14 print(wine.describe())
15 # 找出唯一值
16 print(sorted(wine.quality.unique()))
17 # 计算值的频率
18 print(wine.quality.value_counts())

在 import 语句之后，要做的第一件事就是使用 pandas 的 read_csv 方法将文本文件 winequality-both.csv 读入一个 pandas 数据框。附加参数表示域分隔符为逗号，第一行为列标题。有些列标题中包含空格（例如：fixed acidity），所以在下面一行代码中要使用下划线替换空格。然后，使用 head 函数检查一下标题行和前 5 行数据，确保数据被正确加载。

第 14 行代码使用 pandas 的 describe 函数打印出数据集中每个数值型变量的摘要统计量。这些统计量包括总数、均值、标准差、最小值、第 25 个百分位数、中位数、第 75 个百分位数和最大值。例如，质量评分中有 6497 个观测，评分范围从 3 到 9，平均质量评分为 5.8，标准差为 0.87。

下面一行代码识别出质量列中的唯一值，并以升序打印在屏幕上。输出显示质量列中的唯一值是 3、4、5、6、7、8 和 9。

最终，本节最后一行代码计算出质量列中每个唯一值在数据集中出现的次数，并把它们以降序打印到屏幕上。输出显示，有 2836 个观测的质量评分为 6，2138 个观测的质量评分为 5，1079 个观测的质量评分为 7，216 个观测的质量评分为 4，193 个观测的质量评分为 8，30 个观测的质量评分为 3，5 个观测的质量评分为 9。

7.2.2　分组、直方图与t 检验

前面计算出的统计量是针对整个数据集的，既包括红葡萄酒数据也包括白葡萄酒数据。下面我们分别分析红葡萄酒数据和白葡萄酒数据，看看统计量是否会保持不变：

# 按照葡萄酒类型显示质量的描述性统计量
print(wine.groupby('type')[['quality']].describe().unstack('type'))
# 按照葡萄酒类型显示质量的特定分位数值
print(wine.groupby('type')[['quality']].quantile([0.25, 0.75]).unstack('type'))
# 按照葡萄酒类型查看质量分布
red_wine = wine.loc[wine['type']=='red', 'quality']
white_wine = wine.loc[wine['type']=='white', 'quality']
sns.set_style("dark")
print(sns.distplot(red_wine, \
         norm_hist=True, kde=False, color="red", label="Red wine"))
print(sns.distplot(white_wine, \
         norm_hist=True, kde=False, color="white", label="White wine"))
sns.axlabel("Quality Score", "Density")
plt.title("Distribution of Quality by Wine Type")
plt.legend()
plt.show()
# 检验红葡萄酒和白葡萄酒的平均质量是否有所不同
print(wine.groupby(['type'])[['quality']].agg(['std']))
tstat, pvalue, df = sm.stats.ttest_ind(red_wine, white_wine)
print('tstat: %.3f pvalue: %.4f' % (tstat, pvalue))

本节的第一行代码分别打印出红葡萄酒和白葡萄酒的摘要统计量。groupby 函数使用 type 列中的两个值将数据分为两组。方括号可以生成一个列表，列表中的元素是用来生成输出的列。在本例中，只对质量列应用 describe 函数。这些命令的结果就是生成一列统计量，来自红葡萄酒数据的计算结果和白葡萄酒数据的计算结果是相互垂直地堆叠在一起的。unstack 函数将结果重新排列，这样红葡萄酒和白葡萄酒的统计量就会显示在并排的两列中。

下一行代码和它前面的代码非常相似，但不是使用 describe 函数计算一些描述性统计量，而是使用 quantile 函数对质量列计算第 25 百分位数和第 75 百分位数。

接下来，使用 seaborn 创建一幅统计图，其中有两个直方图，一个表示红葡萄酒，另一个表示白葡萄酒，如图 7-3 所示。

图 7-3：表示两种类型葡萄酒评分分布的密度直方图

红条表示红葡萄酒，白条表示白葡萄酒。因为白葡萄酒数据比红葡萄酒多（白葡萄酒有 4898 条数据，红葡萄酒有 1599 条数据），所以图中显示密度分布，不显示频率分布。从这个统计图可以看出，两种葡萄酒的评分都近似正态分布。与原始数据的摘要统计量相比，直方图更容易看出两种葡萄酒的质量评分的分布。

最后，进行一下 t 检验，判断红葡萄酒和白葡萄酒的平均评分是否有区别。这段代码演示了如何使用 groupby 和 agg 函数来为数据集中的分组计算一系列统计量。在本例中，我们想知道红葡萄酒和白葡萄酒评分的标准差是否相同，所以在 t 检验中可以使用合并方差。t 检验统计量为 -9.69，p 值为 0.00，这说明白葡萄酒的平均质量评分在统计意义上大于红葡萄酒的平均质量评分。

7.2.3　成对变量之间的关系和相关性

前面已经检查了输出变量，下面简单研究一下输入变量。让我们计算一下输入变量两两之间的相关性，并为一些输入变量创建带有回归直线的散点图：

# 计算所有变量的相关矩阵
print(wine.corr())
# 从红葡萄酒和白葡萄酒的数据中取出一个“小”样本来进行绘图
def take_sample(data_frame, replace=False, n=200):
    return data_frame.loc[np.random.choice(data_frame.index, \
    replace=replace, size=n)]
reds_sample = take_sample(wine.loc[wine['type']=='red', :])
whites_sample = take_sample(wine.loc[wine['type']=='white', :])
wine_sample = pd.concat([reds_sample, whites_sample])
wine['in_sample'] = np.where(wine.index.isin(wine_sample.index), 1.,0.)
print(pd.crosstab(wine.in_sample, wine.type, margins=True))
# 查看成对变量之间的关系
sns.set_style("dark")
g = sns.pairplot(wine_sample, kind='reg', plot_kws={"ci": False,\
"x_jitter": 0.25, "y_jitter": 0.25}, hue='type', diag_kind='hist',\
diag_kws={"bins": 10, "alpha": 1.0}, palette=dict(red="red", white="white"),\
markers=["o", "s"], vars=['quality', 'alcohol', 'residual_sugar'])
print(g)
plt.suptitle('Histograms and Scatter Plots of Quality, Alcohol, and Residual\
Sugar', fontsize=14, horizontalalignment='center', verticalalignment='top',\
x=0.5, y=0.999)
plt.show()

corr 函数可以计算出数据集中所有变量两两之间的线性相关性。根据相关系数的符号，从输出中可以知道酒精含量、硫酸盐、pH 值、游离二氧化硫和柠檬酸这些指标与质量是正相关的，相反，非挥发性酸、挥发性酸、残余糖分、氯化物、总二氧化硫和密度这些指标与质量是负相关的。

数据集中有 6000 多个点，所以如果将它们都画在统计图中，就很难分辨出清楚的点。为解决这个问题，我们定义了一个函数 take_sample，用来抽取在统计图中使用的样本点。这个函数使用 pandas 数据框索引和 numpy 的 random.choice 函数随机选择一个行的子集。我们使用这个函数对红葡萄酒和白葡萄酒分别进行抽样，并将抽样所得的两个数据框连接成一个数据框。然后，在 wine 数据框中创建一个新列 in_sample，并使用 numpy 的 where 函数和 pandas 的 isin 函数对这个新列进行填充，填充的值根据此行的索引值是否在抽样数据的索引值中分别设为 1 和 0。最后，我们使用 pandas 的 crosstab 函数来确认 in_sample 列中包含 400 个 1（200 条红葡萄酒数据和 200 条白葡萄酒数据）和 6097 个 0。

seaborn 的 pairplot 函数可以创建一个统计图矩阵。主对角线上的图以直方图或密度图的形式显示了每个变量的单变量分布，对角线之外的图以散点图的形式显示了每两个变量之间的双变量分布，散点图中可以有回归直线，也可以没有。

图 7-4 中的统计图显示了葡萄酒质量、酒精含量和残余糖分之间的关系。红条和红点表示红葡萄酒，白条和白点表示白葡萄酒。因为质量评分都是整数，所以我加上了一点振动，这样更容易看出数据在何处集中。

图 7-4：质量、酒精含量和残余糖分这 3 个变量两两之间的散点图、回归直线和直方图，按葡萄酒类型分类

从这些统计图可以看出，对于红葡萄酒和白葡萄酒来说，酒精含量的均值和标准差是大致相同的，但是，白葡萄酒残余糖分的均值和标准差却大于红葡萄酒残余糖分的均值和标准差。从回归直线可以看出，对于两种类型的葡萄酒，酒精含量增加时，质量评分也随之提高，相反，残余糖分增加时，质量评分则随之降低。这两个变量对白葡萄酒的影响都要大于对红葡萄酒的影响。

7.2.4　使用最小二乘估计进行线性回归

相关系数和两两变量之间的统计图有助于对两个变量之间的关系进行量化和可视化，但是它们不能测量出每个自变量在其他自变量不变时与因变量之间的关系。线性回归可以解决这个问题。

线性回归模型如下：

· ,

对于 i = 1, 2,…, n 个观测和 p 个自变量。

这个模型表示观测 y_i 服从均值为 μ_i 方差为 σ² 的正态分布（高斯分布），其中 μ_i 依赖于自变量，σ² 为一个常数。也就是说，给定了自变量的值之后，我们就可以得到一个具体的质量评分，但在另一天，给定同样的自变量值，我们可能会得到一个和前面不同的质量评分。但是，经过很多天自变量取同样的值（也就是一个很长的周期），质量评分会落在 μ_i±σ 这个范围内。

理解了线性回归模型后，下面让我们使用 statsmodel 包来进行线性回归：

my_formula = 'quality ~ alcohol + chlorides + citric_acid + density\
+ fixed_acidity + free_sulfur_dioxide + pH + residual_sugar + sulphates\
+ total_sulfur_dioxide + volatile_acidity'

lm = ols(my_formula, data=wine).fit()

## 或者，也可以使用广义线性模型（glm）语法进行线性回归
## lm = glm(my_formula, data=wine, family=sm.families.Gaussian()).fit()

print(lm.summary())
print("\nQuantities you can extract from the result:\n%s" % dir(lm))
print("\nCoefficients:\n%s" % lm.params)
print("\nCoefficient Std Errors:\n%s" % lm.bse)
print("\nAdj. R-squared:\n%.2f" % lm.rsquared_adj)
print("\nF-statistic: %.1f   P-value: %.2f" % (lm.fvalue, lm.f_pvalue))
print("\nNumber of obs: %d   Number of fitted values: %d" % (lm.nobs,\
len(lm.fittedvalues)))

第一行代码将一个字符串赋给变量 my_foumula。这个字符串中包含的是类似 R 语言语法的回归公式定义。波浪线（~）左侧的变量 quality 是因变量，波浪线右侧的变量是自变量。

第二行代码使用公式和数据拟合一个普通最小二乘回归模型，并将结果赋给变量 lm。为了演示另外一种拟合方法，下一行代码（注释中的代码）使用广义线性模型（glm）的语法代替普通最小二乘语法，拟合同样的模型。

下面 7 行代码向屏幕上打印模型结果。第一行向屏幕上打印结果的摘要信息。这些摘要信息非常重要，因为它包含了模型系数、系数的标准差和置信区间、修正 R 方、F 统计量等模型详细信息。

下一行代码打印出一个列表，其中包含从模型对象 lm 中提取出的所有数值信息，检查了这个列表之后，我希望提取出模型系数、系数的标准差、修正 R 方、F 统计量和它的 p 值，以及模型拟合值。

接下来的 4 行代码提取出了这些值。lm.params 以一个序列的形式返回模型系数，这样你可以通过定位或名称提取出单个的系数。例如，要提取酒精含量的系数 0.267，你可以使用 lm.params[1] 或 lm.params['alcohol']。同样，lm.bse 以序列的形式返回模型系数的标准差。lm.rsquared_adj 返回修正 R 方，lm.fvalue 和 lm.f_pvalue 分别返回 F 统计量和它的 p 值。最后，lm.fittedvalues 返回拟合值。我没有给出所有的拟合值，而是在观测总数 lm.nobs 后面显示拟合值的数量，以确认拟合值与观测具有同样的数量。输出结果如图 7-5 所示。

图 7-5：具有 11 个特征的葡萄酒质量多元线性回归

7.2.5　系数解释

如果你想使用这个模型弄清楚因变量（葡萄酒质量）和自变量（11 个葡萄酒特性）之间的关系，就应该解释一下模型系数的意义。在这个模型中，某个自变量系数的意义是，在其他自变量保持不变的情况下，这个自变量发生 1 个单位的变化时，导致葡萄酒质量评分发生的平均变化。例如，酒精含量系数的含义就是，从平均意义上来说，如果两种葡萄酒其他自变量的值都相同，那么酒精含量高 1 个单位的葡萄酒的质量评分就要比另一种葡萄酒的质量评分高出 0.27 分。

并不是所有的系数都需要解释。例如，截距系数的意义是当所有自变量的值都为 0 时的期望评分。因为没有任何一种葡萄酒的各种成分都为 0，所以截距系数没有具体意义。

7.2.6　自变量标准化

关于这个模型，还需要注意的一点是，普通最小二乘回归是通过使残差平方和最小化来估计未知的 β 参数值的，这里的残差是指自变量观测值与拟合值之间的差别。因为残差大小是依赖于自变量的测量单位的，所以如果自变量的测量单位相差很大的话，那么将自变量标准化后，就可以更容易对模型进行解释了。对自变量进行标准化的方法是，先从自变量的每个观测值中减去均值，然后再除以这个自变量的标准差。自变量标准化完成以后，它的均值为 0，标准差为 11。

1在 Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, 2007, p.56) 中，Gelman and Hill 建议在数据集中既有连续型自变量又有二值型自变量的情况下，用两倍标准差去除，而不是用一倍标准差。这样的话，标准化自变量一个单位的变化就对应于均值上下一个标准差的变化。因为葡萄酒质量数据集中不包含二值自变量，所以我通过除以一倍标准差将自变量标准化为 z-scores。

使用 wine.describe()，可以看到氯化物的范围是从 0.009 到 0.661，而总二氧化硫的范围是从 6.0 到 440.0。其余各变量的最小值与最大值之间的区别也大致如此。因为各个自变量值的变化范围相差非常悬殊，所以非常应该对自变量进行标准化，看看这样做了之后，能否更容易对结果进行解释。

pandas 在数据框中对变量进行标准化非常容易。你可以对一个观测写一个变换公式，pandas 可以把这个公式扩展到行与列中，来标准化所有变量。以下各行代码使用标准化后的自变量创建了一个新数据框 wine_standardized：

# 创建一个名为dependent_variable的序列来保存质量数据
dependent_variable = wine['quality']

# 创建一个名为independent variables的数据框
# 来保存初始的葡萄酒数据集中除quality、type和in_sample之外的所有变量
independent_variables = wine[wine.columns.difference(['quality', 'type',\
'in_sample'])]

# 对自变量进行标准化
# 对每个变量，在每个观测中减去变量的均值
# 并且使用结果除以变量的标准差
independent_variables_standardized = (independent_variables -\
independent_variables.mean()) / independent_variables.std()

# 将因变量quality作为一列添加到自变量数据框中
# 创建一个带有标准化自变量的
# 新数据集
wine_standardized = pd.concat([dependent_variable, independent_variables\
_standardized], axis=1)

完成了带有标准化自变量的数据集之后，让我们重新进行线性回归，并查看一下摘要统计（见图 7-6）：

lm_standardized = ols(my_formula, data=wine_standardized).fit()
print(lm_standardized.summary())

图 7-6：具有 11 个特征的葡萄酒质量多元线性回归——在回归之前，自变量被标准化为 z-scores

自变量标准化会改变我们对模型系数的解释。现在每个自变量系数的含义是，不同的葡萄酒在其他自变量均相同的情况下，某个自变量相差 1 个标准差，会使葡萄酒的质量评分平均相差多少个标准差。举个例子，酒精含量系数的意义是，从平均意义上说，如果两种葡萄酒其他自变量的值都相同，那么酒精含量高 1 个标准差的葡萄酒的质量评分就要比另一种葡萄酒的质量评分高出 0.32 个标准差。

还是通过 wine.describe() 函数，我们可以看到酒精含量的均值和标准差是 10.5 和 1.2，质量评分的均值和标准差是 5.8 和 0.9。因此，从平均意义上说，如果两种葡萄酒其他的自变量值均相同，那么酒精含量为 11.7（10.5+1.2）的葡萄酒的质量评分就会比酒精含量为均值的葡萄酒的质量评分大 0.32 个标准差。

自变量标准化同样会改变我们对截距的解释。当解释变量被标准化后，截距表示的就是当所有自变量取值为均值时因变量的均值。在我们的模型摘要中，截距系数的意义就是当一种葡萄酒所有的成分都取均值的时候，它的质量评分均值应该是 5.8，标准差为 0.009。

7.2.7　预测

在某些情况下，我们需要使用没有用来拟合模型的新数据进行预测。例如，你会收到关于葡萄酒成分的一个新的观测，并需要根据这些成分预测这种葡萄酒的质量评分。让我们通过选择现有数据集的前 10 个观测，并根据它们的葡萄酒成分预测质量评分，来演示一下如何对新数据做出预测。

需要说明的是，仅出于方便和演示的目的，我们才使用这些已经用于拟合模型的数据。除了这个示例之外，你应该使用未用于拟合模型的数据来评价模型，并用新的观测数据进行预测。记住了这一点之后，下面创建一组“新”观测，并使用它们来预测质量评分：

# 使用葡萄酒数据集中的前10个观测创建10个“新”观测
# 新观测中只包含模型中使用的自变量
new_observations = wine.ix[wine.index.isin(range(10)), \
independent_variables.columns]

# 基于新观测中的葡萄酒特性预测质量评分
y_predicted = lm.predict(new_observations)

# 将预测值保留两位小数并打印到屏幕上
y_predicted_rounded = [round(score, 2) for score in y_predicted]
print(y_predicted_rounded)

变量 y_predicted 中包含着 10 个预测值。为了使输出更简单易懂，我将预测值保留两位小数。在这个示例中，如果我们使用的是真正的新观测，就可以使用这些预测值来评价模型。无论如何，我们得到了一些预测值，可以用来评估或做些别的事情。