《統(tǒng)計分析軟件：使用R與Python》 課件 第7章 R 語言回歸分析

7.1一元回歸模型7.2多元回歸模型7.3回歸模型的擬合優(yōu)度7.4回歸模型診斷7.5模型選擇7.6模型的預(yù)測第7章R語言回歸分析R語言回歸分析一元回歸模型lm函數(shù)，R語言公式回歸模型診斷殘差分析、異常點檢測、共線性檢測（條件數(shù)、VIF、FG檢驗）模型選擇模型的預(yù)測多元回歸模型系數(shù)置信區(qū)間、標(biāo)準(zhǔn)化回歸模型的擬合優(yōu)度需要建立數(shù)學(xué)模型,使得能夠根據(jù)自變量的數(shù)

值預(yù)測因變量的大小,或者解釋因變量的變化.血壓與年齡剎車距離與車速薪金與資歷、教育程度、工作崗位認(rèn)識程度的限制客觀事物的復(fù)雜人們關(guān)心的(因)變量受另外幾個(自)變量的關(guān)聯(lián)性(非因果性)

影響,并且存在眾多隨機(jī)因素.無法分析對象內(nèi)在的因果關(guān)系回歸分析(RegressionAnalysis)的具體步驟收集一組包含因變量和自變量的數(shù)據(jù)；選定因變量與自變量之間的模型，利用數(shù)據(jù)按照最小二乘準(zhǔn)則計算模型中的系數(shù)；利用統(tǒng)計分析方法對不同的模型進(jìn)行比較，找出與數(shù)據(jù)擬合得最好的模型；檢驗得到的模型是否適合于這組數(shù)據(jù)；利用模型對因變量作出預(yù)測或解釋.

7.1

一元回歸模型

R語言模型formula

I就表示其中所有的運算符都是普通意義上的算術(shù)運算符思考題

例：一元回歸分析示例物理學(xué)家James．D．Forbes試圖通過水的沸點來估計海拔高度，他知道通過氣壓計測得的大氣壓可用于得到海拔高度，氣壓越低，高度越高，他測量了17個地方水的沸點（℉）及大氣壓數(shù)據(jù)，并且對數(shù)據(jù)作了簡單的處理，得到了較為明確的數(shù)學(xué)關(guān)系，所提數(shù)據(jù)如下：X<-matrix(c(194.5,20.79,1.3179,131.79,194.3,20.79,1.3179,131.79,197.9,22.40,1.3502,135.02,198.4,22.67,1.3555,135.55,199.4,23.15,1.3646,136.46,199.9,23.35,1.3683,136.83,200.9,23.89,1.3782,137.82,201.1,23.99,1.3800,138.00,201.4,24.02,1.3806,138.06,201.3,24.01,1.3805,138.05,203.6,25.14,1.4004,140.04,204.6,26.57,1.4244,142.44,209.5,28.49,1.4547,145.47,208.6,27.76,1.4434,144.34,210.7,29.04,1.4630,146.30,211.9,29.88,1.4754,147.54,212.2,30.06,1.4780,147.80),ncol=4,byrow=T,dimnames=list(1:17,c("F","h","log","log100")))forbes=as.data.frame(X)

相關(guān)數(shù)據(jù)如下：首先，畫自變量和因變量之間的散點圖，進(jìn)行數(shù)據(jù)探索分析。從散點圖上發(fā)現(xiàn)X和Y的排列基本是在一條直線附近，那么我們可以假設(shè)X和Y的關(guān)系是線性的。R語言實現(xiàn)有三種方式可以實現(xiàn)最小二乘法的簡單線性回歸，假設(shè)數(shù)據(jù)框為forbeslm(forbes$log100~forbes$F)lm(log100~F,data=forbes)attach(forbes)

;lm(log100~F);#處理完后用detach解除綁定>lm.sol<-lm(log100~F,data=forbes)>summary(lm.sol)通過P值(就是Pr那一列)來查看對應(yīng)的解釋變量x的顯著性，通過將p值與0.05進(jìn)行比較，若該值小于0.05，就可以說該變量與被解釋變量存在顯著的相關(guān)性。MultipleR-squared和AdjustedR-squared這兩個值，又被稱為”擬合優(yōu)度“和”修正的擬合優(yōu)度“，是指回歸方程對樣本的擬合程度，這里我們可以看到，修正的擬合優(yōu)度為0.99946，表示擬合程度超過五成，這個值越高越好。F-statistic，也稱為F檢驗，常用于判斷方程整體的顯著性實驗，其p值<2.2e-16，顯然小于0.05，方程在P=0.05的水平上是通過顯著性檢驗的?；貧w方程

可以將得到的直線方程添加到散點圖上。提取模型信息lm()

的返回值是一個模型擬合結(jié)果對象；技術(shù)上就是

"lm"

的一個結(jié)果列表類。關(guān)于擬合模型的信息可以用能調(diào)用對象類

"lm"

的泛型函數(shù)顯示，提取，圖示等等。這包括add1coefeffectskappapredictresidualsaliasdeviancefamilylabelsprintstepanovadrop1formulaplotprojsummary

其中一些常用的泛型函數(shù)可以簡潔描述如下。anova(object_1,

object_2)比較一個子模型和原模型，并且產(chǎn)生方差分析表。coef(object)提取回歸系數(shù)(矩陣)。全稱：coefficients(object).deviance(object)殘差平方和，如有權(quán)重可加權(quán)。formula(object)提取模型信息。plot(object)產(chǎn)生4個圖，顯式殘差，擬合值和一診斷圖。print(object)簡要打印一個對象的內(nèi)容。常常隱式使用。residuals(object)提取殘差(矩陣)，有權(quán)重時可加權(quán)，省略方式：

resid(object)。step(object)通過增加或者減少模型中的項并且保留層次來選擇合適的模型。在逐步搜索過程中，aiC(akaike信息規(guī)范)值最大的模型將會被返回。summary(object)顯示較詳細(xì)的模型擬合結(jié)果。7.2

多元回歸模型

多元回歸分析示例以R語言MASS包中的cement數(shù)據(jù)為例cement數(shù)據(jù)包括美國Portland市某種水泥產(chǎn)生的熱量數(shù)據(jù)，共13組。水泥凝固時放出的熱量y（響應(yīng)變量）和水泥混合物中的四種化學(xué)成分x1,x2,x3,x4.用全部的變量來解釋熱量y，對應(yīng)的R代碼和結(jié)果:可以發(fā)現(xiàn)，F(xiàn)檢驗的p值很小，說明回歸方程是顯著的，但在各個系數(shù)的t檢驗中，p值均大于0.05，說明每個系數(shù)都是不顯著的，回歸方程可能包括不重要的變量或存在共線性等問題，這樣的模型不能直接應(yīng)用?；貧w系數(shù)的置信區(qū)間

其中object參數(shù)為擬合的回歸模型，parm參數(shù)表示計算置信區(qū)間的變量名稱標(biāo)準(zhǔn)化的回歸系數(shù)標(biāo)準(zhǔn)化回歸系數(shù)，它是在對自變量和因變量同時進(jìn)行標(biāo)準(zhǔn)化處理后所得到的回歸系數(shù)，數(shù)據(jù)經(jīng)過標(biāo)準(zhǔn)化處理后消除了量綱、數(shù)量級等差異的影響，使得不同變量之間具有可比性，因此可以用標(biāo)準(zhǔn)化回歸系數(shù)來比較不同自變量對因變量的作用大小。我們用scale()函數(shù)給出標(biāo)準(zhǔn)化回歸系數(shù)的回歸模型代碼R標(biāo)準(zhǔn)化scale函數(shù)

set.seed(1);x<-runif(10)#Manuallyscaling(x-mean(x))/sd(x)ss=scale(x)identical(attr(ss,"scaled:center"),mean(x))identical(attr(ss,"scaled:scale"),sd(x))[1]-0.8717643-0.52873940.11708951.1960620-1.07712101.1644732[1]TRUE[1]TRUE未標(biāo)準(zhǔn)化回歸系數(shù)體現(xiàn)的是自變量變化對因變量的絕對作用大小標(biāo)準(zhǔn)化回歸系數(shù)反映的是不同自變量對因變量的相對作用大小，可以顯示出不同自變量對因變量影響的重要性。如果用標(biāo)準(zhǔn)化回歸系數(shù)構(gòu)建方程，得到的結(jié)論是有偏的，因為此時自變量和因變量的數(shù)據(jù)都發(fā)生了轉(zhuǎn)化，成為了標(biāo)準(zhǔn)化數(shù)據(jù)，因此標(biāo)準(zhǔn)化回歸系數(shù)不能直接用于回歸方程預(yù)測。7.3

回歸模型的擬合優(yōu)度

7.4

回歸模型診斷

擬合好一個回歸模型并不是直接就能夠使用了，而是先要對回歸分析中的假設(shè)（高斯-馬爾可夫條件）以及數(shù)據(jù)進(jìn)行檢驗與分析。檢驗回歸分析中的假設(shè)是否合理。對數(shù)據(jù)的診斷7.4.1殘差的分類

R語言中可利用rstandard()調(diào)用回歸模型的標(biāo)準(zhǔn)化殘差。

7.4.2線性假設(shè)診斷判斷線性假設(shè)，可以對lm()函數(shù)返回的模型對象使用plot()函數(shù)。該函數(shù)可生成評價模型擬合情況的多個圖形，默認(rèn)返回四個圖形?？墒褂胮lot()函數(shù)中的殘差-擬合值（Residualvsfitted）圖判斷回歸模型中線性關(guān)系是否成立。plot(lm.sol,which=1)殘差與預(yù)測值的變化之間不存在明顯的規(guī)律，可以認(rèn)為基本滿足線性關(guān)系。ResidualsvsFitted7.4.3

殘差分析和異常點檢測高斯-馬爾可夫條件中涉及隨機(jī)誤差項的假設(shè)較多，而誤差是無法被觀測的，因此無法直接利用誤差檢驗?zāi)Ｐ图僭O(shè)是否成立。我們可以通過判斷殘差是否符合誤差項的某些性質(zhì)，達(dá)到檢驗?zāi)Ｐ褪欠穹舷鄳?yīng)假設(shè)的目的。在R中，我們可以使用plot()函數(shù)中的正態(tài)Q-Q圖（NormalQQ-plot）判斷Forbes數(shù)據(jù)判斷殘差是否服從正態(tài)分布。plot(lm.sol,which=2)除開12號樣本點外，其他數(shù)據(jù)點分布趨于一條直線,說明殘差近似服從正態(tài)分布。鏈接：不同情形下的NormalQQPlots在R中，我們還可以使用plot()函數(shù)中的尺度-位置（Scale-Location）圖顯示的殘差分布情況，判斷同方差性：若滿足不變方差假設(shè)，那么在位置尺度圖中，水平線周圍的點應(yīng)該隨機(jī)分布。plot(lm.sol,which=3)大部分?jǐn)?shù)據(jù)點呈現(xiàn)出隨機(jī)的分布，但部分?jǐn)?shù)據(jù)（如12號點）的標(biāo)準(zhǔn)化殘差平方根過大，說明原數(shù)據(jù)集中不滿足同方差的假設(shè)。

圖中12號樣本點的標(biāo)準(zhǔn)化殘差較大，屬于離群點，但該點杠桿值和Cook距離不大，故不是高杠桿點。我們可以用plot()函數(shù)中的Cook距離圖進(jìn)一步觀察各點的Cook距離plot(lm.sol,which=4)一般來說，小樣本的Cook距離大于1，大樣本的Cook距離大于4/n算作高杠桿點。對于離群點的處理，可以去掉12號樣本點再擬合可以看出，去掉12號樣本點后回歸方程系數(shù)沒有太大的變化，但是系數(shù)的標(biāo)準(zhǔn)差和殘差標(biāo)準(zhǔn)差有很大的變化，平均減少約3倍左右，R平方的值也有提高。7.4.4多重共線性檢測多重共線性是指在多元回歸模型中，自變量（解釋變量）之間存在一定的線性相關(guān)關(guān)系而使模型估計失真或估計精度不高。在Gauss-Markov條件滿足的條件下，多重共線性不改變估計量的無偏性，但各共線變量參數(shù)的最小二乘法估計值的精度很低，容易使結(jié)果變得不顯著，進(jìn)而無法正確判斷各自變量對因變量的影響。R中常見驗證變量之間是否存在共線性的方法：條件數(shù)法、方差膨脹因子(VIF)、Farrar-Glauber方法

多重共線性問題#ExactCollinearitygen_exact_collin_data=function(num_samples=100){x1=rnorm(n=num_samples,mean=80,sd=10)x2=rnorm(n=num_samples,mean=70,sd=5)x3=2*x1+4*x2+3y=3+x1+x2+rnorm(n=num_samples,mean=0,sd=1)data.frame(y,x1,x2,x3)}set.seed(1234)exact_collin_data=gen_exact_collin_data()當(dāng)我們嘗試使用所有預(yù)測變量在R中擬合回歸模型時會發(fā)生什么？>head(exact_collin_data)yx1x2x31143.487267.9293472.07262427.14922154.097582.7742967.62641439.05423164.359990.8444170.32997466.00874127.731456.5430267.48761386.03655153.472984.2912565.87001435.06256159.656085.0605670.83495456.4609>exact_collin_fit=lm(y~x1+x2+x3,data=exact_collin_data)>summary(exact_collin_fit)Call:lm(formula=y~x1+x2+x3,data=exact_collin_data)Residuals:Min1QMedian3QMax-3.1933-0.67340.10600.58232.7650Coefficients:(1notdefinedbecauseofsingularities)EstimateStd.ErrortvaluePr(>|t|)(Intercept)1.2829981.5371400.8350.406x11.0080510.009633104.649<2e-16***x21.0176640.01874754.284<2e-16***x3NANANANA---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.9624on97degreesoffreedomMultipleR-squared:0.9929, AdjustedR-squared:0.9928F-statistic:6809on2and97DF,p-value:<2.2e-16

>X=cbind(1,as.matrix(exact_collin_data[,-1]))#inverse>solve(t(X)%*%X)Errorinsolve.default(t(X)%*%X):systemiscomputationallysingular:reciprocalconditionnumber=1.44535e-18fit1=lm(y~x1+x2,data=exact_collin_data)fit2=lm(y~x1+x3,data=exact_collin_data)fit3=lm(y~x2+x3,data=exact_collin_data)all.equal(fitted(fit1),fitted(fit2))##[1]TRUE我們看到三個模型的擬合值完全相同。WHY？這是因為x3包含來自x1和x2的所有信息的結(jié)果。只要模型中包含x1或x2之一，x3就可以用于從未包含的變量中恢復(fù)信息。雖然擬合值都相同，但模型的估計系數(shù)卻大不相同。x2的符號在fit2、fit3兩個模型中甚至是相反的！>coef(fit1)(Intercept)x1x21.2829981.0080511.017664>coef(fit2)(Intercept)x1x30.51975030.49921950.2544160>coef(fit3)(Intercept)x2x3-0.2290789-0.99843900.5040257

在(0,10)之間,沒有多重共線性.[10,100)有較強(qiáng)多重共線性.大于100，有嚴(yán)重的多重共線性.例：考慮一個有6個回歸自變量的線性回歸問題，原始數(shù)據(jù)列在表中，這里共有12組數(shù)據(jù)，除第一組外，自變量x1,x2,…,x6的其余11組數(shù)據(jù)滿足線性關(guān)系：x1+x2+x3+x4=10，試用求矩陣條件數(shù)的方法，分析出自變量間存在多重共線性。例：R語言多重共線性collinear<-data.frame(y=c(10.006,9.737,15.087,8.422,8.625,16.289,5.958,9.313,12.96,5.541,8.756,10.937),x1=rep(c(8,0,2,0),c(3,3,3,3)),x2=rep(c(1,0,7,0),c(3,3,3,3)),x3=rep(c(1,9,0),c(3,3,6)),x4=rep(c(1,0,1,10),c(1,2,6,3)),x5=c(0.541,0.13,2.116,-2.397,-0.046,0.365,1.996,0.228,1.38,-0.798,0.257,0.44),x6=c(-0.099,0.07,0.115,0.252,0.017,1.504,-0.865,-0.055,0.502,-0.399,0.101,0.432))先計算出對應(yīng)矩陣的特征值：可以利用定義計算條件數(shù)：e$values[1]/e$values[6][1]2195.908#條件數(shù)大于100有嚴(yán)重的多重共線性。

e$vectors[,6][1]0.4476797190.4211402800.5416891240.5733718720.0060521270.002166594也可直接用kappa()函數(shù)計算條件數(shù):kappa(xx,exact=TRUE)##[1]2196方差膨脹因子(VIF)

Farrar-Glauber方法

用mctest包的omcdiag函數(shù)進(jìn)行檢驗發(fā)現(xiàn)標(biāo)準(zhǔn)化行列式的值為0.001，這是非常小的?？ǚ綑z驗統(tǒng)計量的計算值為56.951，它非常顯著，因此暗示模型中存在多重共線性。這促使我們對多重共線性的位置進(jìn)行下一步Farrar-Glauber檢驗（F檢驗）mctest::imcdiag(lmcoliner)VIF、TOL和Wi列分別為方差膨脹因子、容差和F檢驗的結(jié)果。F檢驗表明:變量“x1”或“x2”或“x3”或“x4”可能是多重共線性的根本原因。正如預(yù)期的那樣，“x1”和“x2,x3,x4”之間高度偏相關(guān)，具有統(tǒng)計顯著性。Farrar-Glauber檢驗指出這幾個變量是所有多重共線性問題的根本原因。我們?nèi)サ簟皒2,x3,x4”，只保留”x1”,對應(yīng)的VIF指標(biāo)如下：VIF的值較小，說明剩余變量中不存在共線性問題

模型診斷-car包實用函數(shù)qqPlot()#分位數(shù)比較圖durbinWastonTest()#誤差的獨立性Durbin-Waston檢驗crPlots()#成分和殘差圖,偏殘差圖，用于檢測因變量與每一個自變量之間的線性關(guān)系是否顯著outlierTest()#Bonferroni離群點檢驗avPlots()#添加的變量圖形，對于每個預(yù)測變量Xk，繪制Xk在其他k-1個預(yù)測變量上回歸的殘差值influencePlot()#回歸影響圖,將離群點、杠桿值和強(qiáng)影響點的信息整合到一幅圖形中ncvTest()#同方差性檢測vif()#檢驗多重共線性7.5模型選擇當(dāng)預(yù)測模型中有許多預(yù)測變量時，模型選擇方法允許自動選擇預(yù)測變量的最佳組合以構(gòu)建最佳的回歸模型。一種很自然的模型選擇策略是窮舉所有可能的預(yù)測變量組合，然后從中選擇最佳的回歸模型。這種方法稱為最佳子集回歸或全子集回歸。另外一個是逐步回歸，即按順序比較包含不同個數(shù)預(yù)測變量組合的線性回歸模型。7.5.1最佳子集回歸最佳子集回歸可用leaps包中的regsubsets()函數(shù)實現(xiàn)參數(shù)nvmax的值表示模型中包含的最大預(yù)測變量的個數(shù)。最好的2變量模型只包含x1和x2（y~x1+x2）最好的3變量模型則是（y~x1+x2+x4）這些最佳模型哪個最好調(diào)整后的R2告訴我們最好的模型是第3個模型即具有3個預(yù)測變量的模型(y~x1+x2+x4)。但是，使用BIC和Cp標(biāo)準(zhǔn)，則應(yīng)該選擇具有2個變量的模型（y~x1+x2）?？筛鶕?jù)實際情況，選擇“最佳”的模型。Allmodelsarewrong,butsomeareuseful7.5.2逐步回歸(stepwiseregression)逐步回歸中，模型會一次添加或者刪除一個變量，直到達(dá)到某個判停準(zhǔn)則為止。向前逐步回歸-forwardstepwise-每次添加一個預(yù)測變量到模型中，直到添加變量不會使模型有所改進(jìn)為止。向后逐步回歸-backwardstepwise-從模型包含所有預(yù)測變量開始，一次刪除一個變量直到降低模型質(zhì)量為止。向前向后逐步回歸（stepwisestepwise，通常稱作逐步回歸），結(jié)合了向前逐步回歸和向后逐步回歸的方法，變量每次進(jìn)入一個，但是每一步中，變量都會被重新評價，對模型沒有貢獻(xiàn)的變量將會被刪除，預(yù)測變量可能會被添加、刪除好幾次，直到獲得最優(yōu)模型為止。逐步回歸

全子集回歸可用leaps包的regsubsets()函數(shù)實現(xiàn)。R自帶的step函數(shù)實現(xiàn)了最基本的逐步回歸模型。MASS包中的stepAIC()函數(shù)可以實現(xiàn)逐步回歸模型（向前、向后和向前向后），依據(jù)的是精確AIC準(zhǔn)則。Call:lm(formula=y~x1+x2+x3+x4,data=cement)Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)62.405470.07100.8910.3991

x11.55110.74482.0830.0708.x20.51020.72380.7050.5009x30.10190.75470.1350.8959x4-0.14410.7091-0.2030.8441---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘

’1Residualstandarderror:2.446on8degreesoffreedomMultipleR-squared:0.9824,AdjustedR-squared:0.9736F-statistic:111.5on4and8DF,p-value:4.756e-07

#解釋變量不顯著考慮所有可用預(yù)測變量的模型作為初始模型Step處理和其他軟件不同，R軟件根據(jù)AIC準(zhǔn)則，自動選擇合適的解釋變量.>summary(lm.step)Call:lm(formula=y~x1+x2+x4,data=cement)Residuals:Min1QMedian3QMax-3.0919-1.80160.25621.28183.8982Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)71.648314.14245.0660.000675***x11.45190.117012.4105.78e-07***x20.41610.18562.2420.051687.x4-0.23650.1733-1.3650.205395---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:2.309on9degreesoffreedomMultipleR-squared:0.9823, AdjustedR-squared:0.9764F-statistic:166.8on3and9DF,p-value:3.323e-08drop1(lm.step)與改進(jìn)：>drop1(lm.step)SingletermdeletionsModel:y~x1+x2+x4DfSumofSqRSSAIC<none>47.9724.974x11820.91868.8860.629x2126.7974.7628.742x419.9357.9025.420刪除變量會改變根據(jù)AIC最小準(zhǔn)則得到的step結(jié)果，根據(jù)drop1分析，刪除x4使AIC值上升最小，為盡量遵守AIC準(zhǔn)則的同時，提高解釋變量的顯著性，刪除x4是“無奈之舉”。drop1(lm.step)與改進(jìn)：newsol=lm(y~x1+x2,data=cement)summary(newsol)Call:lm(formula=y~x1+x2,data=cement)Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)52.577352.2861723.005.46e-10***x11.468310.1213012.112.69e-07***x20.662250.0458514.445.03e-08***Residualstandarderror:2.406on10degreesoffreedomMultipleR-squared:0.9787,AdjustedR-squared:0.9744F-statistic:229.5on2and10DF,p-value:4.407e-09

比較：step()處理后的殘差和可決系數(shù)：Residualstandarderror:2.309AdjustedR-squared:0.9764在MASS包中，stepAIC()函數(shù)實現(xiàn)逐步回歸模型stepAIC(object,scope,scale=0,direction=c("both",

"backward","forward"),trace=1,keep=NULL,steps=1000,use.start=FALSE,k=2,...)“backward”表示從給定模型中依次刪除預(yù)測變量，它產(chǎn)生一系列復(fù)雜性遞減的模型，直到獲得最佳模型；“forward”表示使用數(shù)據(jù)集中的所有可用變量，將預(yù)測變量按順序添加到給定模型中，它產(chǎn)生一系列復(fù)雜性增加的模型，直到獲得最佳模型；“both”（默認(rèn)）表示向前-向后搜索，在每一步，決定是包含還是排除預(yù)測變量，與之前的模式不同，先前被排除/包括的預(yù)測變量可以稍后被包括/排除。stepAIC()函數(shù)默認(rèn)考慮AIC作為度量準(zhǔn)則

第一步返回的結(jié)果是”Step：AIC=29.77”，第二步移除了x3，因為它給出了最低的BIC,此時BIC變?yōu)椤癝tep:AIC=27.23”。最后一步移除了x4,BIC值進(jìn)一步減少，改進(jìn)為: