前言:
線性回歸模型屬於經典的統計學模型,該模型的應用場景是根據已知的變量(即自變量)來預測某個連續的數值變量(即因變量)。例如餐廳根據媒體的營業數據(包括菜譜價格、就餐人數、預訂人數、特價菜折扣等)預測就餐規模或營業額;網站根據訪問的歷史數據(包括新用戶的注冊量、老用戶的活躍度、網站內容的更新頻率等)預測用戶的支付轉化率;醫院根據患者的病歷數據(如體檢指標、藥物復用情況、平時的飲食習慣等)預測某種疾病發生的概率。
由於針對具體操作相關文檔太多,所以本文內容涉及具體操作較少,主要是講方法。
本文內所用到的包:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import scipy.stats as stats
import statsmodels.api as sm
from scipy.stats import chi2_contingency
以經典的剎車距離為例:
ccpp=pd.read_csv('cars.csv')
sns.set(font=getChineseFont(8).get_name())
sns.lmplot(x='speed',y='dist',data=ccpp,
legend_out=False,#將圖例呈現在圖框內
truncate=True#根據實際的數據范圍,對擬合線做截斷操作
)
plt.show()
從散點圖來看,自變量speed與因變量dist之間存在明顯的正相關,即剎車速度越大,剎車距離越長。圖內的陰影部分為擬合線95%的置信區間,每個散點都是盡可能地圍繞在擬合線附近。
通過ols函數求得回歸模型的參數解’y~s’表示簡單線性回歸模型
fit=sm.formula.ols('dist~speed',data=ccpp).fit()
print(fit.params)#Intercept:-17.579095,speed:3.932409
因此,關於剎車距離的簡單線性模型為:
dist=-17.579095+3.932409speed
也就是說,剎車速度每提升一個單位,將促使剎車距離增加3.93個單位。
以利潤表Profit為例,研究影響利潤的因素
表結構如下:
profit=pd.read_csv(‘Profit.csv’,sep=“,”)
fit=sm.formula.ols('Profit~RD_Spend+Administration+Marketing_Spend',data=profit).fit()
print(fit.params)#Intercept:50122.192990,RD_Spend:0.805715,Administration:-0.026816,Marketing_Spend:0.027228
不考慮模型的顯著性和回歸系數的顯著性,得到的回歸模型可以表示為:
Profit=50122.192990+0.805715RD_Spend-0.026816Administration+0.027228Marketing_Spend
但是,在實際應用中,單純的利用ols函數將偏回歸系數計算出來,並構造一個多元線性回歸模型,得到的結果往往不是理想的,這時候需要借助於統計學中的F檢驗法和t檢驗法,完成模型的顯著性檢驗和回歸系數的顯著性檢驗。
步驟如下:
假設認為模型的所有偏回歸系數全為0(即認為沒有一個自變量可以構成因變量的線性組合)
通常在實際的應用中將概率值p與0.05做比較,小於0.05則拒絕原假設,否則接受原假設
模型通過了顯著性檢驗,只能說明模型關於因變量的線性組合是合理的,但不能說明每個自變量對因變量都具有顯著意義,所以還要對模型的回歸系數做顯著性檢驗。
只有當回歸系數通過了t檢驗,才可以認為模型的系數是顯著的
t檢驗的出發點就是驗證每一個自變量是否能夠成為影響因變量的重要因素
t檢驗的原假設是假定第j變量的回歸系數為0,即認為該變量不是因變量的影響因素
t統計量大於理論的t分布值,則拒絕原假設,否則接受原假設;也可以根據概率值P判斷是否需要拒絕原假設
#返回模型概覽
print(fit.summary())
由圖可知:F-statistic:296.0,Prob (F-statistic):4.53e-30,F統計量值為296.0,對應的概率值P遠遠小於0.05,說明應該拒絕原假設,認為模型是顯著的。
在各自變量的t統計中,Administration和Marketing_Spend變量所對應的概率值p大於0.05,說明不能拒絕原假設,認為該變量是不顯著的,無法認定其實影響Profit的重要因素
根據返回的fit模型的概覽信息,由於Administration和Marketing_Spend變量的t檢驗結果是不顯著的,故可以探索其余因變量Profit之間的散點關系,如果確實沒有線性關系,可將其從模型中剔除。
sns.lmplot(x='Administration',y='Profit',data=profit,
legend_out=False,#將圖例呈現在圖框內
fit_reg=False#不顯示擬合曲線
)
sns.lmplot(x='Marketing_Spend',y='Profit',data=profit,
legend_out=False,#將圖例呈現在圖框內
fit_reg=False#不顯示擬合曲線
)
plt.show()
圖中自變量Administration、Marketing_Spend與因變量Profit沒有呈現明顯的線性關系,故可以認為兩者不存在互相依賴關系
#將Administration、Marketing_Spend變量從模型中剔除
fit2 = sm.formula.ols('Profit~RD_Spend',data=profit).fit()
print(fit2.params)#Intercept:49032.899141,RD_Spend :0.854291
print(fit2.summary())#Prob (F-statistic):3.50e-32;P>|t|:0.000
對模型fit重新調整後,得到的新模型fit2仍然通過了顯著性檢驗,而且每個自變量所對應的系數也是通過顯著性檢驗的。
最終得到的模為:
Profit=49032.899141+0.854291RD_Spend
該回歸模型中的系數解釋為:在其他條件不變的情況下RD_Spend每增加一個單位,將使Profit增加0.854291個單位。
繼續使用上面的數據。
#異常值檢驗
outliers=fit2.get_influence()
#高槓桿值點(帽子矩陣)
leverage=outliers.hat_matrix_diag
#diffits值
dffits=outliers.dffits[0]
#學生化殘差
resid_stu=outliers.resid_studentized_external
#cook距離
cook=outliers.cooks_distance[0]
#合並以上4種異常值檢驗的統計量值
concat_result=pd.concat([pd.Series(leverage,name=‘leverage’),pd.Series(dffits,name=‘diffits’),
pd.Series(resid_stu,name=‘resid_stu’),pd.Series(cook,name=‘cook’)],axis=1)
#將上面的concat_result結果與profit數據集合並
raw_outliers=pd.concat([profit,concat_result],axis=1)
前5行數據集打印結果如下:
簡單起見,這裡使用學生化殘差,當學生化殘差大於2時,即認為對應的數據點為異常值
#計算異常值數量的比例
outliers_ratio=sum(np.where(np.abs(raw_outliers.resid_stu)>2,1,0))/raw_outliers.shape[0]
print(outliers_ratio)#0.04
結果顯示,通過學生化殘差識別出了異常值,並且異常值比例為4%。由於異常值比例非常小,故可以考慮將其直接從數據集中刪除,由此繼續建模將會得到更加穩健且合理的模型
#通過篩選的方法,將異常點排除
none_outliers=raw_outliers.loc[np.abs(raw_outliers.resid_stu)<=2,]
#應用無異常值的數據集重新建模
fit3=sm.formula.ols('Profit~RD_Spend',data=profit).fit()
#返回模型的概覽信息
print(fit3.params)#Intercept:49032.899141,RD_Spend:0.854291
print(fit3.summary())
排除異常點之後得到的模型fit3,不管是模型的顯著性檢驗還是系數的顯著性檢驗,各自的概率P值均小於0.05,說明他們均通過顯著性檢驗。
模型fit3為:Profit=49032.899141+0.854291RD_Spend
從fit3模型來看RD_Spend(研發成本)每增加一個單位的成本,所帶來的Profit(收益)提升要明顯比其他的高,所以將更多的成本花費在研發上是個不錯的選擇。
以原始數據profit為例,根據fit3模型重新預測各成本下的收益預測值:
pred=fit3.predict(profit[['RD_Spend']])
#對於實際值與預測值的比較
df=pd.concat([pd.Series(profit.Profit/100,name='real'),pd.Series(pred/100,name='prediction')],axis=1)
df['誤差絕對值']=np.abs((df['real']-df['prediction'])/100)
print(df.head(10))
從結果上看有的預測值比較接近實際值,有的預測測偏離實際值較遠,但從總體上來說,預測值與實際值之間的差異並不是特別大。
以經典的泰坦尼克號數據為例:
titanic=pd.read_csv('Titanic_all.csv',sep=",")
print(titanic.dtypes)
12個變量中涉及5個離散型變量和7個數值型變量。根據知己情況可知,船艙等級Pclass應為離散型變量(3等、2等、1等,並非數值)
查看各變量的缺失比例
print(titanic.isnull().sum(axis=0)/titanic.shape[0])
Cabin缺失比例超過77%,Embarked缺失比例僅為0.22%,二者數據刪除。
在這裡我們需要利用年齡的非缺失數據構造多元線性回歸模型,進而預測缺失比例為19.87%的乘客年齡
如需基於其余變量來預測年齡變量Age,至少有5個變量與年齡無關(乘客編號、姓名、票號信息、座位號信息和登船地點),刪除。
1.刪除無意義的變量
titanic.drop(labels=['PassengerId','Name','Ticket','Cabin','Embarked'],axis=1,inplace=True)
2.啞變量轉換
titanic.Pclass=titanic.Pclass.astype(str)#Pclass變量轉換為字符型變量
#將Pclass和Sex變量作啞變量處理
dummies=pd.get_dummies(data=titanic[['Pclass','Sex']])
#將titanic數據集與dummies數據集進行合並
titanic_new=pd.concat([titanic,dummies],axis=1)
#根據啞變量原則,需要從新生成的0和1變量中提出兩個變量作為參照變量(以性別中的男性為參照變量,以船艙等幾種的Pclass_3作為參照變量)
titanic_new.drop(labels=['Pclass','Sex','Pclass_3','Sex_male'],axis=1,inplace=True)
3.將數據拆分為兩部分
在清洗後的數據集titanic_new中,僅有Age變量存在缺失值
為預測該變量的缺失值,需要將數據集按照年齡是否缺失拆分為兩個數據子集,分別用於線性回歸模型的構造和基於構造好的數據集對其進行預測
missing=titanic_new.loc[titanic.Age.isnull(),]
no_missing=titanic_new.loc[~titanic.Age.isnull(),]
4.構建多元線性回歸模型
#提取出所有的自變量名稱
predictors=no_missing.columns[~no_missing.columns.isin(['Age'])]
#構造多元線性回歸模型的"類"
lm_Class=sm.OLS(endog=no_missing.Age,#指定模型中的因變量
exog=no_missing[predictors]#指定模型中的自變量
)
#基於"類",對模型進行擬合
lm_model=lm_Class.fit()
print(lm_model.summary())
根據F檢驗的結果可知模型是顯著的,但是從t檢驗的結果來看,僅有船艙等級Pclass和性別Sex是通過顯著性檢驗的。
繪制其他變量與年齡之間的散點圖
sns.pairplot(no_missing[['Survived','SibSp','Parch','Fare','Age']])
唯有SibSp與Age之間存在一定趨勢性,將出SibSp之外其他變量從模型中剔除。
基於散點圖結果重新構造多元線性回歸模型
predictors2=['SibSp','Pclass_1','Pclass_2','Sex_female']
lm_Class2=sm.OLS(endog=no_missing.Age,#指定模型中的因變量
exog=no_missing[predictors2]#指定模型中的自變量
)
lm_model2=lm_Class2.fit()
print(lm_model2.summary())
在新的模型中,F檢驗的統計量和t檢驗的統計量所對應的概率P值均小於0.05,說明他們通過了顯著性檢驗。
最後再利用可視化的QQ圖,驗證模型的殘差是否服從正態分布的假設
sns.set(font=getChineseFont(8).get_name())
pp_qq_plot=sm.ProbPlot(lm_model2.resid)
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q圖')
plt.show()
圖形中的散點基本上都圍繞在斜線附近,可以判斷模型lm_model2的殘差服從正態分布,最終證明了lm_model2模型是合理的。
5.未知年齡的預測
pred_Age=lm_model2.predict(exog=missing[predictors2])
#將年齡的預測值替換missing數據集中的Age變量
missing.loc[:,'Age']=pred_Age
#print(missing.head(5))
#將結果插補到原始數據中
先自我介紹一下,小編13年上師交大畢業,曾經在小公司待過,去過華為OPPO等大廠,18年進入阿裡,直到現在。深知大多數初中級java工程師,想要升技能,往往是需要自己摸索成長或是報班學習,但對於培訓機構動則近萬元的學費,著實壓力不小。自己不成體系的自學效率很低又漫長,而且容易碰到天花板技術停止不前。因此我收集了一份《java開發全套學習資料》送給大家,初衷也很簡單,就是希望幫助到想自學又不知道該從何學起的朋友,同時減輕大家的負擔。添加下方名片,即可獲取全套學習資料哦