商務(wù)數(shù)據(jù)分析與統(tǒng)計建模：chap2.2 多元回歸分析R實現(xiàn)

1、 2. 1 多元線性模型多元線性模型 2.2 變量選擇變量選擇 2.3 回歸診斷回歸診斷 2.4 回歸預(yù)測回歸預(yù)測 第第2章多元線性模型章多元線性模型 1 2.1 多元線性模型多元線性模型 1.模型定義模型定義 其中其中x1,xk是非隨機的自變量是非隨機的自變量, y是隨機的因是隨機的因 變量變量, 0是常數(shù)項是常數(shù)項, 1, k是回歸系數(shù)是回歸系數(shù), 是是 隨機誤差項隨機誤差項. 01 1 ,(2.1) kk yxx 2 1.模型定義模型定義 模型模型(2.1)的樣本形式的樣本形式 其中其中i=1,n,表示有表示有n組觀測值組觀測值. 01 1 ,(2.2) iikkii yxx 3 1.模

2、型定義模型定義 模型模型(2.2)的矩陣形式的矩陣形式 其中其中 ,(2.3)yX 4 111101 212212 1 1 1 , 1 k k nnknkn yxx yxx yX yxx 例例2.1(數(shù)據(jù)文件為數(shù)據(jù)文件為eg2.1) 5 主編：費宇 表表2.1 抽樣調(diào)查得到的抽樣調(diào)查得到的36個人的數(shù)據(jù)資料個人的數(shù)據(jù)資料 2021-5-12 yx1x2x3x4x5yx1x2x3x4x5 79220 14010 9811515女女71120 11460 83758女女 79670 13260 98268男男91520 22260 81316男男 186320 81240 9619919男男762

3、20 12510 81012女女 161945 46260 9612019男男74420 12510 811312女女 68420 11460 857212女女121320 23010 651916男男 75320 14010 855915女女77570 17010 646912男男 模型（模型（2.3）的參數(shù)）的參數(shù) 的最小二乘估計的最小二乘估計 的的最小二乘估計最小二乘估計 2.模型的參數(shù)估計和檢驗?zāi)Ｐ偷膮?shù)估計和檢驗 6 1 (),(2.5) TT X XX y 2 .(2.7) 1 T e e nk 2 回歸方程的顯著性檢驗回歸方程的顯著性檢驗 其中其中 是回歸平方和，而是回歸平方和，

4、而 是殘差平方和，拒絕域為是殘差平方和，拒絕域為 2.模型的參數(shù)估計和檢驗?zāi)Ｐ偷膮?shù)估計和檢驗 7 01 :0 / ( ,1),(2.8) /(1) k H SSR k FF k nk SSEnk 2 1 () n ii i SSEyy 2 1 () n i i SSRyy ( ,1)FF k nk 2.模型的參數(shù)估計和檢驗?zāi)Ｐ偷膮?shù)估計和檢驗 8 ( ) (1,1),(2.9) /(1) j j SSESSE FFnk SSEnk ( ) j SSE (1,1) j FFnk 例例2.1續(xù)續(xù)(數(shù)據(jù)文件為數(shù)據(jù)文件為eg2.1) 9 主編：費宇 建立建立y關(guān)于關(guān)于x1、x2、x3和和x4的線性回

5、歸方程的線性回歸方程,并對方程和回歸并對方程和回歸 系數(shù)進(jìn)行顯著性檢驗系數(shù)進(jìn)行顯著性檢驗. 2021-5-12 yx1x2x3x4x5yx1x2x3x4x5 79220 14010 9811515女女71120 11460 83758女女 79670 13260 98268男男91520 22260 81316男男 186320 81240 9619919男男76220 12510 81012女女 161945 46260 9612019男男74420 12510 811312女女 68420 11460 857212女女121320 23010 651916男男 75320 14010 85

6、5915女女77570 17010 646912男男 #eg2.1回歸分析回歸分析 #打開數(shù)據(jù)文件打開數(shù)據(jù)文件eg2.1.xls,選取選取A1:F37區(qū)域區(qū)域,然后復(fù)制然后復(fù)制 data2.1-read.table(clipboard,header=T) #將將eg2.1.xls數(shù)據(jù)讀入到數(shù)據(jù)讀入到data2.1中中 lm.salary|t|) Estimate Std. Error t value Pr(|t|) (Intercept) 48386.0620 11237.2882 4.306 0.000155 (Intercept) 48386.0620 11237.2882 4.306 0

7、.000155 * * * * x1 1.6831 0.1302 12.929 5.01e-14 x1 1.6831 0.1302 12.929 5.01e-14 * * * * x2 -34.5520 130.2602 -0.265 0.792570 x2 -34.5520 130.2602 -0.265 0.792570 x3 -13.0004 13.7882 -0.943 0.353043 x3 -13.0004 13.7882 -0.943 0.353043 x4 808.3223 547.8017 1.476 0.150144 x4 808.3223 547.8017 1.476 0

8、.150144 例例2.1續(xù)續(xù)(數(shù)據(jù)文件為數(shù)據(jù)文件為eg2.1) 2021-5-12 11 回歸分析結(jié)果為回歸分析結(jié)果為 - Signif. codes: 0 Signif. codes: 0 * * * * 0.001 0.001 * * * 0.01 0.01 * * 0.05 . 0.05 . 0.1 1 0.1 1 Residual Residual standard error: 7858 on 31 degrees of freedomstandard error: 7858 on 31 degrees of freedom Multiple R-squared: 0.919, A

9、djusted R-squared: 0.9086 Multiple R-squared: 0.919, Adjusted R-squared: 0.9086 F-statistic: 87.95 on 4 and 31 DF, p-value: 2.2e-16F-statistic: 87.95 on 4 and 31 DF, p-value: 2.2e-16 例例2.1續(xù)續(xù)(數(shù)據(jù)文件為數(shù)據(jù)文件為eg2.1) 2021-5-12 12 回歸方程的回歸方程的F值為值為87.95,相應(yīng)的相應(yīng)的p值為值為2.2 10-16,說明說明 回歸方程是顯著的回歸方程是顯著的;但但t檢驗對應(yīng)的檢驗對應(yīng)的p值

10、則顯示值則顯示:常數(shù)常數(shù) 項和項和x1是顯著的是顯著的,而而x2、x3和和x4不顯著不顯著. 最優(yōu)模型一般滿足最優(yōu)模型一般滿足2個條件個條件 （1）模型反映了變量間的真實關(guān)系）模型反映了變量間的真實關(guān)系 （2）模型包含的變量盡量）模型包含的變量盡量少少 例例2.1（續(xù)）（續(xù)） 建立建立y關(guān)于關(guān)于x1、x2、x3和和x4的線性回歸方程的線性回歸方程,并對方程和回并對方程和回 歸系數(shù)進(jìn)行歸系數(shù)進(jìn)行顯著性檢驗顯著性檢驗（逐步回歸建立（逐步回歸建立“最優(yōu)方程最優(yōu)方程”）. 2.2 變量選擇變量選擇 2021-5-12 13 #假設(shè)假設(shè)eg2.1.xls中的數(shù)據(jù)已經(jīng)讀入到中的數(shù)據(jù)已經(jīng)讀入到data2.1

11、中中, lm.salary-lm(yx1+x2+x3+x4,data=data2.1) #建立全變量建立全變量 回歸方程回歸方程 lm.step-step(lm.salary,direction=both) #用用“一切子集一切子集 回歸法回歸法”來進(jìn)行逐步回歸來進(jìn)行逐步回歸 direction是確定逐步搜索的是確定逐步搜索的 方向方向:both是是“一切子集回歸法一切子集回歸法”, forward是是“向前向前 法法”,backward是是“向后法向后法”,默認(rèn)值是默認(rèn)值是both.所以這個所以這個 回歸過程可以簡寫為回歸過程可以簡寫為lm.step-step(lm.salary) 2.2

12、變量選擇變量選擇 2021-5-12 14 例例2.1(續(xù)續(xù))逐步回歸程序逐步回歸程序 Start: AIC=650.41Start: AIC=650.41 y y x1 + x2 + x3 + x4 x1 + x2 + x3 + x4 Df Sum of Sq RSS AICDf Sum of Sq RSS AIC - x2 1 4.3448e+06 1.9186e+09 648.49- x2 1 4.3448e+06 1.9186e+09 648.49 - x3 1 5.4896e+07 1.9692e+09 649.43- x3 1 5.4896e+07 1.9692e+09 649.4

13、3 1.9143e+09 650.41 1.9143e+09 650.41 - x4 1 1.3445e+08 2.0487e+09 650.85- x4 1 1.3445e+08 2.0487e+09 650.85 - x1 1 1.0323e+10 1.2237e+10 715.19- x1 1 1.0323e+10 1.2237e+10 715.19 Step: AIC=648.49Step: AIC=648.49 y y x1 + x3 + x4 x1 + x3 + x4 Df Sum of Sq RSS AIC Df Sum of Sq RSS AIC - x3 1 6.2078e+

14、07 1.9807e+09 647.64- x3 1 6.2078e+07 1.9807e+09 647.64 2.2 變量選擇變量選擇 2021-5-12 15 例例2.1(續(xù)續(xù))回歸結(jié)果回歸結(jié)果 1.9186e+09 648.49 - x4 1 1.3011e+08 2.0487e+09 648.85 + x2 1 4.3448e+06 1.9143e+09 650.41 - x1 1 1.0341e+10 1.2259e+10 713.26 Step: AIC=647.64 y x1 + x4 Df Sum of Sq RSS AIC 1.9807e+09 647.64 + x3 1 6

15、.2078e+07 1.9186e+09 648.49 + x2 1 1.1527e+07 1.9692e+09 649.43 - x4 1 2.9640e+08 2.2771e+09 650.66 - x1 1 1.1654e+10 1.3635e+10 715.09 2.2 變量選擇變量選擇 2021-5-12 16 Call: lm(formula = y x1 + x4, data = data2.1) Residuals: Min 1Q Median 3Q Max -13632 -4759 -615 1761 25076 Coefficients: Estimate Std. Err

16、or t value Pr(|t|) (Intercept) 42097.165 5265.218 7.995 3.18e-09 * x1 1.631 0.117 13.934 2.22e-15 * x4 1039.260 467.671 2.222 0.0332 * - Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 7747 on 33 degrees of freedom Multiple R-squared: 0.9162, Adjusted R-squared: 0.9111 F-stat

17、istic: 180.4 on 2 and 33 DF, p-value: 2.2e-16 回歸模型匯總信息回歸模型匯總信息:summary(lm.step) 2021-5-12 17 14 42097.165+1.6311039.260yxx 1.殘差分析和異常點探測殘差分析和異常點探測 殘差向量殘差向量e是模型中隨機誤差項是模型中隨機誤差項 的估計，的估計， 殘差分析可以診斷模型的基本假定是否成殘差分析可以診斷模型的基本假定是否成 立。立。 2.3 2.3 模型診斷模型診斷 2021-5-12 18 eyyyX 例例2.2 計算計算例例2.1得到的逐步回歸模型得到的逐步回歸模型lm.ste

18、p的普的普 通殘差和標(biāo)準(zhǔn)化殘差通殘差和標(biāo)準(zhǔn)化殘差,判斷可能存在的異常點判斷可能存在的異常點, 畫出相應(yīng)的殘差散點圖畫出相應(yīng)的殘差散點圖,并直觀判斷模型的基并直觀判斷模型的基 本假定是否成立本假定是否成立. 2.3 模型診斷模型診斷 2021-5-12 19 #假設(shè)由例假設(shè)由例2.1已經(jīng)得到逐步回歸模型已經(jīng)得到逐步回歸模型lm.step y.res-residuals (lm.salary) #計算回歸模型計算回歸模型lm.salary的普通殘差的普通殘差 y.rst-rstandard(lm.step) #計算回歸模型計算回歸模型lm.step的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差 print(y.rst)

19、 #輸出回歸模型輸出回歸模型lm.step的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差y.rst y.fit-predict(lm.step) #計算回歸模型計算回歸模型lm.step的預(yù)測值的預(yù)測值 plot(y.res y.fit) #繪制以普通殘差為縱坐標(biāo)繪制以普通殘差為縱坐標(biāo),預(yù)測值為橫坐標(biāo)的殘差散點圖預(yù)測值為橫坐標(biāo)的殘差散點圖 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測值為橫坐標(biāo)的殘差散點圖預(yù)測值為橫坐標(biāo)的殘差散點圖 2.3 模型診斷模型診斷 2021-5-12 20 采用采用residuals( ),rstandard( )和和rstudent( )來

20、分別計算來分別計算 普通殘差普通殘差,標(biāo)準(zhǔn)化殘差和學(xué)生化殘差標(biāo)準(zhǔn)化殘差和學(xué)生化殘差. 1 2 3 4 5 6 1 2 3 4 5 6 -0.17481171 1.03650457 -1.64144064 -0.17481171 1.03650457 -1.64144064 3.435090883.43509088 -0.69388898 0.21730074 -0.69388898 0.21730074 7 8 9 10 11 12 7 8 9 10 11 12 -0.28221956 -0.57391074 -1.13860082 0.14942833 -0.08074416 -0.7031

21、3518 -0.28221956 -0.57391074 -1.13860082 0.14942833 -0.08074416 -0.70313518 13 14 15 16 17 18 13 14 15 16 17 18 -0.08074416 -1.90879266 0.86519102 -0.04532806 -0.63719129 -0.69529295 -0.08074416 -1.90879266 0.86519102 -0.04532806 -0.63719129 -0.69529295 19 20 21 22 23 24 19 20 21 22 23 24 0.27405524

22、 -0.46453524 0.16499595 -0.07189545 0.22203402 0.41883966 0.27405524 -0.46453524 0.16499595 -0.07189545 0.22203402 0.41883966 25 26 27 28 29 30 25 26 27 28 29 30 0.25844627 -0.12669958 -0.01045782 -0.29460957 -0.15151718 -0.89674431 0.25844627 -0.12669958 -0.01045782 -0.29460957 -0.15151718 -0.89674

23、431 31 32 33 34 35 36 31 32 33 34 35 36 0.47424184 0.64162319 -0.84880014 0.18495911 0.47424184 0.64162319 -0.84880014 0.18495911 3.324239703.32423970 -0.62075537 -0.62075537 2.3 模型診斷模型診斷 2021-5-12 21 回歸模型回歸模型lm.step的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差y.rst如下如下 第第4號和號和35號是異常點號是異常點 2021-5-12 22 圖圖2.1 2.1 例例2.22.2中的普通殘差圖中的普通

24、殘差圖( (左左) )和標(biāo)準(zhǔn)化殘差圖和標(biāo)準(zhǔn)化殘差圖( (右右) ) 例例2.3 通過通過方差穩(wěn)定變換來更新例方差穩(wěn)定變換來更新例2.1得到的逐步得到的逐步 回歸模型回歸模型lm.step,并計算更新后模型的標(biāo)準(zhǔn)化并計算更新后模型的標(biāo)準(zhǔn)化 殘差殘差,畫出相應(yīng)的殘差散點圖畫出相應(yīng)的殘差散點圖,并直觀判斷模型并直觀判斷模型 的基本假定是否成立的基本假定是否成立. 2.3 模型診斷模型診斷 2021-5-12 23 #假設(shè)由例假設(shè)由例2.1已經(jīng)得到逐步回歸模型已經(jīng)得到逐步回歸模型lm.step lm.step_new-update(lm.step,log(.).) #對模型進(jìn)行對數(shù)變換對模型進(jìn)行對數(shù)變

25、換 y.rst-rstandard(lm.step_new) #計算新回歸模型計算新回歸模型lm.step_new的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差 y.fit-predict(lm.step_new) #計算回歸模型計算回歸模型lm.step_new的預(yù)測值的預(yù)測值 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測值為橫坐標(biāo)的殘差散點圖預(yù)測值為橫坐標(biāo)的殘差散點圖 (見圖見圖2-2) 采用對數(shù)變換來解決方差非齊問題采用對數(shù)變換來解決方差非齊問題 2021-5-12 24 對數(shù)變換后：只有第對數(shù)變換后：只有第35號是異常點號是異常點 2021-5-12 25 圖

26、圖2.2 例例2.3中的標(biāo)準(zhǔn)化殘差圖中的標(biāo)準(zhǔn)化殘差圖 lm.salary-lm(log(y)x1+x2+x3+x4,data= data2.1-c(4,35),) #去掉第去掉第4號和號和 35號觀測值再建立全變量回歸方程號觀測值再建立全變量回歸方程 lm.step-step(lm.salary,direction=both) #用用“一切子集回歸法一切子集回歸法”來進(jìn)行來進(jìn)行 逐步回歸逐步回歸 y.rst-rstandard(lm.step) #計算回歸模型計算回歸模型lm.step的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差 y.fit-predict(lm.step) #計算回歸模型計算回歸模型lm.ste

27、p的預(yù)測值的預(yù)測值 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測值為橫坐標(biāo)的殘差散點預(yù)測值為橫坐標(biāo)的殘差散點 圖圖 去掉去掉4 4和和3535號觀測值號觀測值再回歸再回歸 2021-5-12 26 殘差全部落在殘差全部落在-2，2區(qū)域內(nèi)區(qū)域內(nèi) 2021-5-12 27 圖圖2.3 例例2.3中的標(biāo)準(zhǔn)化殘差圖中的標(biāo)準(zhǔn)化殘差圖:去掉去掉4和和35號觀測值號觀測值 2.回歸診斷：一般的方法回歸診斷：一般的方法 殘差分析無法分析模型的影響點殘差分析無法分析模型的影響點,即探測哪些點對模型的即探測哪些點對模型的 推斷有重要影響推斷有重要影響,本節(jié)給出本節(jié)給

28、出的的回歸診斷方法回歸診斷方法,可以診斷模型可以診斷模型 的基本假定是否成立的基本假定是否成立,哪些值是異常點哪些值是異常點,哪些點是強影響點哪些點是強影響點. 在在R中中,函數(shù)函數(shù)plot( )和和influence.measures( )可以用來繪制診斷可以用來繪制診斷 圖和計算診斷統(tǒng)計量圖和計算診斷統(tǒng)計量 2.3 模型診斷模型診斷 2021-5-12 28 例例2.4 對例對例2.3得到的逐步回歸模型得到的逐步回歸模型lm.step_new進(jìn)進(jìn) 行回歸診斷分析行回歸診斷分析. 2.3 模型診斷模型診斷 2021-5-12 29 #假定由例假定由例2.3已經(jīng)獲得模型已經(jīng)獲得模型lm.ste

29、p_new par(mfrow=c(2,2) #在一個在一個22網(wǎng)格中創(chuàng)建網(wǎng)格中創(chuàng)建4個繪圖區(qū)個繪圖區(qū) plot(lm.step_new) #繪制模型診斷圖繪制模型診斷圖 influence.measures(lm.step_new) #計算各個觀測值的診計算各個觀測值的診 斷統(tǒng)計量斷統(tǒng)計量 運行上述程運行上述程進(jìn)行進(jìn)行序可得回歸診斷圖序可得回歸診斷圖(圖圖2.4)和如下和如下 36個觀測值對應(yīng)的診斷統(tǒng)計量的值個觀測值對應(yīng)的診斷統(tǒng)計量的值. 2021-5-12 30 Influence measures of lm(formula = log(y) x1 + x4, data = data2.

30、1) : dfb.1_ dfb.x1 dfb.x4 dffit cov.r cook.d hat inf 1 0.027388 0.05593 -0.057229 -0.08411 1.161 2.42e-03 0.0646 2 0.333715 0.07828 -0.277049 0.35202 1.088 4.11e-02 0.0953 3 0.348810 -3.49426 0.912364 -3.86083 1.255 3.93e+00 0.6054* 4 -0.580142 0.53947 0.390696 1.14942 0.649 3.63e-01 0.1404* 5 -0.07

31、0573 0.03314 0.018032 -0.15068 1.064 7.64e-03 0.0324 . 34 -0.085786 0.13422 0.042539 0.22772 1.218 1.76e-02 0.1287 35 -0.437405 -0.28259 0.631117 0.93919 0.341 2.02e-01 0.0520* 36 -0.053140 0.01007 0.020734 -0.10486 1.096 3.74e-03 0.0312 第第3,4和和35號觀測值可能是異常點和強影響點號觀測值可能是異常點和強影響點. 2021-5-12 31 圖圖2.4 例例

32、2.4的回歸診斷圖的回歸診斷圖 Influence measures ofInfluence measures of lm(formula = log(y) x1 + x4, data = data2.1) : lm(formula = log(y) x1 + x4, data = data2.1) : dfb.1_ dfb.x1 dfb.x4 dffit cov.r cook.d hat dfb.1_ dfb.x1 dfb.x4 dffit cov.r cook.d hat infinf 1 0.027388 0.05593 -0.057229 -0.08411 1.161 2.42e-03

33、 0.06461 0.027388 0.05593 -0.057229 -0.08411 1.161 2.42e-03 0.0646 2 0.333715 0.07828 -0.277049 0.35202 1.088 4.11e-02 0.09532 0.333715 0.07828 -0.277049 0.35202 1.088 4.11e-02 0.0953 3 0.348810 -3.49426 0.912364 -3.86083 1.255 3.93e+00 0.60543 0.348810 -3.49426 0.912364 -3.86083 1.255 3.93e+00 0.6054* * 4 -0.580142 0.53947 0.390696 1.14942 0.649 3.63e-01 0.14044 -0.580142 0.53947 0.390696 1.14942 0.649 3.63e-01