利用 python 進行線性迴歸
線性迴歸也被稱為最小二乘法迴歸（Linear Regression, also called Ordinary Least-Squares (OLS) Regression）。它的數學模型是這樣的：
y = a+ b* x＋e
其中，a 被稱為常數項或截距；b 被稱為模型的迴歸係數或斜率；e 為誤差項。a 和 b 是模型的引數。
當然，模型的引數只能從樣本資料中估計出來：
y'= a' + b'* x
我們的目標是選擇合適的引數，讓這一線性模型最好地擬合觀測值。擬合程度越高，模型越好。
那麼，接下來的問題就是，我們如何判斷擬合的質量呢？
這一線性模型可以用二維平面上的一條直線來表示，被稱為迴歸線。
模型的擬合程度越高，也即意味著樣本點圍繞回歸線越緊密。
如何計算樣本點與迴歸線之間的緊密程度呢？
高斯和勒讓德找到的方法是：被選擇的引數，應該使算出來的迴歸線與觀測值之差的平房和最小。用函式表示為
這被稱為最小二乘法。最小二乘法的原理是這樣的：當預測值和實際值距離的平方和最小時，就選定模型中的兩個引數（a 和 b）。這一模型並不一定反映解釋變數和反應變數真實的關係。但它的計算成本低；相比複雜模型更容易解釋。
模型估計出來後，我們要回答的問題是：
我們的模型擬合程度如何？或者說，這個模型對因變數的解釋力如何？（R2）
整個模型是否能顯著預測因變數的變化？（F 檢驗）
每個自變數是否能顯著預測因變數的變化？（t 檢驗）
首先回答第一個問題。為了評估模型的擬合程度如何，我們必須有一個可以比較的基線模型。
如果讓你預測一個人的體重是多少？在沒有任何額外資訊的情況下，你可能會用平均值來預測，儘管會存在一定誤差，但總比瞎猜好。
現在，如果你知道他的身高資訊，你的預測值肯定與平均值不一樣。額外資訊相比平均值更能準確地預測被預測的變數的能力，就代表模型的解釋力大小。
上圖中，SSA 代表由自變數 x 引起的 y 的離差平方和，即迴歸平方和，代表迴歸模型的解釋力；SSE 代表由隨機因素引起的 y 的離差平方和，即剩餘平方和，代表迴歸模型未能解釋的部分；SST 為總的離差平方和，即我們僅憑 y 的平均值去估計 y 時所產生的誤差。
用模型能夠解釋的變異除以總的變異就是模型的擬合程度：
R2=SSA/SST=1-SSE
R2（R 的平方）也被稱為決定係數或判定係數。
第二個問題，我們的模型是否顯著預測了 y 的變化？
假設 y 與 x 的線性關係不明顯，那麼 SSA 相對 SSE 佔有較大的比例的概率則越小。換句話說，在 y 與 x 無線性關係的前提下，SSA 相對 SSE 的佔比越高的概率是越小的，這會呈現一定的概率分佈。統計學家告訴我們它滿足 F 分佈，就像這樣：
如果 SSA 相對 SSE 佔比較大的情況出現了，比如根據 F 分佈，這個值出現的概率小於 5%。那麼，我們最好是拒絕 y 與 x 線性關係不顯著的原始假設，認為二者存在顯著的線性關係較為合適。
第三個問題，每個自變數是否能顯著預測因變數的變化？換句話說，迴歸係數是否顯著？
迴歸係數的顯著性檢驗是圍繞回歸係數的抽樣分佈（t 分佈）來進行的，推斷過程類似於整個模型的檢驗過程，不贅言。
實際上，對於只有一個自變數的一元線性模型，模型的顯著性檢驗和迴歸係數的檢驗是一致的，但對於多元線性模型來說，二者就不能等價了。
利用 statsmodels 進行最小二乘迴歸
＃匯入相應模組
In [1]: import numpy as np
In [2]: import pandas as pd
In [3]: import statsmodels.api as sm
＃將資料匯入 pandas 的 dataframe 物件，第一列（年份）作為行標籤
In [4]: df=pd.read_csv('/Users/xiangzhendong/Downloads/vincentarelbundock-Rdatasets-1218370/csv/datasets/longley.csv', index_col=0)
＃檢視頭部資料
In [5]: df.head()
Out[5]:
GNP.deflator GNP Unemployed Armed.Forces Population Year \1947 83.0 234.289 235.6 159.0 107.608 1947
1948 88.5 259.426 232.5 145.6 108.632 1948
1949 88.2 258.054 368.2 161.6 109.773 1949
1950 89.5 284.599 335.1 165.0 110.929 1950
1951 96.2 328.975 209.9 309.9 112.075 1951
Employed
1947 60.323
1948 61.122
1949 60.171
1950 61.187
1951 63.221
＃設定預測變數和結果變數，用 GNP 預測 Employed
In [6]: y=df.Employed ＃結果變數
In [7]: X=df.GNP ＃預測變數
＃為模型增加常數項，即迴歸線在 y 軸上的截距
In [8]: X=sm.add_constant(X)
＃執行最小二乘迴歸，X 可以是 numpy array 或 pandas dataframe（行數等於資料點個數，列數為預測變數個數），y 可以是一維陣列（numpy array）或 pandas series
In [10]: est=sm.OLS(y,X)
使用 OLS 物件的 fit() 方法來進行模型擬合
In [11]: est=est.fit()
＃檢視模型擬合的結果
In [12]: est.summary()
Out[12]:
＃檢視最終模型的引數
In [13]: est.params
Out[13]:
const 51.843590
GNP 0.034752
dtype: float64
＃選擇 100 個從最小值到最大值平均分佈（equally spaced）的資料點
In [14]: X_prime=np.linspace(X.GNP.min(), X.GNP.max(),100)[:,np.newaxis]
In [15]: X_prime=sm.add_constant(X_prime)
＃計算預測值
In [16]: y_hat=est.predict(X_prime)
In [17]: plt.scatter(X.GNP, y, alpha=0.3) ＃畫出原始資料
＃分別給 x 軸和 y 軸命名
In [18]: plt.xlabel("Gross National Product")
In [19]: plt.ylabel("Total Employment")
In [20]: plt.plot(X_prime[:,1], y_hat, 'r', alpha=0.9) ＃添加回歸線，紅色
多元線性迴歸（預測變數不止一個）
我們用一條直線來描述一元線性模型中預測變數和結果變數的關係，而在多元迴歸中，我們將用一個多維（p）空間來擬合多個預測變數。下面表現了兩個預測變數的三維圖形：商品的銷量以及在電視和廣播兩種不同媒介的廣告預算。
數學模型是：
Sales = beta_0 + beta_1＊TV + beta_2＊Radio
圖中，白色的資料點是平面上的點，黑色的資料點事平面下的點。平面的顏色是由對應的商品銷量的高低決定的，高是紅色，低是藍色。
利用 statsmodels 進行多元線性迴歸
In [1]: import pandas as pd
In [2]: import numpy as np
In [3]: import statsmodels.api as sm
In [4]: df_adv=pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv',index_col=0)
In [6]: X=df_adv[['TV','Radio']]
In [7]: y=df_adv['Sales']
In [8]: df_adv.head()
Out[8]:
TV Radio Newspaper Sales
1 230.1 37.8 69.2 22.1
2 44.5 39.3 45.1 10.4
3 17.2 45.9 69.3 9.3
4 151.5 41.3 58.5 18.5
5 180.8 10.8 58.4 12.9
In [9]: X=sm.add_constant(X)
In [10]: est=sm.OLS(y,X).fit()
In [11]: est.summary()
Out[11]:
你也可以使用 statsmodels 的 formula 模組來建立多元迴歸模型
In [12]: import statsmodels.formula.api as smf
In [13]: est=smf.ols(formula='Sales ~ TV + Radio',data=df_adv).fit()
處理分類變數
性別或地域都屬於分類變數。
In [15]: df= pd.read_csv('http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/SAheart.data', index_col=0)
In [16]: X=df.copy()
利用 dataframe 的 pop 方法將 chd 列單獨提取出來
In [17]: y=X.pop('chd')
In [18]: df.head()
Out[18]:
sbp tobacco ldl adiposity famhist typea obesity alcohol \
row.names
1 160 12.00 5.73 23.11 Present 49 25.30 97.20
2 144 0.01 4.41 28.61 Absent 55 28.87 2.06
3 118 0.08 3.48 32.28 Present 52 29.14 3.81
4 170 7.50 6.41 38.03 Present 51 31.99 24.26
5 134 13.60 3.50 27.78 Present 60 25.99 57.34
age chd
row.names
1 52 1
2 63 1
3 46 0
4 58 1
5 49 1
In [19]: y.groupby(X.famhist).mean()
Out[19]:
famhist
Absent 0.237037
Present 0.500000
Name: chd, dtype: float64
In [20]: import statsmodels.formula.api as smf
In [21]: df['famhist_ord']=pd.Categorical(df.famhist).labels
In [22]: est=smf.ols(formula="chd ~ famhist_ord", data=df).fit()
分類變數的編碼方式有許多，其中一種編碼方式是虛擬變數編碼（dummy－encoding），就是把一個 k 個水平的分類變數編碼成 k－1 個二分變數。在 statsmodels 中使用 C 函式實現。
In [24]: est=smf.ols(formula="chd ~ C(famhist)", data=df).fit()
In [26]: est.summary()
Out[26]:
處理互動作用
隨著教育年限（education）的增長，薪酬 (wage) 會增加嗎？這種影響對男性和女性而言是一樣的嗎？
這裡的問題就涉及性別與教育年限的互動作用。
換言之，教育年限對薪酬的影響是男女有別的。
＃匯入相關模組
In [1]: import pandas as pd
In [2]: import numpy as np
In [4]: import statsmodels.api as sm
＃匯入資料，存入 dataframe 物件
In [5]: df=pd.read_csv('/Users/xiangzhendong/Downloads/pydatafromweb/wages.csv')
In [6]: df[['Wage','Education','Sex']].tail()
Out[6]:
Wage Education Sex
529 11.36 18 0
530 6.10 12 1
531 23.25 17 1
532 19.88 12 0
533 15.38 16 0
由於性別是一個二分變數，我們可以繪製兩條迴歸線，一條是 sex＝0（男性），一條是 sex＝1（女性）
＃繪製散點圖
In [7]: plt.scatter(df.Education,df.Wage, alpha=0.3)
In [9]: plt.xlabel('education')
In [10]: plt.ylabel('wage')
＃linspace 的作用是生成從最小到最大的均勻分佈的 n 個數
In [17]: education_linspace=np.linspace(df.Education.min(), df.Education.max(),100)
In [12]: import statsmodels.formula.api as smf
In [13]: est=smf.ols(formula='Wage ~ Education + Sex', data=df).fit()
In [18]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]0, 'r')
In [19]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]1, 'g')
作者：wyrover
連結：https://www.jianshu.com/p/e55a8c9e4b56
來源：簡書
著作權歸作者所有。商業轉載請聯絡作者獲得授權，非商業轉載請註明出處。
以上兩條線是平行的。這是因為分類變數隻影響迴歸線的截距，不影響斜率。
接下來我們可以為迴歸模型增加互動項來探索互動效應。也就是說，對於兩個類別，迴歸線的斜率是不一樣的。
In [32]: plt.scatter(df.Education,df.Wage, alpha=0.3)
In [33]: plt.xlabel('education')
In [34]: plt.ylabel('wage')
＃使用＊代表我們的迴歸模型中除了互動效應，也包括兩個變數的主效應；如果只想看互動效應，可以用：代替，但通常不會只看互動效應
In [35]: est=smf.ols(formula='Wage ~ Sex*Education', data=df).fit()
In [36]: plt.plot(education_linspace, est.params[0]+est.params[1]0+est.params[2]education_linspace+est.params[3]0education_linspace, 'r')
In [37]: plt.plot(education_linspace, est.params[0]+est.params[1]1+est.params[2]education_linspace+est.params[3]1education_linspace, 'g')
參考資料：
DataRobot | Ordinary Least Squares in Python
DataRoboe | Multiple Regression using Statsmodels
AnalyticsVidhya | 7 Types of Regression Techniques you should know!
維基百科 | 最小二乘法
作者：wyrover
連結：https://www.jianshu.com/p/e55a8c9e4b56
來源：簡書
著作權歸作者所有。商業轉載請聯絡作者獲得授權，非商業轉載請註明出處。