1. 程式人生 > >R建模之迴歸(一)

R建模之迴歸(一)

3種常見的迴歸模型:

線性迴歸(預測連續型變數比如嬰兒出生體重),邏輯迴歸(預測二元變數比如過低出生體重與正常出生體重),泊松分佈(計數比如每年或每個國家過低出生體重嬰兒人數)

我們以gamlss.data包提供的usair資料集進行研究,US空氣汙染資料集。我們希望預測根據城市面積(以人口規模/千人為統計依據)估計的空氣汙染程度(這裡也就是資料集中的x3),空氣汙染以每立方米空氣中二氧化硫的含量(毫克)為指標(也就是資料集中的y)

首先檢視一下資料集的資料結構:


library(gamlss.data)           #匯入資料集所在的包
data(usair)
mode1.0<-lm(y~x3,data=usair)   #線性擬合函式
summary(mode1.0)               #獲取描述性統計量

上述程式碼執行結果如下:


擬合的方程x3係數為0.02,常數項係數為:17.87,將回歸結果視覺化:

plot(y~x3,data=usair,cex.lab=1.5)  #y和x3的點圖,並設定座標軸標籤的縮放倍數
abline(mode1.0,col="red",lwd=2.5)  #新增擬合結果直線,並設定顏色和線條寬度
legend('bottomright',legend='y~x3',lty=1,col='red',lwd=2.5,title='Regression line')  #新增圖例,第一個引數表式圖例的位置,第二個表示內容,第三個線的型別,第四個顏色第五個寬度,第六個圖例的標題


圖中直線很好地擬合了資料。我們使用最小二乘法來實現最佳擬合,所以該模型也被稱為普通最小二乘迴歸(OLS)

那麼上面的直線是如何得到的呢?是通過使殘差和最小來找到最佳的擬合曲線,殘差指的是觀測值(散點圖中的原始資料點)與預測理論值(迴歸直線上對應的值)的差值:

usair$prediction<-predict(mode1.0)    #預測結果儲存在prediction變數裡面
usair$residual<-resid(mode1.0)        #殘差
plot(y~x3,data=usair,cex.lab=1.5)     #畫原始資料影象
abline(mode1.0,col='red',lwd=2.5)       #擬合曲線
segments(usair$x3,usair$y,usair$x3,usair$prediction,col='blue',lty=2)       #對線段作圖,
legend('bottomright',legend=c('y~x3','residuals'),lty=c(1,2),col=c('red','blue'),lwd=2.5,title='Regression line')  #新增圖例

以上就是對於單變數的擬合。另一方面,如果我們希望將人口數與工業化的影響區分考慮,建立更復雜的模型,需要引入變數x2,該變數為工人數大於20的工廠個數。可以有兩種方法對模型進行更新:

mode1.1<-lm(y~x3+x2,data=usair)
summary(mode1.1)
mode1.1<-update(mode1.0,.~.,+x2)
summary(mode1.1)


從圖中可以看到,x3的係數是-0.05661,而第一個模型中x3的係數為0.02,這說明,當引入x2這一變數後,y和x3的關係變成了負相關,這表示人口數每上升1000,空氣中so2濃度下降0.05661單位,該結論具有統計顯著效應。當我們固定工業化水平,人口規模和空氣汙染之間的關聯就變為負相關,這意味著城市工業化水平不變,城市規模增加,則汙染物擴散的面積也增加了。由此,可以得到結論,變數x2可能是一個混雜因子,它歪曲了y和x3之間的關聯關係。

基於上述模型,我們可以預測變數的響應值:比如:預測一個擁有300 000居民和160個工廠,平均每個工廠人數都超過20的城市二氧化硫濃度:


可以看出效果還是不錯的。

如果有兩個預測變數,那麼迴歸直線就是一個平面了,使用scatterplot3d展示:


library(scatterplot3d)
plot3d<-scatterplot3d(usair$x3,usair$x2,usair$y,pch=19,type='h',       #繪製3d圖形
     highlight.3d=TRUE,main='3-D Scatterplot')
 plot3d$plane3d(mode1.1,lty='solid',col='red')

一般是看二位投影比較好:

par(mfrow=c(1,2))                                  #畫一行兩列的圖形
 plot(y~x3,data=usair,main='2D projection for x3')
 abline(mode1.1,col='red',lwd=2.5)                   #添加回歸直線
plot(y~x2,data=usair,main='2D projection for x2' )
abline(lm(y~x2+x3,data=usair),col='red',lwd=2.5)      
#添加回歸直線

5.3 模型假定 

以下:(完全引用)

  採用標準預測技術的線性迴歸模型能夠對輸出變數、預測變數以及兩者之間的相關性進行假設:
  1)Y是一連續變數(非二元變數,定類變數或定序變數)
  2)錯誤(殘差)具有統計無關性;
  3)在Y和每個X之間具有隨機的線性關聯
  4)固定X,Y服從正態分佈

  5)不論X的固定值是多少,Y擁有相同的方差

假定2)的特例發生在趨勢分析中,如果我們使用時間作為預測變數,由於連續年份是相互依賴的,因此相互之間的錯誤不會獨立無關。

 假定3)的特例是變數之間相互關係不是完全線性相關,而是背離了趨勢直線。假定4和假定5要求Y的條件分佈服從正態分佈,假定5也被稱為同方差假定。如果違背了該假定,線性迴歸模型則存在異方差性。

寫了這麼多,我自己也不是很理解,但是下面的程式出來就能明白一些了:

library(Hmisc)
library(ggplot2)
library(gridExtra)
set.seed(7) 
x <-sort(rnorm(1000,10,100))[26:975]
y <-x*500+rnorm(950,5000,20000)
df <-data.frame(x=x,y=y,cuts=factor(cut2(x,g=5)),resid=resid(lm(y~x)))
scatterP1<-ggplot(df,aes(x=x,y=y))+
           geom_point(aes(colour=cuts,fill=cuts),shape=1,show.legend=FALSE)+
           geom_smooth(method=lm,level=0.99)
plot_left<-ggplot(df,aes(x=y,fill=cuts))+
           geom_density(alpha=.5)+coord_flip()+scale_y_reverse()
plot_right<- ggplot(data=df,aes(x=resid,fill=cuts))+
           geom_density(alpha=.5)+coord_flip()
grid.arrange(plot_left,scatterP1,plot_right,ncol=3,nrow=1,widths=c(1,3,1))

下面介紹去掉孤立點的方法,去掉孤立點後有可能使得假定成立,可以使用gvlma包來快速驗證上面五個模型假定是否都能滿足,

library(gvlma)
gvlma(mode1.1)

五個假定只有三個能滿足。當我們把第31個觀測點去掉,再對R方進行分析:


從結果可以看出去掉31個觀測值,擬合效果變得更好。

下面簡要介紹評價迴歸線的擬合效果(AIC)

      儘管我們知道結果所得的趨勢直線在所有可能的線性直線中擬合效果最佳,但是我們並不知道其對於實際資料的擬合效果。這個時候可以通過驗證零假設來求得迴歸引數的顯著性,零假設假定迴歸引數等於零。

F檢驗適用於迴歸引數確實為零的情況。p檢驗適用於一般情況的迴歸顯著性檢驗,p值小於0.05即說明迴歸直線具有顯著性。

除了顯著的F值,殘差是對擬合誤差的度量。R平方值將以上這些值統一為一個度量。

R平方值是指:響應變數的變異中有多少百分比是有迴歸模型解釋的;數學表示是預測Y值得方差除以觀測Y值的方差。

mode1.2明顯比mod1.0R平方大,同種第一個是R平方值,第二個是R平方調整值,但是R平方調整值不能作為非巢狀模型篩選的標準。如果使用者的模型為非巢狀模型,可以考慮使用赤池資訊量準則(AIC),優先考慮AIC值小的那一個,如果兩個模型的差小於2,這兩個模型差別不大:

summary(mode1.3 <-update(mode1.2,.~.-x2+x1))$coefficients
summary(mode1.4 <-update(mode1.2,.~.-x3+x1))$coefficients
AIC(mode1.3,mode1.4)       #AIC準則,赤池資訊量準則