Question

下面来分析客户流失数据集。首先，将数据读入一个数据框，然后格式化列标题，为数据框 churn 创建一个新的数值型二值变量，并检查数据框中前面几行数据。为了完成这些操作，需要创建一个新脚本 customer_churn.py，并输入以下各行代码：

#!/usr/bin/env python3
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
churn = pd.read_csv('churn.csv', sep=',', header=0)
churn.columns = [heading.lower() for heading in \
churn.columns.str.replace(' ', '_').str.replace("\'", "").str.strip('?')]
churn['churn01'] = np.where(churn['churn'] == 'True.', 1., 0.)
print(churn.head())

import 语句之后的第一行代码将数据读入数据框 churn。下一行代码两次使用 replace 函数将列标题中的空格替换成下划线，并删除嵌入的单引号。请注意第二个 replace 函数，其中第一个参数的两个双引号中间是一个反斜线加一个单引号，逗号后面的第二个参数就是一对双引号。这行代码还使用 strip 函数除去了列标题 Churn? 末尾的问号。最后，这行代码使用列表生成式将所有列标题转换为小写。

下一行代码创建一个新列 churn01，并使用 numpy 的 where 函数根据 churn 这一列中的值用 1 或 0 来填充它。churn 这一列中的值不是 True 就是 False，所以如果 churn 中的值是 True，那么 churn01 中的值就是 1，如果 churn 中的值是 False，那么 churn01 中的值就是 0。

最后一行代码使用 head 函数显示标题行和前 5 个数据行，这样可以检查一下数据加载以及列标题的格式化是否正确。

在将数据加载到数据框中之后，可以通过计算流失客户和未流失客户的描述性统计量，来看看这两组客户有什么区别。下面的代码按照 churn 这一列中的值将数据分成了两组：已流失的客户和未流失的客户。然后为每个分组中的一些特定的列计算 3 个统计量：总数、均值和标准差：

# 为分组数据计算描述性统计量
print(churn.groupby(['churn'])[['day_charge', 'eve_charge', 'night_charge',\
'intl_charge', 'account_length', 'custserv_calls']].agg(['count', 'mean',\
'std']))

下面的代码演示了如何为不同的变量计算不同类型的多个统计量。这行代码为 4 个变量计算均值和标准差，为 2 个变量计算总数、最小值和最大值。我们还是按照 churn 列分组，所以代码分别为流失客户和未流失客户计算这些统计量：

# 为不同的变量计算不同的统计量
print(churn.groupby(['churn']).agg({'day_charge' : ['mean', 'std'],
                   'eve_charge' : ['mean', 'std'],
                   'night_charge' : ['mean', 'std'],
                   'intl_charge' : ['mean', 'std'],
                   'account_length' : ['count', 'min', 'max'],
                   'custserv_calls' : ['count', 'min', 'max']}))

下一段代码对客户服务通话次数这部分数据进行了摘要分析，先按照一个新变量 total_charges 中的值使用等宽分箱法将数据分成 5 个组，然后为每个分组计算 5 个统计量：总数、最小值、均值、最大值和标准差。为了完成这些操作，第一行代码创建一个新变量 total_charges，表示白天、傍晚、夜间和国际通话费用的总和。下一行代码使用 cut 函数按照等宽分箱法将 total_charges 分成 5 组。然后定义一个函数 get_stats，为每个分组返回一个统计量字典。下一行代码按照 5 个 total_charges 分组将客户服务通话次数也分成同样的 5 组。最后，在分组数据上应用 get_stats 函数，为 5 个分组计算统计量：

# 创建total_charges
# 将其分为5组，并为每一组计算统计量
churn['total_charges'] = churn['day_charge'] +churn['eve_charge'] +\
churn['night_charge'] + churn['intl_charge']
factor_cut = pd.cut(churn.total_charges, 5, precision=2)
def get_stats(group):
    return {'min' : group.min(), 'max' : group.max(),
              'count' : group.count(), 'mean' : group.mean(),
              'std' : group.std()}
grouped = churn.custserv_calls.groupby(factor_cut)
print(grouped.apply(get_stats).unstack())

和前一段代码相似，下一段代码也是用 5 个统计量对客户服务通话数据进行了摘要分析。但是，这段代码使用 qcut 函数通过等深分箱法（按照分位数进行划分）将 account_length 分成了 4 组，而不是使用等宽分箱法对数据进行的分组：

# 将account_length按照分位数进行分组
# 并为每个分位数分组计算统计量
factor_qcut = pd.qcut(churn.account_length, [0., 0.25, 0.5, 0.75, 1.])
grouped = churn.custserv_calls.groupby(factor_qcut)
print(grouped.apply(get_stats).unstack())
通过分位数对account_length 进行划分，可以保证

通过分位数对 account_length 进行划分，可以保证每个分组中包含数目大致相同的观测。前一段代码中通过等宽分箱法得到的每个分组中包含的观测数目是不一样的。qcut 函数使用一个整数或一个分位数数组来设定分位数的数量，所以你可以使用整数 4 来代替 [0., 0.25, 0.5, 0.75, 1.] 设定 4 等分，或使用 10 来设定 10 等分。

下面一段代码演示了如何使用 pandas 的 get_dummies 函数来创建二值指标变量，并将这些变量添加进数据框。前两行代码为 intl_plan 列和 vmail_plan 列创建二值指标变量，并使用原来的变量名作为新列的前缀。下一行使用 join 命令将 churn 列和新的二值指标列合并，并将结果赋给一个新的文本框 churn_with_dummies。这个新文本框有 5 列：churn、intl_plan_no、intl_plan_yes、vmail_plan_no 和 vmail_plan_yes。

# 为intl_plan和vmail_plan创建二值（虚拟）指标变量
# 并将它们与新数据框中的churn列连接起来
intl_dummies = pd.get_dummies(churn['intl_plan'], prefix='intl_plan')
vmail_dummies = pd.get_dummies(churn['vmail_plan'], prefix='vmail_plan')
churn_with_dummies = churn[['churn']].join([intl_dummies, vmail_dummies])
print(churn_with_dummies.head())

下面这段代码演示了如何将一列按照四分位数进行划分，为每个四分位数创建二值指标变量，并将新列添加到原来的数据框中。qcut 函数将 total_charges 列按照四分位数进行划分，并使用 qcut_names 中的名称对每个四分位数进行标记。get_dummies 函数为四分位数创建 4 个二值指标变量，并使用 total_charges 作为新列的前缀。最终结果为 4 个新的虚拟变量：total_charges_1st_quartile、total_charges_2nd_quartile、total_charges_3rd_quartile、total_charges_4th_quartile。join 函数将这 4 个变量追加到数据框 churn 中。

# 将total_charges按照分位数分组，为每个分位数分组创建一个二值指标变量
# 并将它们加入到churn数据框中
qcut_names = ['1st_quartile', '2nd_quartile', '3rd_quartile', '4th_quartile']
total_charges_quartiles = pd.qcut(churn.total_charges, 4, labels=qcut_names)
dummies = pd.get_dummies(total_charges_quartiles, prefix='total_charges')
churn_with_dummies = churn.join(dummies)
print(churn_with_dummies.head())

最后一段代码创建了 3 个数据透视表。第一行代码在对 total_charges 列按照流失情况和客户服务通话次数进行透视转换（或行列分组）之后，计算每组的均值。结果是一长列数值，表示每个流失情况和客户服务通话次数组合的平均总费用。第二行代码表示对结果重新格式化，使用流失情况作为行，客户服务通话次数作为列。最后，第三行代码使用客户服务通话次数作为行，流失情况作为列，演示了指定要计算的统计量、处理缺失值和是否显示边际值的方法：

# 创建透视表
print(churn.pivot_table(['total_charges'], index=['churn', 'custserv_calls']))
print(churn.pivot_table(['total_charges'], index=['churn'],\
columns=['custserv_calls']))
print(churn.pivot_table(['total_charges'], index=['custserv_calls'],\
columns=['churn'], aggfunc='mean', fill_value='NaN', margins=True ))

7.3.1　逻辑斯蒂回归

在这个数据集中，因变量是一个二值变量，表示客户是否已经流失并不再是公司客户。线性回归不适合这种情况，因为它可能会生成小于 0 或大于 1 的预测结果，这在概率上是没有意义的。因为因变量是一个二值变量，所以需要将预测值限制在 0 和 1 之间。逻辑斯蒂回归可以满足这个要求。

逻辑斯蒂回归模型如下所示。

·

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

等价于：

逻辑斯蒂回归通过使用逻辑函数（或称逻辑斯蒂函数）的反函数估计概率的方式来测量自变量和二值型因变量之间的关系。这个函数可以将连续值转换为 0 和 1 之间的值，这是个必要条件，因为预测值表示概率，而概率必须在 0 和 1 之间。这样，逻辑斯蒂回归预测的就是某种结果的概率，比如客户流失概率。

逻辑斯蒂回归通过一种能够实现极大似然估计的迭代算法来估计未知的 β 参数值。

逻辑斯蒂回归的语法与线性回归有一点区别。对于逻辑斯蒂回归，需要分别设置因变量和自变量，而不是将它们写在一个公式中：

dependent_variable = churn['churn01']
independent_variables = churn[['account_length', 'custserv_calls',\
'total_charges']]
independent_variables_with_constant = sm.add_constant(independent_variables,\
prepend=True)

logit_model = sm.Logit(dependent_variable, independent_variables_with_constant)\
.fit()

print(logit_model.summary())
print("\nQuantities you can extract from the result:\n%s" % dir(logit_model))
print("\nCoefficients:\n%s" % logit_model.params)
print("\nCoefficient Std Errors:\n%s" % logit_model.bse)

第一行代码创建一个变量 dependent_varible 并赋给它 churn01 列中的一系列值。

同样，第二行代码设定了用作自变量的 3 列，并将它们赋给变量 independent_variables。

然后，我们使用 statsmodels 的 add_constant 函数向输入变量中加入一列 1。

下一行代码拟合逻辑斯蒂模型，并将拟合结果赋给变量 logit_model。

最后 4 行代码向屏幕上打印出模型中具体的结果。第一行代码向屏幕上打印模型的摘要信息。这个摘要信息非常重要，因为其中包含了模型系数、系数标准差和置信区间、伪 R 方等模型详细信息。

下一行代码打印出一个列表，其中包含从模型对象 logit_model 中提取出的所有数值信息。检查了这个列表之后，提取出模型系数和它们的标准差。

接下来的 2 行代码提取出了这些值。logit_model.params 以一个序列的形式返回模型系数，这样便可以通过定位或名称提取出单个模型系数。同样，logit_model.bse 以一个序列的形式返回系数标准差。输出结果如图 7-7 所示。

图 7-7：使用 3 个变量对客户流失进行多元逻辑斯蒂回归

7.3.2　系数解释

对逻辑斯蒂回归系数的解释不像线性回归那么直观，因为逻辑斯蒂函数的反函数是条曲线，这说明自变量一个单位的变化所造成的因变量的变化不是一个常数。

因为逻辑斯蒂函数的反函数是一条曲线，所以必须选择使用哪个函数值来评价自变量对成功概率的影响。和线性回归一样，截距系数的意义是当所有自变量为 0 时成功的概率。有些时候 0 是没有意义的，所以另外一种方式是当自变量都取均值时，看看函数的值有何意义：

def inverse_logit(model_value):
    from math import exp
    return (1.0 / (1.0 + exp(-model_value)))

at_means = float(logit_model.params[0]) + \
    float(logit_model.params[1])*float(churn['account_length'].mean()) + \
    float(logit_model.params[2])*float(churn['custserv_calls'].mean()) + \
    float(logit_model.params[3])*float(churn['total_charges'].mean())

print("Probability of churn at mean values: %.2f" % inverse_logit(at_means))

第一段代码定义了一个函数 inverse_logit，将线性模型得出的连续预测值转换为 0 和 1 之间的概率值。

第二段代码计算出当所有自变量取均值时观测的预测值。4 个 logit_model.params[…] 值是模型系数，3 个 churn['…'].mean() 分别是账户时间、客户服务通话次数和总费用的均值。给定了模型系数和各变量的均值，方程变为 -7.2+(0.001*101.1)+(0.444*1.6)+(0.073*59.5)，所以变量 at_means 的值为 -2.068。

最后一行代码在屏幕上打印出逻辑斯蒂函数的反函数在 at_means 处的值，保留两位小数。-2.068 处的反函数值为 0.112，所以当账户时间、客户服务通话次数和总费用均取均值时，客户流失的概率就是 11.2%。

同样，要计算某个自变量一个单位的变化造成的因变量的变化，可以通过计算当某个自变量从均值发生一个单位的变化时，成功概率发生了多大的变化。

例如，下面来计算一下当客户服务通话次数在均值的基础上发生一个单位的变化时，对客户流失概率造成的影响：

cust_serv_mean = float(logit_model.params[0]) + \
    float(logit_model.params[1])*float(churn['account_length'].mean()) + \
    float(logit_model.params[2])*float(churn['custserv_calls'].mean()) + \
    float(logit_model.params[3])*float(churn['total_charges'].mean())

cust_serv_mean_minus_one = float(logit_model.params[0]) + \
    float(logit_model.params[1])*float(churn['account_length'].mean()) + \
    float(logit_model.params[2])*float(churn['custserv_calls'].mean()-1.0) + \
    float(logit_model.params[3])*float(churn['total_charges'].mean())

print("Probability of churn when account length changes by 1: %.2f" % \
(inverse_logit(cust_serv_mean) - inverse_logit(cust_serv_mean_minus_one)))

第一段代码与计算 at_means 的代码相同。第二段代码也与其基本相同，区别在于将客户服务通话次数的均值减去了 1。

最后一行代码打印出了两个逻辑斯蒂函数反函数值的差，其中一个是所有自变量为均值时的反函数值，另一个是两个自变量为均值，客户服务通话次数为均值减 1 时的反函数值。

在这个示例中，cust_serv_mean 的值与 at_means 相同，都是 -2.068。cust_serv_mean_minus_one 的值是 -2.512。-2.068 的反函数值减去 -2.512 的反函数值的结果是 0.0372，所以在均值附近减少一次客户服务通话就对应着客户流失概率提高 3.7 个百分点。

7.3.3　预测

与前面的葡萄酒质量预测一样，你也可以使用这个拟合模型来对“新”观测进行预测：

# 在churn数据集中
# 使用前10个观测创建10个“新”观测
new_observations = churn.ix[churn.index.isin(range(10)),\
independent_variables.columns]
new_observations_with_constant = sm.add_constant(new_observations, prepend=True)

# 基于新观测的账户特性
# 预测客户流失可能性
y_predicted = logit_model.predict(new_observations_with_constant)

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

同样，变量 y_predicted 中包含着 10 个预测值。为了使输出更简单易懂，可以将预测值保留两位小数。现在我们得到了可用的预测值，如果使用的是真正的新观测，那么就可以使用这些预测值来评价模型了。