1、如何認識可視化？
圖形總是比數據更加醒目、直觀。解決統計回歸問題，無論在分析問題的過程中，還是在結果的呈現和發表時，都需要可視化工具的幫助和支持。
需要指出的是，雖然不同繪圖工具包的功能、效果會有差異，但在常用功能上相差並不是很大。與選擇哪種繪圖工具包相比，更重要的是針對不同的問題，需要思考選擇什么方式、何種圖形去展示分析過程和結果。換句話說，可視化只是手段和形式，手段要為目的服務，形式要為內容服務，這個關系一定不能顛倒了。
因此，可視化是伴隨着分析問題、解決問題的過程而進行思考、設計和實現的，而且還會影響問題的分析和解決過程：
-
可視化工具是數據探索的常用手段
回歸分析是基於數據的建模，在導入數據后首先要進行數據探索，對給出的或收集的數據有個大概的了解，主要包括數據質量探索和數據特征分析。數據准備中的異常值分析，往往就需要用到箱形圖（Boxplot）。對於數據特征的分析，經常使用頻率分布圖或頻率分布直方圖（Hist），餅圖（Pie）。
-
分析問題需要可視化工具的幫助
對於問題中變量之間的關系，有些可以通過定性分析來確定或猜想，需要進一步的驗證，有些復雜關系難以由分析得到，則要通過對數據進行初步的相關分析來尋找線索。在分析問題、嘗試求解的過程中，雖然可以得到各種統計量、特征值，但可視化圖形能提供更快捷、直觀、豐富的信息，對於發現規律、產生靈感很有幫助。
-
解題過程需要可視化工具的支持
在解決問題的過程中，也經常會希望盡快獲得初步的結果、總體的評價，以便確認解決問題的思路和方法是否正確。這些情況下，我們更關心的往往是繪圖的便捷性，圖形的表現效果反而是次要的。
-
可視化是結果發布的重要內容
問題解決之后需要對結果進行呈現或發表，這時則需要結合表達的需要，特別是表達的邏輯框架，設計可視化的方案，選擇適當的圖形種類和形式，准備圖形數據。在此基礎上，才談得上選擇何種繪圖工具包，如何呈現更好的表現效果。
歡迎關注 Youcans 原創系列，每周更新數模筆記
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計回歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法
2、StatsModels 繪圖工具包 （Graphics）
Statsmodels 本身支持繪圖功能（Graphics），包括擬合圖（Fit Plots）、箱線圖（Box Plots）、相關圖（Correlation Plots）、函數圖（Functional Plots）、回歸圖（Regression Plots）和時間序列圖（Time Series Plots）。
Statsmodels 內置繪圖功能 Graphics 的使用似乎並不流行，網絡上的介紹也不多。分析其原因，一是 Graphics 做的並不太好用，文檔和例程不友好，二是學習成本高：能用通用的可視化包實現的功能，何必還要花時間去學習一個專用的 Graphics？
下面是 Statsmodels 官方文檔的例程，最簡單的單變量線性回歸問題，繪制樣本數據散點圖和擬合直線圖。Graphics 提供了將擬合與繪圖合二為一的函數 qqline()，但是為了繪制出樣本數據則要調用 Matplotlib 的 matplotlib.pyplot.scatter()，所以...
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.graphics.gofplots import qqline
foodexp = sm.datasets.engel.load(as_pandas=False)
x = foodexp.exog
y = foodexp.endog
ax = plt.subplot(111)
plt.scatter(x, y)
qqline(ax, "r", x, y)
plt.show()
右圖看起來有點象 Seaborn中的 relplot，但把官方文檔研究了半天也沒搞明白，只好直接分析例程和數據，最后的結論是：基本沒啥用。
這大概就是更多用戶直接選擇 Python 的可視化工具包進行繪圖的原因吧。最常用的當屬 Matplotlib 無疑，而在統計回歸分析中 Seaborn 繪圖工具包則更好用更炫酷。
3、Matplotlib 繪圖工具包
Matplotlib 繪圖包就不用介紹了。Matplotlib 用於 Statsmodels 可視化，最大的優勢在於Matplotlib 誰都會用，實現統計回歸的基本圖形的也很簡單。如果需要復雜的圖形，炫酷的效果，雖然 Matplotlib 原理上也能實現，但往往需要比較繁瑣的數據准備，並不常用的函數和參數設置。既然學習成本高，出錯概率大，就沒必要非 Matplotlib 不可了。
Matplotlib 在統計回歸問題中經常用到的是折線圖、散點圖、箱線圖和直方圖。這也是 Matplotlib 最常用的繪圖形式，本系列文中也有相關例程，本文不再具體介紹相關函數的用法。
例如，在本系列《Python學習筆記-StatsModels 統計回歸（2）線性回歸》的例程和附圖，不僅顯示了原始檢測數據、理論模型數據、擬合模型數據，而且給出了置信區間的上下限，看起來還是比較“高級”的。但是，如果把置信區間的邊界線隱藏起來，圖形馬上就顯得不那么“高級”，比較“平常”了——這就是選擇什么方式、何種圖形進行展示的區別。
由此所反映的問題，還是表達的邏輯和數據的准備：要表達什么內容，為什么要表達這個內容，有沒有相應的數據？問題的關鍵並不是什么工具包或什么函數，更不是什么顏色什么線性，而是有沒有置信區間上下限的數據。
如果需要復雜的圖形，炫酷的效果，雖然 Matplotlib 原理上也能實現，但往往需要比較繁瑣的數據准備，使用並不常用的函數和參數設置。學習成本高，出錯概率大，就沒必要非 Matplotlib 不可了。
4、Seaborn 繪圖工具包
Seaborn 是在 Matplotlib 上構建的，支持 Scipy 和 Statamodels 的統計模型可視化，可以實現：
- 賞心悅目的內置主題及顏色主題
- 展示和比較 一維變量、二維變量 各變量的分布情況
- 可視化 線性回歸模型中的獨立變量和關聯變量
- 可視化 矩陣數據，通過聚類算法探究矩陣間的結構
- 可視化 時間序列，展示不確定性
- 復雜的可視化，如在分割區域制圖
Seaborn 繪圖工具包以數據可視化為中心來挖掘與理解數據，本身就帶有一定的統計回歸功能，而且簡單好用，特別適合進行定性分析、初步評價。
下圖給出了幾種常用的 Seaborn 圖形，分別是帶擬合線的直方圖（distplot）、箱線圖（boxplot）、散點圖（scatterplot）和回歸圖（regplot），后文給出了對應的程序。
實際上，這些圖形用 StatsModels Graphics、Matplotlib 也可以繪制，估計任何繪圖包都可以實現。那么，為什么還要推薦 Seaborn 工具包，把這些圖歸入 Seaborn 的實例呢？我們來看看實現的例程就明白了：簡單，便捷，舒服。不需要數據准備和變換處理，直接調用變量數據，自帶回歸功能；不需要復雜的參數設置，直接給出舒服的圖形，自帶圖形風格設計。
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
fig1, axes = plt.subplots(2, 2, figsize=(10, 8)) # 創建一個 2行 2列的畫布
sns.distplot(df['price'], bins=10, ax=axes[0, 0]) # axes[0,1] 左上圖
sns.boxplot(df['price'], df['sales'], data=df, ax=axes[0, 1]) # axes[0,1] 右上圖
sns.scatterplot(x=df['advertise'], y=df['sales'], ax=axes[1, 0]) # axes[1,0] 左下圖
sns.regplot(x=df['difference'], y=df['sales'], ax=axes[1, 1]) # axes[1,1] 右下圖
plt.show()
5、多元回歸案例分析（Statsmodels）
5.1 問題描述
數據文件中收集了 30個月本公司牙膏銷售量、價格、廣告費用及同期的市場均價。
（1）分析牙膏銷售量與價格、廣告投入之間的關系，建立數學模型；
（2）估計所建立數學模型的參數，進行統計分析；
（3）利用擬合模型，預測在不同價格和廣告費用下的牙膏銷售量。
* 本問題及數據來自：姜啟源、謝金星，數學模型（第 3版），高等教育出版社。
5.2 問題分析
本案例在《Python學習筆記-StatsModels 統計回歸（3）模型數據的准備》中就曾出現，文中還提到該文的例程並不是最佳的求解方法和結果。這是因為該文例程是直接將所有給出的特征變量（銷售價格、市場均價、廣告費、價格差）都作為自變量，直接進行線性回歸。謝金星老師說，這不科學。科學的方法是先分析這些特征變量對目標變量（銷量）的影響，然后選擇能影響目標的特征變量，或者對特征變量進行適當變換（如：平方、對數）后，再進行線性回歸。以下參考視頻教程中的解題思路進行分析。
- 觀察數據分布特征
案例問題的數據量很小，數據完整規范，實際上並不需要進行數據探索和數據清洗，不過可以看一下數據的分布特性。例程和結果如下，我是沒看出什么名堂來，與正態分布的差距都不小。
# 數據探索：分布特征
fig1, axes = plt.subplots(2, 2, figsize=(10, 8)) # 創建一個 2行 2列的畫布
sns.distplot(dfData['price'], bins=10, ax=axes[0,0]) # axes[0,1] 左上圖
sns.distplot(dfData['average'], bins=10, ax=axes[0,1]) # axes[0,1] 右上圖
sns.distplot(dfData['advertise'], bins=10, ax=axes[1,0]) # axes[1,0] 左下圖
sns.distplot(dfData['difference'], bins=10, ax=axes[1,1]) # axes[1,1] 右下圖
plt.show()
- 觀察數據間的相關性
既然將所有特征變量都作為自變量直接進行線性回歸不科學，就要先對每個自變量與因變量的關系進行考察。
# 數據探索：相關性
fig2, axes = plt.subplots(2, 2, figsize=(10, 8)) # 創建一個 2行 2列的畫布
sns.regplot(x=dfData['price'], y=dfData['sales'], ax=axes[0,0])
sns.regplot(x=dfData['average'], y=dfData['sales'], ax=axes[0,1])
sns.regplot(x=dfData['advertise'], y=dfData['sales'], ax=axes[1,0])
sns.regplot(x=dfData['difference'], y=dfData['sales'], ax=axes[1,1])
plt.show()
單變量線性回歸圖還是很有價值的。首先上面兩圖（sales-price，sales-average）的數據點分散，與回歸直線差的太遠，說明與銷量的相關性小——謝金星老師講課中也是這樣分析的。其次下面兩圖（sales-advertise，sales-difference）的線性度較高，至少比上圖好多了，回歸直線和置信區間也反映出線性關系。因此，可以將廣告費（advertise）、價格差（difference）作為自變量建模進行線性回歸。
進一步地，有人觀察散點圖后認為銷量與廣告費的關系（sales-advertise）更接近二次曲線，對此也可以通過回歸圖對 sales 與 advertise 進行高階多項式回歸擬合，結果如下圖。
-
建模與擬合
- 模型 1：將所有特征變量都作為自變量直接進行線性回歸，這就是《（3）模型數據的准備》中的方案。
- 模型 2：選擇價格差（difference）、廣告費（advertise）作為自變量建模進行線性回歸。
- 模型 3：選擇價格差（difference）、廣告費（advertise）及廣告費的平方項作為作為自變量建模進行線性回歸。
下段給出了使用不同模型進行線性回歸的例程和運行結果。對於這個問題的分析和結果討論，謝金星老師在視頻中講的很詳細，網絡上也有不少相關文章。由於本文主要講可視化，對結果就不做詳細討論了。
6、Python 例程（Statsmodels）
6.1 問題描述
數據文件中收集了 30個月本公司牙膏銷售量、價格、廣告費用及同期的市場均價。
（1）分析牙膏銷售量與價格、廣告投入之間的關系，建立數學模型；
（2）估計所建立數學模型的參數，進行統計分析；
（3）利用擬合模型，預測在不同價格和廣告費用下的牙膏銷售量。
6.2 Python 程序
# LinearRegression_v4.py
# v4.0: 分析和結果的可視化
# 日期：2021-05-08
# Copyright 2021 YouCans, XUPT
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt
import seaborn as sns
# 主程序
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
def main():
# 讀取數據文件
readPath = "../data/toothpaste.csv" # 數據文件的地址和文件名
dfOpenFile = pd.read_csv(readPath, header=0, sep=",") # 間隔符為逗號，首行為標題行
# 准備建模數據：分析因變量 Y(sales) 與 自變量 x1~x4 的關系
dfData = dfOpenFile.dropna() # 刪除含有缺失值的數據
sns.set_style('dark')
# 數據探索：分布特征
fig1, axes = plt.subplots(2, 2, figsize=(10, 8)) # 創建一個 2行 2列的畫布
sns.distplot(dfData['price'], bins=10, ax=axes[0,0]) # axes[0,1] 左上圖
sns.distplot(dfData['average'], bins=10, ax=axes[0,1]) # axes[0,1] 右上圖
sns.distplot(dfData['advertise'], bins=10, ax=axes[1,0]) # axes[1,0] 左下圖
sns.distplot(dfData['difference'], bins=10, ax=axes[1,1]) # axes[1,1] 右下圖
plt.show()
# 數據探索：相關性
fig2, axes = plt.subplots(2, 2, figsize=(10, 8)) # 創建一個 2行 2列的畫布
sns.regplot(x=dfData['price'], y=dfData['sales'], ax=axes[0,0])
sns.regplot(x=dfData['average'], y=dfData['sales'], ax=axes[0,1])
sns.regplot(x=dfData['advertise'], y=dfData['sales'], ax=axes[1,0])
sns.regplot(x=dfData['difference'], y=dfData['sales'], ax=axes[1,1])
plt.show()
# 數據探索：考察自變量平方項的相關性
fig3, axes = plt.subplots(1, 2, figsize=(10, 4)) # 創建一個 2行 2列的畫布
sns.regplot(x=dfData['advertise'], y=dfData['sales'], order=2, ax=axes[0]) # order=2, 按 y=b*x**2 回歸
sns.regplot(x=dfData['difference'], y=dfData['sales'], order=2, ax=axes[1]) # YouCans, XUPT
plt.show()
# 線性回歸：分析因變量 Y(sales) 與 自變量 X1(Price diffrence)、X2(Advertise) 的關系
y = dfData['sales'] # 根據因變量列名 list，建立 因變量數據集
x0 = np.ones(dfData.shape[0]) # 截距列 x0=[1,...1]
x1 = dfData['difference'] # 價格差，x4 = x1 - x2
x2 = dfData['advertise'] # 廣告費
x3 = dfData['price'] # 銷售價格
x4 = dfData['average'] # 市場均價
x5 = x2**2 # 廣告費的二次元
x6 = x1 * x2 # 考察兩個變量的相互作用
# Model 1：Y = b0 + b1*X1 + b2*X2 + e
# # 線性回歸：分析因變量 Y(sales) 與 自變量 X1(Price diffrence)、X2(Advertise) 的關系
X = np.column_stack((x0,x1,x2)) # [x0,x1,x2]
Model1 = sm.OLS(y, X) # 建立 OLS 模型: Y = b0 + b1*X1 + b2*X2 + e
result1 = Model1.fit() # 返回模型擬合結果
yFit1 = result1.fittedvalues # 模型擬合的 y 值
prstd, ivLow, ivUp = wls_prediction_std(result1) # 返回標准偏差和置信區間
print(result1.summary()) # 輸出回歸分析的摘要
print("\nModel1: Y = b0 + b1*X + b2*X2")
print('Parameters: ', result1.params) # 輸出：擬合模型的系數
# # Model 2：Y = b0 + b1*X1 + b2*X2 + b3*X3 + b4*X4 + e
# 線性回歸：分析因變量 Y(sales) 與 自變量 X1~X4 的關系
X = np.column_stack((x0,x1,x2,x3,x4)) #[x0,x1,x2,...,x4]
Model2 = sm.OLS(y, X) # 建立 OLS 模型: Y = b0 + b1*X1 + b2*X2 + b3*X3 + e
result2 = Model2.fit() # 返回模型擬合結果
yFit2 = result2.fittedvalues # 模型擬合的 y 值
prstd, ivLow, ivUp = wls_prediction_std(result2) # 返回標准偏差和置信區間
print(result2.summary()) # 輸出回歸分析的摘要
print("\nModel2: Y = b0 + b1*X + ... + b4*X4")
print('Parameters: ', result2.params) # 輸出：擬合模型的系數
# # Model 3：Y = b0 + b1*X1 + b2*X2 + b3*X2**2 + e
# # 線性回歸：分析因變量 Y(sales) 與 自變量 X1、X2 及 X2平方（X5）的關系
X = np.column_stack((x0,x1,x2,x5)) # [x0,x1,x2,x2**2]
Model3 = sm.OLS(y, X) # 建立 OLS 模型: Y = b0 + b1*X1 + b2*X2 + b3*X2**2 + e
result3 = Model3.fit() # 返回模型擬合結果
yFit3 = result3.fittedvalues # 模型擬合的 y 值
prstd, ivLow, ivUp = wls_prediction_std(result3) # 返回標准偏差和置信區間
print(result3.summary()) # 輸出回歸分析的摘要
print("\nModel3: Y = b0 + b1*X1 + b2*X2 + b3*X2**2")
print('Parameters: ', result3.params) # 輸出：擬合模型的系數
# 擬合結果繪圖
fig, ax = plt.subplots(figsize=(8,6)) # YouCans, XUPT
ax.plot(range(len(y)), y, 'b-.', label='Sample') # 樣本數據
ax.plot(range(len(y)), yFit3, 'r-', label='Fitting') # 擬合數據
# ax.plot(range(len(y)), yFit2, 'm--', label='fitting') # 擬合數據
ax.plot(range(len(y)), ivUp, '--',color='pink',label="ConfR") # 95% 置信區間 上限
ax.plot(range(len(y)), ivLow, '--',color='pink') # 95% 置信區間 下限
ax.legend(loc='best') # 顯示圖例
plt.title('Regression analysis with sales of toothpaste')
plt.xlabel('period')
plt.ylabel('sales')
plt.show()
return
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
if __name__ == '__main__':
main()
6.3 程序運行結果：
OLS Regression Results
==============================================================================
Dep. Variable: sales R-squared: 0.886
Model: OLS Adj. R-squared: 0.878
Method: Least Squares F-statistic: 105.0
Date: Sat, 08 May 2021 Prob (F-statistic): 1.84e-13
Time: 22:18:04 Log-Likelihood: 2.0347
No. Observations: 30 AIC: 1.931
Df Residuals: 27 BIC: 6.134
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.4075 0.722 6.102 0.000 2.925 5.890
x1 1.5883 0.299 5.304 0.000 0.974 2.203
x2 0.5635 0.119 4.733 0.000 0.319 0.808
==============================================================================
Omnibus: 1.445 Durbin-Watson: 1.627
Prob(Omnibus): 0.486 Jarque-Bera (JB): 0.487
Skew: 0.195 Prob(JB): 0.784
Kurtosis: 3.486 Cond. No. 115.
==============================================================================
Model1: Y = b0 + b1*X + b2*X2
Parameters:
const 4.407493
x1 1.588286
x2 0.563482
OLS Regression Results
==============================================================================
Dep. Variable: sales R-squared: 0.895
Model: OLS Adj. R-squared: 0.883
Method: Least Squares F-statistic: 74.20
Date: Sat, 08 May 2021 Prob (F-statistic): 7.12e-13
Time: 22:18:04 Log-Likelihood: 3.3225
No. Observations: 30 AIC: 1.355
Df Residuals: 26 BIC: 6.960
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 8.0368 2.480 3.241 0.003 2.940 13.134
x1 1.3832 0.288 4.798 0.000 0.791 1.976
x2 0.4927 0.125 3.938 0.001 0.236 0.750
x3 -1.1184 0.398 -2.811 0.009 -1.936 -0.300
x4 0.2648 0.199 1.332 0.195 -0.144 0.674
==============================================================================
Omnibus: 0.141 Durbin-Watson: 1.762
Prob(Omnibus): 0.932 Jarque-Bera (JB): 0.030
Skew: 0.052 Prob(JB): 0.985
Kurtosis: 2.885 Cond. No. 2.68e+16
==============================================================================
Model2: Y = b0 + b1*X + ... + b4*X4
Parameters:
const 8.036813
x1 1.383207
x2 0.492728
x3 -1.118418
x4 0.264789
OLS Regression Results
==============================================================================
Dep. Variable: sales R-squared: 0.905
Model: OLS Adj. R-squared: 0.894
Method: Least Squares F-statistic: 82.94
Date: Sat, 08 May 2021 Prob (F-statistic): 1.94e-13
Time: 22:18:04 Log-Likelihood: 4.8260
No. Observations: 30 AIC: -1.652
Df Residuals: 26 BIC: 3.953
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 17.3244 5.641 3.071 0.005 5.728 28.921
x1 1.3070 0.304 4.305 0.000 0.683 1.931
x2 -3.6956 1.850 -1.997 0.056 -7.499 0.108
x3 0.3486 0.151 2.306 0.029 0.038 0.659
==============================================================================
Omnibus: 0.631 Durbin-Watson: 1.619
Prob(Omnibus): 0.729 Jarque-Bera (JB): 0.716
Skew: 0.203 Prob(JB): 0.699
Kurtosis: 2.362 Cond. No. 6.33e+03
==============================================================================
Model3: Y = b0 + b1*X1 + b2*X2 + b3*X2**2
Parameters:
const 17.324369
x1 1.306989
x2 -3.695587
x3 0.348612
版權說明：
YouCans 原創作品
Copyright 2021 YouCans, XUPT
Crated：2021-05-08
歡迎關注 Youcans 原創系列，每周更新數模筆記
Python數模筆記-PuLP庫（1）線性規划入門
Python數模筆記-PuLP庫（2）線性規划進階
Python數模筆記-PuLP庫（3）線性規划實例
Python數模筆記-Scipy庫（1）線性規划問題
Python數模筆記-StatsModels 統計回歸（1）簡介
Python數模筆記-StatsModels 統計回歸（2）線性回歸
Python數模筆記-StatsModels 統計回歸（3）模型數據的准備
Python數模筆記-StatsModels 統計回歸（4）可視化
Python數模筆記-Sklearn （1）介紹
Python數模筆記-Sklearn （2）聚類分析
Python數模筆記-Sklearn （3）主成分分析
Python數模筆記-Sklearn （4）線性回歸
Python數模筆記-Sklearn （5）支持向量機
Python數模筆記-NetworkX（1）圖的操作
Python數模筆記-NetworkX（2）最短路徑
Python數模筆記-NetworkX（3）條件最短路徑
Python數模筆記-模擬退火算法（1）多變量函數優化
Python數模筆記-模擬退火算法（2）約束條件的處理
Python數模筆記-模擬退火算法（3）整數規划問題
Python數模筆記-模擬退火算法（4）旅行商問題