線性回歸
https://sanwen8.cn/p/3cbCi2d.html
理解什麽是線性回歸
線性回歸也被稱為最小二乘法回歸（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-121****70/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‘)
以上兩條線是平行的。這是因為分類變量只影響回歸線的截距，不影響斜率。
接下來我們可以為回歸模型增加交互項來探索交互效應。也就是說，對於兩個類別，回歸線的斜率是不一樣的。
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]*0*education_linspace, ‘r‘)
In [37]: plt.plot(education_linspace, est.params[0]+est.params[1]*1+est.params[2]*education_linspace+est.params[3]*1*education_linspace, ‘g‘)
參考資料：
-
DataRobot | Ordinary Least Squares in Python
-
DataRoboe | Multiple Regression using Statsmodels
-
AnalyticsVidhya | 7 Types of Regression Techniques you should know!
-
維基百科 | 最小二乘法
線性回歸