作者：Ben，著有《Python数据科学：技术详解与商业实践》、《用商业案例学R语言数据挖掘》、《胸有成竹-数据分析的SAS EG进阶》作者。2005年进入数据科学领域，先后在亚信、德勤、百度等企业从事电信、金融行业数据挖掘工作。专注于消费信贷和财富管理方面的数据治理、客户智能与风险智能领域。

配套视频教程：Python数据科学-技术详解与商业实践(第二期)课程链接: https://edu.hellobi.com/course/280

电信公司希望针对客户的信息预测其流失可能性，数据存放在“telecom_churn.csv”中。

1、分析思路：

在对客户流失与否的影响因素进行模型研究之前，首先对各解释变量与被解释变量进行两变量独立性分析，以初步判断影响流失的因素，进而建立客户

流失预测模型

主要变量说明如下：

#subscriberID="个人客户的ID"

#churn="是否流失：1=流失";

#Age="年龄"

#incomeCode="用户居住区域平均收入的代码"

#duration="在网时长"

#peakMinAv="统计期间内最高单月通话时长"

#peakMinDiff="统计期间结束月份与开始月份相比通话时长增加数量"

#posTrend="该用户通话时长是否呈现出上升态势：是=1"

#negTrend="该用户通话时长是否呈现出下降态势：是=1"

#nrProm="电话公司营销的数量"

#prom="最近一个月是否被营销过：是=1"

#curPlan="统计时间开始时套餐类型：1=最高通过200分钟；2=300分钟；3=350分钟；4=500分钟"

#avPlan="统计期间内平均套餐类型"

#planChange="统计结束时和开始时套餐的变化：正值代表套餐档次提升，负值代表下降，0代表不变"

#posPlanChange="统计期间是否提高套餐：1=是"

#negPlanChange="统计期间是否降低套餐：1=是"

#call_10086="拨打10086的次数"

2、作业安排：

2.1 基础知识：

1)作了一次营销活动，营销了1000人 。事后统计结果，120人购买，其余人没有购买。请分别用矩估计法、极大似然估计法计算这个随机事件分布的参数(提示：该随机事件服从伯努利分布)

2)推导线性回归参数估计的最小二乘、矩估计、极大似然估计，推导逻辑回归的极大似然估计公式。线性回归和逻辑回归的极大似然法哪个可以得到

显性的公式解，哪个需要使用迭代法求解？解释极大似然法求解过程中用到的牛顿迭代法、随机梯度法的做法。

3)解释统计学习算法中超参的概念，请问目前统计方法中学习的线性回归、逻辑回归中涉及超参了吗？岭回归和Laso算法中的超参分别是什么？超参的作用是什么？统计学习算法中如何确定最优超参的取值？

4)比较统计分析法和统计学习(即机器学习)得到最优模型的思路。

5)二分类模型中(比如逻辑回归)的评估模型优劣的决策类和排序类评估指标分别包括哪些指标?

*2.1.1答案请参见附件中的PDF

2.2 案例解答步骤如下：

1)两变量分析：检验该用户通话时长是否呈现出上升态势(posTrend)对流失(churn) 是否有预测价值

2)首先将原始数据拆分为训练和测试数据集，使用训练数据集建立在网时长对流失的逻辑回归，使用测试数据集制作混淆矩阵(阈值为0.5)，提供准

确性、召回率指标，提供ROC曲线和AUC。

3)使用向前逐步法从其它备选变量中选择变量，构建基于AIC的最优模型，绘制ROC曲线，同时检验模型的膨胀系数。

4)使用岭回归和Laso算法重建第三步中的模型，使用交叉验证法确定惩罚参数(C值)。并比较步骤四中Laso算法得到的模型和第三步得到的模型的差异

。

2.2.1参考代码如下:

# In[1]:

import os

import numpy as np

from scipy import stats

import pandas as pd

import statsmodels.api as sm

import statsmodels.formula.api as smf

import matplotlib.pyplot as plt

#pd.set_option('display.max_columns', None)

os.chdir(r"D:\Python_Training\script_Python\8logistic\HW")

# 导入数据和数据清洗

# In[5]:

churn = pd.read_csv(r'telecom_churn.csv', skipinitialspace=True)

churn.head()

#1两变量分析：检验该用户通话时长是否呈现出上升态势(posTrend)对流失(churn) 是否有预测价值

# ## 分类变量的相关关系

#

# 交叉表

# In[6]:

cross_table = pd.crosstab(churn.posTrend,

churn.churn, margins=True)

cross_table

# 列联表

# In[7]:

def percConvert(ser):

return ser/float(ser[-1])

cross_table.apply(percConvert, axis=1)

# In[8]:

print('''chisq = %6.4f

p-value = %6.4f

dof = %i

expected_freq = %s''' %stats.chi2_contingency(cross_table.iloc[:2, :2]))

#2首先将原始数据拆分为训练和测试数据集，使用训练数据集建立在网时长对流失的逻辑回归，使用测试数据集制作混淆矩阵

#(阈值为0.5)，提供准确性、召回率指标，提供ROC曲线和AUC。

# ## 逻辑回归

# In[9]:

churn.plot(x='duration', y='churn', kind='scatter')

# •随机抽样，建立训练集与测试集

# In[10]:

train = churn.sample(frac=0.7, random_state=1234).copy()

test = churn[~ churn.index.isin(train.index)].copy()

print(' 训练集样本量: %i \n 测试集样本量: %i' %(len(train), len(test)))

# In[11]:

lg = smf.glm('churn ~ duration', data=train,

family=sm.families.Binomial(sm.families.links.logit)).fit()

lg.summary()

# 预测

# In[19]:

train['proba'] = lg.predict(train)

test['proba'] = lg.predict(test)

test['proba'].head()

# In[12]:

# ## 模型评估

#

# 设定阈值

# In[20]:

test['prediction'] = (test['proba'] > 0.5).astype('int')

# 混淆矩阵

# In[22]:

pd.crosstab(test.churn, test.prediction, margins=True)

# - 计算准确率

# In[23]:

acc = sum(test['prediction'] == test['churn']) /np.float(len(test))

print('The accurancy is %.2f' %acc)

# In[25]:

for i in np.arange(0.1, 0.9, 0.1):

prediction = (test['proba'] > i).astype('int')

confusion_matrix = pd.crosstab(prediction,test.churn,

margins = True)

precision = confusion_matrix.ix[0, 0] /confusion_matrix.ix['All', 0]

recall = confusion_matrix.ix[0, 0] / confusion_matrix.ix[0, 'All']

Specificity = confusion_matrix.ix[1, 1] /confusion_matrix.ix[1,'All']

f1_score = 2 * (precision * recall) / (precision + recall)

print('threshold: %s, precision: %.2f, recall:%.2f ,Specificity:%.2f , f1_score:%.2f'%(i, precision, recall, Specificity,f1_score))

# - 绘制ROC曲线

# In[27]:

import sklearn.metrics as metrics

fpr_test, tpr_test, th_test = metrics.roc_curve(test.churn, test.proba)

fpr_train, tpr_train, th_train = metrics.roc_curve(train.churn, train.proba)

plt.figure(figsize=[3, 3])

plt.plot(fpr_test, tpr_test, 'b--')

plt.plot(fpr_train, tpr_train, '')

plt.title('ROC curve')

plt.show()

# In[28]:

print('AUC = %.4f' %metrics.auc(fpr_test, tpr_test))

#3使用向前逐步法从其它备选变量中选择变量，构建基于AIC的最优模型，绘制ROC曲线，同时检验模型的膨胀系数。

# In[14]:

#- 多元逻辑回归

# 向前法

def forward_select(data, response):

remaining = set(data.columns)

remaining.remove(response)

selected = []

current_score, best_new_score = float('inf'), float('inf')

while remaining:

aic_with_candidates=[]

for candidate in remaining:

formula = "{} ~ {}".format(

response,' + '.join(selected + [candidate]))

aic = smf.glm(

formula=formula, data=data,

family=sm.families.Binomial(sm.families.links.logit)

).fit().aic

aic_with_candidates.append((aic, candidate))

aic_with_candidates.sort(reverse=True)

best_new_score, best_candidate=aic_with_candidates.pop()

if current_score > best_new_score:

remaining.remove(best_candidate)

selected.append(best_candidate)

current_score = best_new_score

print ('aic is {},continuing!'.format(current_score))

else:

print ('forward selection over!')

break

formula = "{} ~ {} ".format(response,' + '.join(selected))

print('final formula is {}'.format(formula))

model = smf.glm(

formula=formula, data=data,

family=sm.families.Binomial(sm.families.links.logit)

).fit()

return(model)

# In[16]:

candidates = ['churn','duration','AGE','edu_class','posTrend','negTrend','nrProm','prom','curPlan','avgplan','planChange','incomeCode','feton','peakMinAv','peakMinDiff','call_10086']

data_for_select = train[candidates]

lg_m1 = forward_select(data=data_for_select, response='churn')

lg_m1.summary()

# Seemingly wrong when using 'statsmmodels.stats.outliers_influence.variance_inflation_factor'

# In[17]:

def vif(df, col_i):

from statsmodels.formula.api import ols

cols = list(df.columns)

cols.remove(col_i)

cols_noti = cols

formula = col_i + '~' + '+'.join(cols_noti)

r2 = ols(formula, df).fit().rsquared

return 1. / (1. - r2)

# In[18]:

exog = train[candidates].drop(['churn'], axis=1)

for i in exog.columns:

print(i, '\t', vif(df=exog, col_i=i))

#posTrend,negTrend;curPlan,avgplan有明显的共线性问题,剔除其中两个后重新建模.

#%%

#4)使用岭回归和Laso算法重建第三步中的模型，使用交叉验证法确定惩罚参数(C值)。并比较步骤四中Laso算法得到的模型和第三步得到的模型的差异

from sklearn.linear_model import LogisticRegression

from sklearn.preprocessing import StandardScaler

candidates = ['duration','AGE','edu_class','posTrend','negTrend','nrProm','prom','curPlan','avgplan','planChange','incomeCode','feton','peakMinAv','peakMinDiff','call_10086']

#data_for_select = churn[candidates]

scaler = StandardScaler() # 标准化

X = scaler.fit_transform(churn[candidates])

y = churn['churn']

#%%

from sklearn import linear_model

from sklearn.svm import l1_min_c

cs = l1_min_c(X, y, loss='log') * np.logspace(0, 4)

print("Computing regularization path ...")

#start = datetime.now()

clf = linear_model.LogisticRegression(C=1.0, penalty='l1', tol=1e-6)

coefs_ = []

for c in cs:

clf.set_params(C=c)

clf.fit(X, y)

coefs_.append(clf.coef_.ravel().copy())

#print("This took ", datetime.now() - start)

coefs_ = np.array(coefs_)

plt.plot(np.log10(cs), coefs_)

ymin, ymax = plt.ylim()

plt.xlabel('log(C)')

plt.ylabel('Coefficients')

plt.title('Logistic Regression Path')

plt.axis('tight')

plt.show()

#%%

cs = l1_min_c(X, y, loss='log') * np.logspace(0, 4)

import matplotlib.pyplot as plt #可视化模块

from sklearn.cross_validation import cross_val_score # K折交叉验证模块

k_scores = []

clf = linear_model.LogisticRegression(penalty='l1')

#藉由迭代的方式来计算不同参数对模型的影响，并返回交叉验证后的平均准确率

for c in cs:

clf.set_params(C=c)

scores = cross_val_score(clf, X, y, cv=10, scoring='roc_auc') #http://scikit-learn.org/stable/modules/model_evaluation.html

k_scores.append([c,scores.mean(),scores.std()])

#%%

#可视化数据

#%%

data=pd.DataFrame(k_scores)#将字典转换成为数据框

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.plot(np.log10(data[0]), data[1],'b')

ax1.set_ylabel('Mean ROC(Blue)')

ax1.set_xlabel('log10(cs)')

ax2 = ax1.twinx()

ax2.plot(np.log10(data[0]), data[2],'r')

ax2.set_ylabel('Std ROC Index(Red)')

#得到合理的C为 np.exp(-1.9)

#%%

#重新实现Laso算法

from sklearn.linear_model import LogisticRegression

from sklearn.preprocessing import StandardScaler

candidates = ['duration','AGE','edu_class','posTrend','negTrend','nrProm','prom','curPlan','avgplan','planChange','incomeCode','feton','peakMinAv','peakMinDiff','call_10086']

#data_for_select = churn[candidates]

scaler = StandardScaler() # 标准化

X = scaler.fit_transform(churn[candidates])

y = churn['churn']

#%%

from sklearn import linear_model

clf = linear_model.LogisticRegression(C=np.exp(-1.9), penalty='l1')

clf.fit(X, y)

clf.coef_

#%%