Python数据分析基础【建模篇】

葡萄酒质量

判定红葡萄酒和白葡萄酒评分的标准差是否相同
“type”列用来区分这行数据是红葡萄酒还是白葡萄酒的数据

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.formula.api import ols, glm# 使用pandas 的read_csv 方法读入一个pandas 数据框,附加参数表示域分隔符为逗号，第一行为列标题
wine = pd.read_csv('winequality-both.csv', sep=',', header=0)
# 有些列标题中包含空格（例如：fixed acidity），所以要使用下划线替换空格
wine.columns = wine.columns.str.replace(' ', '_')
red_wine = wine.loc[wine['type']=='red', 'quality']
white_wine = wine.loc[wine['type']=='white', 'quality']
tstat, pvalue, df = sm.stats.ttest_ind(red_wine, white_wine)
print('tstat: %.3f pvalue: %.4f' % (tstat, pvalue))

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

# 计算所有变量的相关矩阵
print(wine.corr())
# 从红葡萄酒和白葡萄酒的数据中取出一个“小”样本来进行绘图
# 数据集中有6000多个点，所以如果将它们都画在统计图中，就很难分辨出清楚的点。因此定义了一个函数take_sample，用来抽取在统计图中使用的样本点
def take_sample(data_frame, replace=False, n=200):
# np.random.choice(5, 3)和np.random.randint(0,5,3)意思相同，表示从[0,5)之间随机以等概率选取3个数
#np.random.choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0])表示分别以p=[0.1, 0, 0.3, 0.6, 0]的概率从[0,1,2,3,4]这四个数中选取3个数
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，并使用numpy的where函数和pandas的isin 函数对这个新列进行填充
#填充的值根据此行的索引值是否在抽样数据的索引值中分别设为1和 0
wine['in_sample'] = np.where(wine.index.isin(wine_sample.index), 1.,0.)
# 使用pandas 的crosstab 函数来确认in_sample 列中包含400 个1（200 条红葡萄酒数据和200 条白葡萄酒数据）和6097 个0。
print(pd.crosstab(wine.in_sample, wine.type, margins=True))# 查看成对变量之间的关系
# seaborn 的pairplot 函数可以创建一个统计图矩阵。主对角线上的图以直方图或密度图的形式显示了每个变量的单变量分布，对角线之外的图以散点图的形式显示了每两个变量之间的双变量分布，散点图中可以有回归直线，也可以没有
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()

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

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

yi∼N(μi,σ2)y_{i} \sim N\left(\mu_{i}, \sigma^{2}\right)yi∼N(μi,σ2)
μi=β0+β1xi1+β2xi2+⋯+βpxip\mu_{i}=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{p} x_{i p}μi=β0+β1xi1+β2xi2+⋯+βpxip

对于i = 1, 2, …, n 个观测和p 个自变量。
这个模型表示观测yi 服从均值为μi 方差为σ2 的正态分布（高斯分布），其中μi 依赖于自变量，σ2 为一个常数。也就是说，给定了自变量的值之后，我们就可以得到一个具体的质量评分，但在另一天，给定同样的自变量值，我们可能会得到一个和前面不同的质量评分。但是，经过很多天自变量取同样的值（也就是一个很长的周期），质量评分会落在μi±σ 这个范围内
Statsmodels 在计量的简便性上是远远不及 Stata 等软件的，但它的优点在于可以与 Python 的其他的任务（如 NumPy、Pandas）有效结合，提高工作效率。下面让使用statsmodel 包来进行线性回归：

# 将一个字符串赋给变量my_foumula，该字符串中包含的是类似R语言语法的回归公式定义；波浪线（~）左侧的变量quality 是因变量，波浪线右侧的变量是自变量
my_formula = 'quality ~ alcohol + chlorides + citric_acid + density\
+ fixed_acidity + free_sulfur_dioxide + pH + residual_sugar + sulphates\
+ total_sulfur_dioxide + volatile_acidity'
# 使用公式和数据拟合一个普通最小二乘回归模型，并将结果赋给变量lm
lm = ols(my_formula, data=wine).fit()
## 或者，也可以使用广义线性模型（glm）语法进行线性回归
## lm = glm(my_formula, data=wine, family=sm.families.Gaussian()).fit()
#向屏幕上打印结果的摘要信息，包含了模型系数、系数的标准差和置信区间、修正R方、F统计量等模型详细信息
print(lm.summary())
print("\nQuantities you can extract from the result:\n%s" % dir(lm))
# 打印出一个列表，其中包含从模型对象lm 中提取出的所有数值信息，检查了这个列表之后
#希望提取出模型系数、系数的标准差、修正R方、F统计量和它的p值，以及模型拟合值。
# lm.params 以一个序列的形式返回模型系数，这样可以通过定位或名称提取出单个的系数
# 例如，要提取酒精含量的系数0.267，可以使用lm.params[1] 或lm.params['alcohol']
print("\nCoefficients:\n%s" % lm.params)
# lm.bse 以序列的形式返回模型系数的标准差
print("\nCoefficient Std Errors:\n%s" % lm.bse)
# lm.rsquared_adj 返回修正R方
print("\nAdj. R-squared:\n%.2f" % lm.rsquared_adj)
# lm.fvalue 和lm.f_pvalue 分别返回F统计量和它的p值
print("\nF-statistic: %.1f P-value: %.2f" % (lm.fvalue, lm.f_pvalue))
# lm.fittedvalues 返回拟合值
print("\nNumber of obs: %d Number of fitted values: %d" % (lm.nobs,\
len(lm.fittedvalues)))

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

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

逻辑斯蒂回归

在这个数据集中，因变量是一个二值变量，表示客户是否已经流失并不再是公司客户。线性回归不适合这种情况，因为它可能会生成小于0 或大于1 的预测结果，这在概率上是没有意义的。因为因变量是一个二值变量，所以需要将预测值限制在0 和1 之间。逻辑斯蒂回归可以满足这个要求。
Pr⁡(yi=1)=logit⁡−1(β0+β1xi1+β2xi2+⋯+βpxip)\operatorname{Pr}\left(y_{i}=1\right)=\operatorname{logit}^{-1}\left(\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{p} x_{i p}\right)Pr(yi=1)=logit−1(β0+β1xi1+β2xi2+⋯+βpxip)
对于i = 1, 2, …, n 个观测和p 个输入变量
等价于：
Pr⁡(yi=1)=pilogit⁡(pi)=(β0+β1xi1+β2xi2+⋯+βpxip)\begin{aligned} &\operatorname{Pr}\left(y_{i}=1\right)=p_{i}\\ &\operatorname{logit}\left(p_{i}\right)=\left(\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{p} x_{i p}\right) \end{aligned}Pr(yi=1)=pilogit(pi)=(β0+β1xi1+β2xi2+⋯+βpxip)

逻辑斯蒂回归通过使用逻辑函数（或称逻辑斯蒂函数）的反函数估计概率的方式来测量自变量和二值型因变量之间的关系。这个函数可以将连续值转换为0 和1 之间的值，这是个必要条件，因为预测值表示概率，而概率必须在0 和1 之间。这样，逻辑斯蒂回归预测的就是某种结果的概率，比如客户流失概率。
逻辑斯蒂回归通过一种能够实现极大似然估计的迭代算法来估计未知的β 参数值。
逻辑斯蒂回归的语法与线性回归有一点区别。对于逻辑斯蒂回归，需要分别设置因变量和自变量，而不是将它们写在一个公式中：

# 创建一个变量dependent_varible 并赋给它churn01 列中的一系列值
dependent_variable = churn['churn01']
# 设定了用作自变量的3 列，并将它们赋给变量independent_variables
independent_variables = churn[['account_length', 'custserv_calls',\
'total_charges']]
# 使用statsmodels 的add_constant 函数向输入变量中加入一列1
independent_variables_with_constant = sm.add_constant(independent_variables,\
prepend=True)
# 拟合逻辑斯蒂模型，并将拟合结果赋给变量logit_model
logit_model = sm.Logit(dependent_variable, independent_variables_with_constant)\
.fit()# 向屏幕上打印模型的摘要信息。这个摘要信息非常重要，因为其中包含了模型系数、系数标准差和置信区间、伪R方等模型详细信息
print(logit_model.summary())
# 打印出一个列表，其中包含从模型对象logit_model 中提取出的所有数值信息。检查了这个列表之后，提取出模型系数和它们的标准差。
print("\nQuantities you can extract from the result:\n%s" % dir(logit_model))
# logit_model.params 以一个序列的形式返回模型系数，这样便可以通过定位或名称提取出单个模型系数。同样，logit_model.bse 以一个序列的形式返回系数标准差。
print("\nCoefficients:\n%s" % logit_model.params)
print("\nCoefficient Std Errors:\n%s" % logit_model.bse)