本實戰案共分為五個部分：商業背景、指標設計、部署環境、數據準備、回歸分析，其中回歸分析包括：模型構建、模型診斷、模型結果、模型應用。

??數據STUDIO投稿?·?作者｜玄武

1.商業背景

眾所周知，移動通信市場已經日趨飽和，增加規模已經變得異常艱難，通信運營商互挖墻角已經成為家常便事。很多消費者，今天還是中國移動的客戶，明天只要中國電信給點好處，就變成中國電信的客戶，后天一看中國聯通推出打折促銷活動，又變成中國聯通的客戶，再過幾天，中國移動稍微關懷一下，又重新回到中國移動的懷抱。在這樣一個周而復始的拉鋸戰中，通信運營商耗盡了有限的營銷資源，客戶也沒有得到實質性的好處，因為更換運營商其實也是一種消耗。此時，增強客戶的忠誠度，提升公司的盈利能力，對通信運營商來說，就變得非常重要。

對于兵家必爭的高校大學生市場，某通信運營商針對大學生群體，推出了校園網計劃來提升客戶忠誠度。按照校園網運營規則，如果一名大學生希望加入校園網，他首先必須是該運營商的用戶，否則無法參加該計劃，此外，還得由現有校園網用戶進行邀請。作為回報，校園網內所有通話資費非常便宜，而且數據流量優惠也非常巨大，但與網外朋友通信資費照舊。所以，為了降低自身資費，現有校園網成員都有很大動力邀請朋友加入校園網。這樣的話，大量的日常溝通將發生在校園網內，這不僅降低了通信資費，還享受了更好的溝通服務，同時，已經加入校園網的成員則發現很難離開，因為大部分朋友以及主要通信社交網絡都在校園網中，個人一旦離開，如果還想保持跟過去一樣的溝通強度，成本將非常昂貴。

那么，通信運營商的付出與回報又如何呢？為了深度“套牢”在校大學生客戶，運營商有著重要的付出，也就是降低資費。此外，為了迅速擴張校園網，鼓勵大家推薦新客戶，運營商對推薦者有一定的獎勵，比如獎勵話費或者流量甚至現金。

那么，通信運營商付出這么多，希望得到的回報是什么呢？

第一個回報，高忠誠度、低離網率。既然用戶把好友都拉進了校園網，那么用戶的主要通信社交網絡都被校園網覆蓋，這能否增加忠誠度，降低離網率，進而間接降低客戶的獲取和維護成本？

第二個回報，總利潤不降反升。雖然下調了資費，但是，資費的下調很可能會刺激消費量的上升，從而使得最終總利潤不降反升。

但是，運營一段時間后發現，好消息是離網率確實下降不少、總利潤也有所上升，壞消息是總利潤上升低于預期。因為有些校園網用戶邀請了很多低端客戶進來，這些低端客戶的總消費量并沒有因為入網而有任何上升，相反，由于資費的下降，他們給公司貢獻的收入卻大幅下降。但是，也有些校園網用戶邀請了很多優質客戶進來，相比入網前，這些優質客戶的溝通更加密切，因此，盡管單位時長的資費水平下降很多，但是他們對企業的總利潤貢獻卻上升不少。

這說明，不同的客戶作為推薦者，能帶來的被推薦者的價值是不一樣的，也就是說，并非每個人都能推薦有價值的客戶，甚至有些推薦者帶來的客戶的貢獻是負的。因為，一名被推薦進來的客戶，他對校園網貢獻的大小，除依賴自身特征外，還極大地依賴于推薦者，如果推薦者是話癆，被推薦者很可能就是那個被話癆的對象，如果推薦者是個游戲高手，被推薦者很可能就是他游戲中的豬隊友。因此，有必要研究，帶來低價值客戶的推薦者與帶來高價值客戶的推薦者之間有沒有系統性差異？如果能夠掌握此規律，就可把有限的獎勵資源，有針對性地投放到那些能為企業帶來高價值客戶的推薦者身上。

因此，需要詳細研究：什么樣的推薦者能夠帶來高價值客戶，什么樣的推薦者帶來的客戶是低價值客戶？

2.指標設計

這個時候，我們就需要設計一個指標來衡量推薦者的價值，并且這個指標必須對業務有好的指導意義。那么，推薦者的價值應當通過什么指標來評估呢？這個指標就是我們研究的因變量。我相信，大家可以找出不少指標來評估推薦者的價值，比如，推薦者所推薦的所有客戶，在加入校園網后，其絕對和相對收入的增長或者絕對和相對利潤的增長，在這里，我們使用“某推薦者當月推薦的所有客戶加入校園網后次月的利潤環比增長率”為評估指標，也就是為我們研究的因變量。

在確定因變量之后，我們則需要考慮有哪些因素會影響著推薦者的價值，也就是需要尋找自變量。在實際工作中，我們有大量的有用指標，能夠詳細地刻畫推薦者的方方面面。比如，可以考慮消費者的消費行為，主要包括該用戶在各項通信及增值業務上的花費。再比如，還可以考慮消費者的通話特征，包括該用戶的通話時長、頻率、時間等，甚至還可以將通話時長拆分成主叫、被叫、本地、長途、漫游等?？偠灾?#xff0c;實際工作中可考慮的指標很多，它們都有助于更好地描述推薦者，它們都可以成為自變量。在這里，純粹為了簡單起見，只考慮這么三個自變量。

第一個自變量，月通話總量

月通話總量，指該推薦者進行推薦的當月的通話總時長，以百分鐘計。很顯然，這是個重要的自變量，它直接刻畫了用戶的活躍程度。由于校園網提供非常優惠的通話資費，因此對那些通話總量高的用戶很有吸引力，假設高通話量用戶的好友也可能是高通話量用戶，那么，具有高通話量特征的推薦者就更有可能帶來優質客戶。

第二個自變量，大網占比

大網指的是該運營商的通信網絡，大網占比，就是在該推薦者推薦當月的通話總時長中，發生在該運營商網內的通話總時長占比。這個占比越高，說明該推薦者的通話，越多地發生在該運營商生態體系之內，因此，他具有更好的基礎為校園網計劃拉新客戶。

第三個自變量，小網占比

小網指的是校園網，小網占比，就是在該推薦者推薦當月發生在該運營商網內的通話總時長中，發生在校園網內的通話總時長占比。

大網占比，衡量了推薦者全部通話中有多少發生在該運營商網內，小網占比衡量的則是發生在該運營商網內的通話時長中，又有多少發生在校園網。由于加入校園網的前提條件是，消費者為該運營商的客戶，所以，一般來說，每個推薦者能夠發展多少個校園網用戶是有上限的，這個上限就是該推薦者所有發生在該運營商網絡內部的社交關系。如果一名推薦者的小網占比很高，就說明該推薦者可被推薦的社會關系網絡中的絕大部分已經加入了校園網，因此，該推薦者沒有更多的被推薦對象，所以推測他能為企業帶來的價值也許不高。

3.部署環境

3.1 導入依賴庫

包括數據處理模塊、機器學習模塊、統計分析模塊等等。

import?numpy?as?np import?pandas?as?pd import?statsmodels.api?as?sm from?scipy?import?stats from?itertools?import?combinations from?sklearn.cluster?import?KMeans??#?K均值聚類 from?pandas.plotting?import?lag_plot,autocorrelation_plot??#?滯后殘差圖、自相關系數折線圖 from?statsmodels.stats.diagnostic?import?acorr_breusch_godfrey?#序列自相關拉格朗日乘數檢驗(BG檢驗、LM檢驗) from?statsmodels.graphics.tsaplots?import?plot_acf,?plot_pacf?#自相關圖、偏自相關圖 from?sklearn.preprocessing?import?StandardScaler??#?數據標準化 from?IPython.core.display_functions?import?display

3.2 設置可視化

主要針對數據可視化部分，進行全局樣式變量進行配置。

import?seaborn?as?sns sns.set_style('white')??#?當出現pycharm風格錯亂時，可使用此語句修正 from?matplotlib?import?pyplot?as?plt plt.rcParams['xtick.bottom']=True??#?設置橫軸刻度線可見 plt.rcParams['ytick.left']=True??#?設置縱軸刻度線可見 plt.rcParams['figure.autolayout']=True??#自動tight布局，解決多個plots出現重疊問題 plt.rcParams['axes.unicode_minus']=False??#?解決負號顯示問題 plt.rcParams['font.sans-serif']=['Microsoft?YaHei']??#?解決中文顯示問題 plt.rcParams['figure.dpi']=80??#?設備圖形分辨率 %matplotlib?inline

3.3定義函數

繪制分組箱線圖。

def?grouping_boxplot(grouping_variable,grouping_variable_label,box_variable,box_variable_label,data): #?繪制分組箱線圖boxplot_data=[data[data[grouping_variable]==value][box_variable].dropna()?for?value?in?data[grouping_variable].value_counts().index]??#?構造用于制作箱形圖的數據plt.boxplot(boxplot_data,??#?數據結構為列表套列表，一個列表繪制一個箱線圖labels=data[grouping_variable].value_counts().index,??#?每個箱線圖在橫軸的標簽patch_artist=False,??#?設置是否往箱體填充顏色，默認不填充(False)showmeans=True,??#?在箱子內標出均值所在位置，箱體內部的橫線為中位數widths=0.5,??#?設置箱子的寬度，默認為0.5showfliers=False,??#?設置是否顯示異常值，默認顯示異常值(True)medianprops={'linestyle':?'-',?'color':?'black'}??#?設置中位數線的類型和顏色)plt.xlabel(grouping_variable_label)plt.ylabel(box_variable_label)

步進法(stepwise)選擇進入模型的自變量。

上下滑動查看更多源碼def?stepwise_select_variable(x_var,y_var,pvalue_in=0.05,pvalue_out=0.1): #?步進法(stepwise)選擇進入模型的自變量： #?如果p值小于pvalue_in，則往模型中放入此變量； #?如果p值大于pvalue_out，則從模型中移出此變量variables_in=[]??#?用于存儲模型中的變量while?True:??#?直接設置為True，則讓while永遠循環，直到遇到break才跳出循環flag=False#?---------#?挑選p值最小且小于pvalue_in的自變量variables_out=list(set(x_var.columns)-set(variables_in))??#?用于存儲未進入模型的變量，#?其值為矩陣x的所有變量列減去模型內的變量列（set是一個集合，無序且內容不能重復）pvalues=pd.Series(index=variables_out,dtype='float64')??#?創建一個series（索引為未進入模型的變量的名稱），#?用于存儲自變量回歸系數(不含截距)的p值for?new_column?in?variables_out:??#?遍歷未進入模型的變量x_var_cons=sm.add_constant(x_var[variables_in+[new_column]])??#?確定自變量，并在最左邊加上一列全為1的數據，使得模型矩陣中包含截距model_stepwise=sm.OLS(y_var,x_var_cons).fit()??#?用未標準化數據擬合模型：如自變量為x_var_cons，#?則擬合含截距模型；如自變量為x_var，則擬合不含截距模型pvalues[new_column]=model_stepwise.pvalues[new_column]??#?將剛進入模型的變量的回歸系數的p值存儲到pvalues列表（不含截距的p值）sig_pvalue=pvalues.min()??#?獲取pvalues列表中的最小值if?sig_pvalue<pvalue_in:flag=True?#?由于存在小于pvalue_in的p值，所以循環得繼續，不可跳出循環，因此，設置flag為Truesig_variable=pvalues.index[pvalues.argmin()]??#?獲取最小p值對應的自變量（argmin返回最小值的坐標，#?如遇多個最小值，則返回第1個最小值的坐標）variables_in.append(sig_variable)??#?將最小p值對應的自變量添加到用于儲存模型中的變量的variables_in列表#?--------#?制作系數表格（含變量名稱、非標準化回歸系數、t值、p值）x_var_cons=sm.add_constant(x_var[variables_in])model_stepwise=sm.OLS(y_var,x_var_cons).fit()??#?用未標準化數據擬合模型：如自變量為x_var_cons，#?則擬合含截距模型；如自變量為x_var，則擬合不含截距模型print('variable in:?{:10} Adj.R-squared：?{:5}?方程顯著性檢驗p值：?{:5}'.format(sig_variable,model_stepwise.rsquared_adj,model_stepwise.f_pvalue))coefficient=(pd.DataFrame({'未標準化回歸系數':model_stepwise.params}).assign(t值=model_stepwise.tvalues)??#?增加“t值"列.assign(顯著性=model_stepwise.pvalues)??#?增加”顯著性“列)display(coefficient)#?---------#?剔除p值最大且大于pvalue_out的自變量x_var_cons=sm.add_constant(x_var[variables_in])??#?確定自變量，并在最左邊加上一列全為1的數據，使得模型矩陣中包含截距model_stepwise=sm.OLS(y_var,x_var_cons).fit()??#?用未標準化數據擬合模型：如自變量為x_var_cons，則擬合含截距模型；#?如自變量為x_var，則擬合不含截距模型pvalues=model_stepwise.pvalues.iloc[1:]??#?除截距p值外的所有自變量的p值notsig_pvalue=pvalues.max()??#?獲取pvalues列表中的最大值if?notsig_pvalue>pvalue_out:flag=True?#?由于存在大于pvalue_out的p值，所以循環得繼續，不可跳出循環，因此，設置flag為Truenotsig_variable=pvalues.index[pvalues.argmax()]??#?獲取最大p值對應的自變量（argmax返回最大值的坐標，如遇多個最大值，#?則返回第1個最大值的坐標）variables_in.remove(notsig_variable)??#?將最大p值對應的自變量移出用于儲存模型中的變量的variables_in列表#?-------#?制作系數表格（含變量名稱、非標準化回歸系數、t值、p值）x_var_cons=sm.add_constant(x_var[variables_in])model_stepwise=sm.OLS(y_var,x_var_cons).fit()??#?用未標準化數據擬合模型：如自變量為x_var_cons，則擬合含截距模型；#?如自變量為x_var，則擬合不含截距模型print('variable out:?{:10} Adj.R-squared：?{:5}?方程顯著性檢驗p值：?{:5}'.format(notsig_variable,model_stepwise.rsquared_adj,model_stepwise.f_pvalue))coefficient=(pd.DataFrame({'未標準化回歸系數':model_stepwise.params}).assign(t值=model_stepwise.tvalues)??#?增加“t值"列.assign(顯著性=model_stepwise.pvalues)??#?增加”顯著性“列)display(coefficient)if?not?flag:??#?當flag為False時（既沒有需要進入模型的變量，也沒有需要剔除的變量），跳出循環breakreturn?variables_in

4.數據準備

4.1原始數據導入

df=pd.read_csv('../../數據/移動通信網絡的客戶價值分析.csv',encoding='gbk')?? #?數據讀取中文不能正確顯示時，使用encoding #?完整代碼及數據集獲取：@公眾號：數據STUDIO 后臺回復【data】 display(df)

4.2數據記錄去重

if?df.duplicated().sum()>0:display(df[df.duplicated()])??#?顯示重復數據記錄print('數據集中存在以上{}條重復數據記錄，現已刪除。'.format(df[df.duplicated()].shape[0]))df.drop_duplicates(inplace=True)??#?刪除所有字段都完全重復的數據并立即生效display(df) else:print('數據集中不存在重復數據記錄，無需去重。')數據集中不存在重復數據記錄，無需去重。

5.探索分析

5.1數據的描述分析

df.describe(percentiles=(0.01,0.25,0.5,0.75,0.99)).T

結果解讀：

每個變量都有1123條數據記錄，不存在缺失值，未發現明顯的極端值，所有字段都不存在違背業務邏輯的數據記錄。

第一，就因變量“利潤環比增長率”而言，無論是樣本均值還是中位數，都還不錯，一個是19.3%，另一個是18.7%，這顯示校園網計劃獲得初步成功，推薦者確實為校園網帶來了正的相對利潤。但是，從標準差來看，差異性非常大，高達為13.2%，有的推薦者所推薦客戶的利潤環比增長率上升非常大，最大值達到99.3%，而有的推薦者所推薦客戶的利潤環比增長率下降非常大，最大跌幅高達49.8%。正是有這么大的差異，所以本研究就顯得非常必要。

第二，就自變量“月通話總量”而言，推薦者推薦當月的通話總量在258分鐘左右，最大值為360分鐘，最小值為78分鐘，顯示無異常數據。

第三，就自變量“大網占比”而言，其平均值為84.6%，中位數為89.8%，說明所有推薦者都是該運營商深度套牢用戶，因為他們主要通信社交圈都發生在該運營商的網絡內。

第四，就自變量“小網占比”而言，其平均值為25.2%，中位數為19.8%，說明推薦者的通信社交圈為校園網覆蓋的比例并不高，還有很大的發展潛力。

5.2數據的分布形態

查看數據分布情況，有助于特征工程根據數據分布選擇合適的數據處理辦法（包括缺失值、異常值處理，連續特征離散化），還有助于深入了解用戶行為。

對于連續數據，當偏度系數等于0時，數據呈左右對稱分布；當偏度系數絕對值大于等于1時，數據呈嚴重偏斜分布；當偏度系數絕對值大于等于0.5并且小于1時，數據呈中等偏斜分布；當偏度系數絕對值大于0并且小于0.5時，數據呈輕微偏斜分布。

對于分類數據，主要觀察柱狀圖的左右對稱性。

本數據集全為連續數據，沒有分類數據。

#?偏度系數分析 skw_analysis=(pd.DataFrame({'偏度系數':df.skew(numeric_only=True)})??#?只計算數值型字段的偏度系數.assign(偏斜程度=df.skew(numeric_only=True).to_frame()[0].apply(lambda?z:'嚴重'?if?abs(z)>=1?else?'中等'?if?abs(z)>=0.5?else?'輕微'))??#?新增“偏斜程度”列.assign(偏斜方向=df.skew(numeric_only=True).to_frame()[0].apply(lambda?z:'左偏'?if?z<0?else?'對稱'?if?z==0?else?'右偏'))??#?新增“偏斜方向”列).sort_values(by='偏度系數',ascending=True)??#?按“偏度系數”變量升序排序 ) display(skw_analysis)#?數據分布可視化 plt.figure(figsize=(16,8)) rows?=?2 cols?=?4 gs=plt.GridSpec(rows,cols) for?col?in?range(cols):colnames?=?['利潤環比增長率',?'月通話總量',?'大網占比',?'小網占比']label?=?colnames[col]?if?colnames[col]?!=?'月通話總量'?else?colnames[col]?+?'（分鐘）'plt.subplot(gs[0,col])plt.hist(df[colnames[col]],edgecolor='white',bins=100)plt.xlabel(colnames[col]?if?colnames[col]?!=?'月通話總量'?else?colnames[col]?+?'（分鐘）')plt.ylabel('頻數')if?col?==?2:plt.subplot(gs[1,?col],?sharey=plt.subplot(gs[1,3]))else:plt.subplot(gs[1,?col])plt.boxplot(df[col],?labels=label,showmeans=True,?showfliers=True)

結果解讀：

• 第一，字段“大網占比”數據呈嚴重左偏分布，存在極小值。

• 第二，字段“小網占比”數據呈中等右偏分布，存在極大值。

• 第三，字段“月通話總量”、“利潤環比增長率”數據分別呈輕微左偏、輕微右偏分布，分布形態近似對稱，但是都存在極大值和極小值，后續在數據清洗時需進行處理（本次商業數據分析暫時不進行數據清洗），不然會影響模型的穩定性。

5.3數據的正態檢驗

檢驗數據是否服從正態分布的方法比較多，不同檢驗方法對樣本量的敏感度不一樣：樣本量n<50時，優先使用Shapiro-Wilk檢驗；50≤樣本量n<5000時，酌情使用W檢驗及K-S檢驗；5000≤樣本量n時，建議使用K-S檢驗。

不過，我更推薦使用D’Agostino and Pearson omnibus normality test，因為這是一種通用和強大的正態性檢驗方法，其基本思想：首先，計算偏度和峰度以便在不對稱和形狀方面量化分布離正態分布的距離：然后，計算這些值中的每一個與正態分布的預期值之間的差異，并基于這些差異的總和，計算各P值。

stats.normaltest(df)NormaltestResult( statistic=array( [101.1282882 , 66.29959736, 318.22271185, 124.7065164 ]), pvalue=array([1.09716321e-22, 4.01074783e-15,7.92165846e-70, 8.32384597e-28]))

結果解讀：

由于所有P值都小于0.05，所以都拒絕原假設，認為這4個變量都不服從正態分布。

樣本數據不服從正態分布是常態，但只要其分布不存在明顯的偏態或取值不存在異常值，一般無需對數據進行處理，如果一定想處理的話，一般采用取自然對數的方法進行處理。

5.4數據的相關性分析

#?繪制配對圖 sns.pairplot(df)?? #sns.pairplot(df,x_vars=['月通話總量','大網占比','小網占比'],y_vars='利潤環比增長率',height=4,aspect=1)#?計算pearson相關系數 df.corr(method='pearson')

結果解讀：

• 第一，配對圖顯示，利潤環比增長率與月通話總量存在線性相關。

• 第二，相關系數顯示，利潤環比增長率與月通話總量存在顯著線性正相關（r=-0.7379），大網占比與小網占比存在顯著負相關（r=-0.6557）。

5.5數據的箱線圖分析

由于本案例的4個變量都是連續變量，因此，優先選擇散點圖來嘗試探索變量間的大概關系(正相關、負相關)，但是，對于大多數實際問題而言，數據噪聲很大，很難從散點圖上清晰地發現規律，所以，經常使用分組箱線圖進行探索。

對于本案例，先將"利潤環比增長率"因變量按升序排序，把其前27%的數據劃為低價值組（303個樣本），把其后27%的數據劃為高價值組（303個樣本），然后對每一個自變量做分組箱形圖，分組箱形圖可以直觀的展示出因變量與自變量之間的關系。

df=df.assign(推薦者價值=df.利潤環比增長率.map(lambda?z:'高'?if?z>=df.利潤環比增長率.nlargest(303).min()?else?'低'?if?z<=df.利潤環比增長率.nsmallest(303).max()?else?None))#df.boxplot(column=['大網占比','小網占比'],by='推薦者價值',layout=(1,3),figsize=(12,4),grid=False)??#?縱坐標的量綱需一致，否則箱子會變形，因為所有圖形共享一個縱坐標 plt.figure(figsize=(12,4)) gs=plt.GridSpec(1,3) plt.subplot(gs[0,0]) grouping_boxplot('推薦者價值','推薦者價值','月通話總量','月通話總量（分鐘）',data=df) plt.subplot(gs[0,1]) grouping_boxplot('推薦者價值','推薦者價值','大網占比','大網占比',data=df) plt.subplot(gs[0,2],sharey=plt.subplot(gs[0,1]))??#?共享gs[0,1]的縱坐標 grouping_boxplot('推薦者價值','推薦者價值','小網占比','小網占比',data=df)

結果解讀：

• 第一，高價值推薦者的月通話總量比低價值推薦者高。

• 第二，高價值推薦者的大網占比比低價值推薦者高。

• 第三，高價值推薦者的小網占比比低價值推薦者低。

6.回歸分析

6.1模型構建

6.1.1無交互效應模型

1、輸入法篩選自變量

x_enter=df[['月通話總量',?'大網占比',?'小網占比']]??#?確定自變量數據 y_enter=df.利潤環比增長率??#?確定因變量數據 X_enter=sm.add_constant(x_enter)??#?加上一列全為1的數據，使得模型矩陣中包含截距 model_enter=sm.OLS(y_enter,X_enter).fit()??#?用未標準化數據擬合模型：#?X大寫則擬合含截距模型，#?x小寫則擬合不含截距項模型，#?標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 print(model_enter.summary())

結果解讀：

自變量“小網占比”未通過顯著性檢驗（p=0.638>0.05），最終建模時需刪除“小網占比”自變量。

2、步進法篩選自變量

x_step=df[['月通話總量',?'大網占比',?'小網占比']]??#?確定自變量全集數據 y_step=df.利潤環比增長率??#?確定因變量數據 x_step=df[stepwise_select_variable(x_step,y_step)]??#?使用步進法篩選出的自變量，stepwise_select_variable函數的功能是使用步進法篩選自變量，并返回最終篩選出的自變量名稱 X_step=sm.add_constant(x_step)??#?加上一列全為1的數據，使得模型矩陣中包含截距 model_step=sm.OLS(y_step,X_step).fit()??#?用未標準化數據擬合模型：X大寫則擬合含截距模型，x小寫則擬合不含截距項模型，標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 print(model_step.summary())variable in: 月通話總量 Adj.R-squared：0.54414892458074 方程顯著性檢驗p值：1.1532634227093457e-193variable in: 大網占比 Adj.R-squared：0.5562037611861936 方程顯著性檢驗p值：9.74806646667517e-199

結果解讀：

最終只有“月通話總量”和“大網占比”進入模型，“小網占比”被剔除出模型，和使用“輸入法”選擇自變量的結果一致。

3、子集法篩選自變量

x_all=df[['月通話總量',?'大網占比',?'小網占比']]??#?確定包含所有自變量的數據 y_all=df.利潤環比增長率??#?確定因變量數據 subset_score=pd.DataFrame(columns=['independent_variable','Rsquare','Rsquare_adj','AIC','BIC']) i=0??#?用于確定往subse_score數據集的第幾行添加數據 for?x_quantity?in?range(1,x_all.shape[1]+1):??#?確定自變量組合中自變量的數量for?x_com?in?combinations(x_all.columns,x_quantity):??#?遍歷自變量數量為x_quantity的各種自變量組合x=df[list(x_com)]X=sm.add_constant(x)model_all=sm.OLS(y_all,X).fit()subset_score.loc[i]=[x_com,model_all.rsquared,model_all.rsquared_adj,model_all.aic,model_all.bic]i+=1

通過Rsquare_adj選擇自變量子集。

subset_score.loc[subset_score.Rsquare_adj.argmax(),'independent_variable']? #?選擇rsquare_adj取值最大的子集進入最終模型，argmax返回最大值對應的索引號('月通話總量', '大網占比')

通過AIC準則選擇自變量子集。

subset_score.loc[subset_score.AIC.argmin(),'independent_variable']?? #?選擇rsquare_adj取值最大的子集進入最終模型，argmax返回最大值對應的索引號('月通話總量', '大網占比')

通過BIC準則選擇自變量子集。

subset_score.loc[subset_score.BIC.argmin(),'independent_variable']?? #?選擇rsquare_adj取值最大的子集進入最終模型，argmax返回最大值對應的索引號('月通話總量', '大網占比')

結果解讀：

經過對自變量的所有子集進行擬合，使用Rsquare_adj、AIC準則、BIC準則對模型進行篩選，所篩選出的自變量都是“月通話總量”和“大網占比”兩個。

結論：

“輸入法”、“步進法”、“子集法”等三種自變量篩選法所篩選出的自變量都是“月通話總量”和“大網占比”兩個，因此，最終進入模型的自變量為“月通話問題”和“大網占比”。

6.1.2有交互效應模型

x_interact=df[['月通話總量','大網占比']].assign(月通話總量x大網占比=df.月通話總量*df.大網占比) X_interact=sm.add_constant(x_interact) y_interact=df.利潤環比增長率 model_interact=sm.OLS(y_interact,X_interact).fit() print(model_interact.summary())

結果解讀：

由于交互項“月通話總量x大網占比”不顯著（p=0.590>0.05），因此，不存在交互效應。

6.2模型診斷

6.2.1多重共線性診斷

#?計算VIF IDV=['月通話總量','大網占比','小網占比']??#?確定自變量名稱全集，IDV=independent?variable vifs=[] for?variable?in?IDV:x_vif=df[list(set(IDV)-{variable})]X_vif=sm.add_constant(x_vif)y_vif=df[variable]model_vif=sm.OLS(y_vif,X_vif).fit()??#?用未標準化數據擬合模型：X大寫則擬合含截距模型，x小寫則擬合不含截距項模型，標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致vif=1/(1-model_vif.rsquared)vifs.append(vif) VIFS=pd.DataFrame(index=IDV).assign(VIF=vifs).sort_values(by='VIF',ascending=False)??#?按各變量vif值降序排序 display(VIFS)

結果解讀：

當VIF<5時，回歸方程存在輕度多重共線性；當5≤VIF<10時，回歸方程存在較嚴重的多重共線性；當10≤VIF時，回歸方程存在嚴重的多重共線性。

因此，自變量之間不存在多重共線性問題。

6.2.2方差不齊性診斷

1、殘差散點圖法

x_hetero=df[['月通話總量',?'大網占比',?'小網占比']]??#?確定自變量數據 y_hetero=df.利潤環比增長率??#?確定因變量數據 X_hetero=sm.add_constant(x_hetero)??#?加上一列全為1的數據，使得模型矩陣中包含截距 model_hetero=sm.OLS(y_hetero,X_hetero).fit()??#?用未標準化數據擬合模型：X大寫則擬合含截距模型，x小寫則擬合不含截距項模型，標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 df=(df.assign(resid=model_hetero.resid)??#?新增“resid”列，表示未標準化殘差.assign(sdresid=model_hetero.outlier_test().student_resid)??#?新增”sdresid“列，表示學生化刪除后殘差(等價于model_hetero.get_influence().resid_studentized_external) ) sns.pairplot(df,x_vars=['月通話總量','大網占比','小網占比'],y_vars=['resid','sdresid'],height=4,aspect=1)??#?繪制殘差散點圖

結果解讀：

散點圖顯示，未標準化殘差及學生化刪除后殘差與三個自變量都不存在明顯的線性關系，因此，不存在方差不齊性。

2、等級相關系數法

IDV=['月通話總量','大網占比','小網占比']??#?確定自變量名稱全集，IDV=independent?variable for?variable?in?IDV:print('未標準化殘差與“{}”變量的spearman等級相關系數為{}，p值為{}。'.format(variable,stats.spearmanr(abs(df.resid),df[variable])[0],stats.spearmanr(abs(df.resid),df[variable])[1]))print('學生化刪除后殘差與“{}”變量的spearman等級相關系數為{}，p值為{}。'.format(variable,stats.spearmanr(abs(df.sdresid),df[variable])[0],stats.spearmanr(abs(df.sdresid),df[variable])[1]))未標準化殘差與“月通話總量”變量的spearman等級相關系數為 0.031471496792663256，p值為0.2920047768383987。 學生化刪除后殘差與“月通話總量”變量的spearman等級相關系數為 0.031522051641127945，p值為0.29122989772460245。 未標準化殘差與“大網占比”變量的spearman等級相關系數為 -0.004794315991555318，p值為0.8724988457644585。 學生化刪除后殘差與“大網占比”變量的spearman等級相關系數為 -0.0049131032889456515，p值為0.8693677868252431。 未標準化殘差與“小網占比”變量的spearman等級相關系數為 0.01220849987834211，p值為0.6827737839962207。 學生化刪除后殘差與“小網占比”變量的spearman等級相關系數為 0.012965747385529707，p值為0.6642651943572604。

結果解讀：

未標準化殘差、學生化刪除后殘差與三個自變量的spearman等級相關系數均不顯著（p值都大于0.05），因此，認為不存在方差不齊性。

3、BP檢驗

sm.stats.diagnostic.het_breuschpagan(model_hetero.resid,exog_het=model_hetero.model.exog)(3.4984809632531215, 0.32095918822451, 1.1656378996396506, 0.3216685113072825)

結果解讀：

第一個值是LM統計量，第二個值是LM統計量對應的p值，結果說明接受殘差方差為常數的原假設（p>0.05）；第三個值是F統計量，用于檢驗殘差平方與自變量之間是否獨立，如果獨立則說明殘差方差齊性，第四個值為F統計量對應的p值，說明殘差項滿足方差齊性（p>0.05）。

6.2.3序列相關性診斷

1、自相關系數和偏自相關系數

x_autocorr=df[['月通話總量',?'大網占比',?'小網占比']]??#?確定自變量數據 y_autocorr=df.利潤環比增長率??#?確定因變量數據 X_autocorr=sm.add_constant(x_autocorr)??#?加上一列全為1的數據，使得模型矩陣中包含截距 model_autocorr=sm.OLS(y_autocorr,X_autocorr).fit()??#?用未標準化數據擬合模型：X大寫則擬合含截距模型，x小寫則擬合不含截距項模型，標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 #?計算自相關系數和偏自相關系數 pd.DataFrame({'自相關系數':sm.tsa.stattools.acf(model_autocorr.resid,nlags=20),'偏自相關系數':sm.tsa.stattools.pacf(model_autocorr.resid,nlags=20,method='ywm')})??#?采用Yule-Walker方法，結果與SPSS結果一致上下滑動查看更多#?繪制自相關系數折線圖 plt.figure(figsize=(12,4)) autocorrelation_plot(model_autocorr.resid) plt.xticks(list(range(23))) plt.title('自相關系數折線圖') plt.xlabel('滯后階數') plt.ylabel('自相關系數')

結果解讀：

自相關系數的取值區間為[-0.059~0.079]，都非常小，故認為不存在序列相關性。

#?繪制自相關圖和偏自相關圖 fig,axes=plt.subplots(nrows=1,ncols=2,figsize=(16,5),dpi=80) plot_acf(model_autocorr.resid,ax=axes[0],lags=20)??#?繪制自相關圖 axes[0].set_title('自相關圖')??#?設置自相關圖標題，也可不設置，采用默認值Autocorrelation axes[0].set_ylabel('自相關系數') axes[0].set_xlabel('滯后階數') axes[0].set_xticks(list(range(21)))plot_pacf(model_autocorr.resid,lags=20,ax=axes[1],method='ywm')??#?采用Yule-Walker方法，結果與SPSS結果一致 axes[1].set_title('偏自相關圖')?#設置偏自相關圖標題，也可不設置，采用默認值Partial?Autocorrelation axes[1].set_ylabel('偏自相關系數') axes[1].set_xlabel('滯后階數') axes[1].set_xticks(list(range(21)))

結果解讀：

自相關圖和偏自相關圖，反映了殘差序列的各階自相關系數的大小，該圖的高度值對應的是各階自相關系數的值，藍色區域是95%置信區間，這兩條界線是檢測自相關系數是否為0時所使用的判別標準：當代表自相關系數的柱條超過這兩條界線時，可以認定自相關系數顯著不為0。

觀察上圖可知，基本上所有的點都落在95%置信區間內，所以初步判斷不存在序列相關性。

2、滯后殘差散點圖

#?繪制t期殘差與其滯后1-4階的圖 def?draw_resid():plt.figure(figsize=(16,8))n,?rows,?cols?=?0,?2,?4gs=plt.GridSpec(rows,cols)for?row?in?range(rows):for?col?in?range(cols):n?+=?1plt.subplot(gs[row,col])lag_plot(model_autocorr.resid,lag=n)??#?lag表示滯后期數draw_resid()

結果解讀：

由圖可知，散點未呈現任何線性關系，故認為不存在序列相關性。

3、殘差時序圖

plt.scatter(df.index,model_autocorr.resid) plt.xlabel('索引號') plt.ylabel('殘差')

結果解讀：

由圖可知，散點未呈現任何規律性，故認為不存在序列相關性。

4、Durbin-Watson檢驗

print('Durbin-Watson值為{}。'.format(sm.stats.stattools.durbin_watson(model_autocorr.resid)))? #?計算Durbin-Watson值Durbin-Watson值為1.9906362533785131。

結果解讀：

樣本量太多，無法查尋DW檢驗表，故無法通過DW值判斷序列相關性，需使用LM檢驗。不過，由于DW值趨近于2，根據自相關系數的計算公式，可知自相關系數趨近于0，認為不存在序列相關性（一階）。

5、拉格朗日乘數檢驗(又稱為LM檢驗、BG檢驗)

Davidson and MacKinnon(1993)建議，把殘差中因滯后而缺失的項用其期望值0來代替，以保持樣本容量為n。acorr_breusch_godfrey函數采用的是Davidson-MacKinnon方法。

acorr_breusch_godfrey函數的重要參數：

• res：回歸結果，對該模型的殘差進行自相關檢驗，此處應當填寫模型名稱（如model_autocorr）

• nlags：滯后階數

acorr_breusch_godfrey函數的返回結果：

• lm_statistic：LM統計量值

• lm_pvalue：LM統計量的p值，若p值小于顯著性水平，則拒絕無自相關性的原假設，即存在自相關性

• f_statistic：F統計量值（resid_lag1、resid_lag2、resid_lag3、……、resid_lagk聯合顯著的F檢驗統計量）

• f_pvalue：F統計量對應的p值，若p值小于顯著性水平，則拒絕無自相關性的原假設，即存在自相關性

bg_result=acorr_breusch_godfrey(model_autocorr,nlags=20)??#?拉格朗日乘數檢驗，nlags為滯后階數 bg_test_output=pd.Series(bg_result[0:4],index=['lm_statistic','lm_pvalue','f_statistic','f_pvalue']) display(bg_test_output)lm_statistic 24.898569 lm_pvalue 0.205339 f_statistic 1.245947 f_pvalue 0.207554 dtype: float64

結果解讀：

兩個檢驗的p值都大于0.05，故認為不存在序列相關性。

6.3模型結果

6.3.1非標準化回歸模型

由于：

• 第一，根據“輸入法”、“步進法”、“子集法”自變量篩選方法的分析結論，可知最終模型只包括”月通話總量“和”大網占比“兩個自變量。

• 第二，根據模型的多重共線性、方差齊次性、序列相關性診斷結論，可知模型不存在多重共線性、方差齊次性、序列相關性等問題。

因此：

在最終的回歸模型當中，只納入”月通話總量“和”大網占比“兩個自變量，剔除“小網占比”自變量。

x_end_variables=['月通話總量','大網占比']??#?確定最終進入模型的自變量的名稱 #y_end_variable=['利潤環比增長率']??#?確定最終進入模型的因變量的名稱 x_end=df[x_end_variables]??#?確定最終進入模型的自變量數據 y_end=df['利潤環比增長率']??#?確定最終進入模型的因變量數據 scaler=StandardScaler()??#?建模：創建數據標準化模型 x_end_std=scaler.fit_transform(x_end)??#?標準化x：fit-用數據訓練模型；# transform-用訓練好的模型對數據進行轉換；#?fit_transform-用數據訓練模型，并用訓練好的模型對數據進行轉換（二合一） y_end_std=scaler.fit_transform(np.array(y_end).reshape(-1,1))??#?標準化y X_end=sm.add_constant(x_end)??#?加上一列全為1的數據，使得模型矩陣中包含截距 X_end_std=sm.add_constant(x_end_std)??#?加上一列全為1的數據，使得模型矩陣中包含截距 model_end=sm.OLS(y_end,X_end).fit()??#?用未標準化數據擬合模型：#?X大寫則擬合含截距模型，#?x小寫則擬合不含截距項模型，#?標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 model_end_std=sm.OLS(y_end_std,X_end_std).fit()??#?用標準化數據擬合模型：#?X大寫則擬合含截距模型，#?x小寫則擬合不含截距項模型，#?標準化與未標準化模型的x大小寫需一致，否則檢驗統計量取值將不一致 print(model_end.summary())??#?非標準化回歸模型摘要

6.3.2標準化回歸模型

print(model_end_std.summary())??#?標準化回歸模型摘要

6.3.3結果呈現表格

#?繪制結果呈現表格 pd.DataFrame({'coef':round(model_end.params,4),'Std.Error':round(model_end.bse,4),'tvalue':round(model_end.tvalues,3),'pvalue':round(model_end.pvalues,3),'lower?95%?CI':round(model_end.conf_int(0.05)[0],3),'upper?95%?CI':round(model_end.conf_int(0.05)[1],3),'std?coef':model_end_std.params.round(4),'VIF':VIFS.loc[x_end_variables].VIF})

6.4模型應用

6.4.1客戶群體細分

根據模型結果，最有價值的推薦者應具備兩個特征：第一，月通話總量高，這意味著推薦者是個“話癆”；第二，大網占比高，這意味著推薦者的主要移動通信社交圈都發生在該運營商的網絡內。

對于這兩個特征，使用聚類分析法將“月通話總量”變量分為“高”和“低”兩類，將“大網占比”變量也分為“高”和“低”兩類，于是便將客戶群細分為4類：第一類，高月通話總量、高大網占比；第二類，高月通話總量、低大網占比；第三類，低月通話總量、高大網占比；第四類，低月通話總量、低大網占比。

1、單變量細分

KMeans算法的基本原理

KMeans算法名稱中的K代表類別數量，Means代表每個類別內樣本的均值，所以KMeans算法又稱為K-均值算法。KMeans算法以距離作為樣本間相似度的度量標準，將距離相近的樣本分配至同一個類別。樣本間距離的計算方式可以是歐氏距離、曼哈頓距離、余弦相似度等，KMeans算法通常采用歐氏距離來度量各樣本間的距離。

KMeans算法的核心思想是對每個樣本點計算到各個中心點的距離，并將該樣本點分配給距離最近的中心點代表的類別，一次迭代完成后，根據聚類結果更新每個類別的中心點，然后重復之前操作再次迭代，直到前后兩次分類結果沒有差別。

基本操作步驟：

• 步驟1：隨機選擇K個中心點。

• 步驟2：依據歐氏距離度量相似度，將每個樣本點都分配給最近的中心點。

• 步驟3：根據聚類結果，更新每個類別的中心點。

• 步驟4：重復步驟2，將每個樣本點重新分配給距離最近的中心點。

• 步驟5：重復步驟3，更新中心點。

• 步驟6：重復步驟4，判斷聚類結果是否與上次結果一致，若一致則算法終止，否則回到步驟5。

kms=KMeans(n_clusters=2,random_state=20220923)?? #?建模：設置模型參數：n_clusters表示將樣本聚為2類；# random_state設置種子參數，使得每次運行代碼時選取的初始中心點都一樣。#?增加”月通話總量分類“和”大網占比分類“新列 df['月通話總量分類']=kms.fit(np.array(df.月通話總量).reshape(-1,1)).labels_?? #?擬合：用數據擬合（訓練）剛剛所建模型，labels_屬性返回聚類結果 df['大網占比分類']=kms.fit(np.array(df.大網占比).reshape(-1,1)).labels_?? #?擬合：用數據擬合（訓練）剛剛所建模型，labels_屬性返回聚類結果 display(df)#?確定分類代碼的業務含義 print(df.groupby('月通話總量分類').月通話總量.mean()) print(df.groupby('大網占比分類').大網占比.mean())月通話總量分類 0 2.888206 1 2.259800 Name: 月通話總量, dtype: float64 大網占比分類 0 0.615438 1 0.915795 Name: 大網占比, dtype: float64

結果解讀：

• 第一，在對”月通話總量“變量進行分類時，0代表”高“，1代表”低“。

• 第二，在對”大網占比“變量進行分類時，0代表”低“，1代表”高“。

#?把數字0和1換成具有業務含義的”高“和”低“ df.月通話總量分類.replace([0,1],['高','低'],inplace=True) df.大網占比分類.replace([0,1],['低','高'],inplace=True)

2、雙變量細分

#?計算各類的數量 crosstab_quantity=pd.crosstab(df.大網占比分類,df.月通話總量分類) display(crosstab_quantity)#?計算各類的占比 crosstab_percent=pd.crosstab(df.大網占比分類,df.月通話總量分類,normalize=True) display(crosstab_percent)#?繪制客戶群體細分表 crosstab=pd.DataFrame(index=['大網占比_低','大網占比_高'],columns=['月通話總量_低','月通話總量_高']) crosstab.iloc[0,0]='{}({:.2%})[D.最劣質客戶-放棄]'.format(crosstab_quantity.iloc[0,0],crosstab_percent.iloc[0,0]) crosstab.iloc[0,1]='{}({:.2%})[C.高潛力客戶-策反]'.format(crosstab_quantity.iloc[0,1],crosstab_percent.iloc[0,1]) crosstab.iloc[1,0]='{}({:.2%})[B.待激活客戶-激活]'.format(crosstab_quantity.iloc[1,0],crosstab_percent.iloc[1,0]) crosstab.iloc[1,1]='{}({:.2%})[A.高價值客戶-套牢]'.format(crosstab_quantity.iloc[1,1],crosstab_percent.iloc[1,1]) display(crosstab)plt.scatter(df[(df.大網占比分類=='高')&(df.月通話總量分類=='高')].大網占比,df[(df.大網占比分類=='高')&(df.月通話總量分類=='高')].月通話總量,label='高價值客戶',c='blue',marker='o') plt.scatter(df[(df.大網占比分類=='低')&(df.月通話總量分類=='低')].大網占比,df[(df.大網占比分類=='低')&(df.月通話總量分類=='低')].月通話總量,label='最劣質客戶',c='red',marker='^') plt.scatter(df[(df.大網占比分類=='高')&(df.月通話總量分類=='低')].大網占比,df[(df.大網占比分類=='高')&(df.月通話總量分類=='低')].月通話總量,label='待激活客戶',c='black',marker='*') plt.scatter(df[(df.大網占比分類=='低')&(df.月通話總量分類=='高')].大網占比,df[(df.大網占比分類=='低')&(df.月通話總量分類=='高')].月通話總量,label='高潛力客戶',c='green',marker='x') plt.xlabel('大網占比') plt.ylabel('月通話總量（百分鐘）') plt.legend()

6.4.2營銷策略建議

根據前面的模型分析結果可知，最有價值的推薦者應當具備兩個特征：第一，月通話總量高，這意味著推薦者是個話癆；第二，大網占比高，這意味著推薦者的主要通信社交圈都發生在該運營商的網絡內。

基于這兩個特征，運用K均值聚類分析法，將自變量”月通話總量“、”大網占比“分別分成2組，并構建客戶群體細分表，于是，就把推薦者分成4類，對每類推薦者，我們可以制定針對類營銷策略。

對于群體A，對應的是最優質的高價值客戶，占樣本比為40.52%，其特征是高通話量、高大網占比。因此，對這部分客戶的營銷策略是向其提供最優質的營銷資源和最好的客戶服務，將他們深度的套牢，同時，還應當向他們提供最激進的獎勵機制，鼓勵他們將大網中的好友發展為校園網用戶。

對于群體B，對應的是待激活客戶，占樣本比為36.15%，其特征是低通話總量、高大網占比。這部分客戶的優勢是大網占比高，這說明他們的通信社交圈被同一個運營商大量覆蓋，他們的的劣勢是通話總量不高，活躍度低。如果通過合理的營銷策略刺激他們提高通話量，有可能其通話對象的通話量也會被大大提高。

對于群體C，對應的是高法力客戶，占樣本10.95%，其特征是高通話總量、低大網占比。這部分客戶的優勢是通話量高，活躍度高，劣勢是大網占比低。這說明他們有大量高價值通信對象未被該運營商覆蓋。因此，可考慮為這批客戶設立激進的推薦獎勵機制，鼓勵他對身邊好友進行策反，邀請他身邊好友加入該運營商，以提高大網占比。

對于群體D，對應的是最劣質客戶，占樣本12.38%，其特征是低通話總量、低大網占比。他們既不活躍，又有大量的通信對象未被該運營商覆蓋，是最劣質的一批用戶，在營銷資源有限的情況下，他們或許可以被暫時放棄。

-?END -

對比Excel系列圖書累積銷量達15w冊，讓你輕松掌握數據分析技能，可以點擊下方鏈接進行了解選購：

總結

以上是生活随笔為你收集整理的移动通信客户价值数据挖掘分析实战的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。