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

1、 2. 1 多元線性模型多元線性模型 2.2 變量選擇變量選擇 2.3 回歸診斷回歸診斷 2.4 回歸預(yù)測(cè)回歸預(yù)測(cè) 第第2章多元線性模型章多元線性模型 1 2.1 多元線性模型多元線性模型 1.模型定義模型定義 其中其中x1,xk是非隨機(jī)的自變量是非隨機(jī)的自變量, y是隨機(jī)的因是隨機(jī)的因 變量變量, 0是常數(shù)項(xiàng)是常數(shù)項(xiàng), 1, k是回歸系數(shù)是回歸系數(shù), 是是 隨機(jī)誤差項(xiàng)隨機(jī)誤差項(xiàng). 01 1 ,(2.1) kk yxx 2 1.模型定義模型定義 模型模型(2.1)的樣本形式的樣本形式 其中其中i=1,n,表示有表示有n組觀測(cè)值組觀測(cè)值. 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 主編：費(fèi)宇 表表2.1 抽樣調(diào)查得到的抽樣調(diào)查得到的36個(gè)人的數(shù)據(jù)資料個(gè)人的數(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ù) 的最小二乘估計(jì)的最小二乘估計(jì) 的的最小二乘估計(jì)最小二乘估計(jì) 2.模型的參數(shù)估計(jì)和檢驗(yàn)?zāi)Ｐ偷膮?shù)估計(jì)和檢驗(yàn) 6 1 (),(2.5) TT X XX y 2 .(2.7) 1 T e e nk 2 回歸方程的顯著性檢驗(yàn)回歸方程的顯著性檢驗(yàn) 其中其中 是回歸平方和，而是回歸平方和，

4、而 是殘差平方和，拒絕域?yàn)槭菤埐钇椒胶停芙^域?yàn)?2.模型的參數(shù)估計(jì)和檢驗(yàn)?zāi)Ｐ偷膮?shù)估計(jì)和檢驗(yàn) 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ù)估計(jì)和檢驗(yàn)?zāi)Ｐ偷膮?shù)估計(jì)和檢驗(yàn) 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 主編：費(fèi)宇 建立建立y關(guān)于關(guān)于x1、x2、x3和和x4的線性回

5、歸方程的線性回歸方程,并對(duì)方程和回歸并對(duì)方程和回歸 系數(shù)進(jìn)行顯著性檢驗(yàn)系數(shù)進(jìn)行顯著性檢驗(yà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àn)對(duì)應(yīng)的檢驗(yàn)對(duì)應(yīng)的p值

10、則顯示值則顯示:常數(shù)常數(shù) 項(xiàng)和項(xiàng)和x1是顯著的是顯著的,而而x2、x3和和x4不顯著不顯著. 最優(yōu)模型一般滿足最優(yōu)模型一般滿足2個(gè)條件個(gè)條件 （1）模型反映了變量間的真實(shí)關(guān)系）模型反映了變量間的真實(shí)關(guān)系 （2）模型包含的變量盡量）模型包含的變量盡量少少 例例2.1（續(xù)）（續(xù)） 建立建立y關(guān)于關(guān)于x1、x2、x3和和x4的線性回歸方程的線性回歸方程,并對(duì)方程和回并對(duì)方程和回 歸系數(shù)進(jìn)行歸系數(shù)進(jìn)行顯著性檢驗(yàn)顯著性檢驗(yà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.所以這個(gè)所以這個(gè) 回歸過程可以簡(jiǎn)寫為回歸過程可以簡(jiǎn)寫為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.殘差分析和異常點(diǎn)探測(cè)殘差分析和異常點(diǎn)探測(cè) 殘差向量殘差向量e是模型中隨機(jī)誤差項(xiàng)是模型中隨機(jī)誤差項(xiàng) 的估計(jì)，的估計(jì)， 殘差分析可以診斷模型的基本假定是否成殘差分析可以診斷模型的基本假定是否成 立。立。 2.3 2.3 模型診斷模型診斷 2021-5-12 18 eyyyX 例例2.2 計(jì)算計(jì)算例例2.1得到的逐步回歸模型得到的逐步回歸模型lm.ste

18、p的普的普 通殘差和標(biāo)準(zhǔn)化殘差通殘差和標(biāo)準(zhǔn)化殘差,判斷可能存在的異常點(diǎn)判斷可能存在的異常點(diǎn), 畫出相應(yīng)的殘差散點(diǎn)圖畫出相應(yīng)的殘差散點(diǎn)圖,并直觀判斷模型的基并直觀判斷模型的基 本假定是否成立本假定是否成立. 2.3 模型診斷模型診斷 2021-5-12 19 #假設(shè)由例假設(shè)由例2.1已經(jīng)得到逐步回歸模型已經(jīng)得到逐步回歸模型lm.step y.res-residuals (lm.salary) #計(jì)算回歸模型計(jì)算回歸模型lm.salary的普通殘差的普通殘差 y.rst-rstandard(lm.step) #計(jì)算回歸模型計(jì)算回歸模型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) #計(jì)算回歸模型計(jì)算回歸模型lm.step的預(yù)測(cè)值的預(yù)測(cè)值 plot(y.res y.fit) #繪制以普通殘差為縱坐標(biāo)繪制以普通殘差為縱坐標(biāo),預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖 2.3 模型診斷模型診斷 2021-5-12 20 采用采用residuals( ),rstandard( )和和rstudent( )來

20、分別計(jì)算來分別計(jì)算 普通殘差普通殘差,標(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號(hào)和號(hào)和35號(hào)是異常點(diǎn)號(hào)是異常點(diǎn) 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,并計(jì)算更新后模型的標(biāo)準(zhǔn)化并計(jì)算更新后模型的標(biāo)準(zhǔn)化 殘差殘差,畫出相應(yīng)的殘差散點(diǎn)圖畫出相應(yīng)的殘差散點(diǎn)圖,并直觀判斷模型并直觀判斷模型 的基本假定是否成立的基本假定是否成立. 2.3 模型診斷模型診斷 2021-5-12 23 #假設(shè)由例假設(shè)由例2.1已經(jīng)得到逐步回歸模型已經(jīng)得到逐步回歸模型lm.step lm.step_new-update(lm.step,log(.).) #對(duì)模型進(jìn)行對(duì)數(shù)變換對(duì)模型進(jìn)行對(duì)數(shù)變

25、換 y.rst-rstandard(lm.step_new) #計(jì)算新回歸模型計(jì)算新回歸模型lm.step_new的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差 y.fit-predict(lm.step_new) #計(jì)算回歸模型計(jì)算回歸模型lm.step_new的預(yù)測(cè)值的預(yù)測(cè)值 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)圖 (見圖見圖2-2) 采用對(duì)數(shù)變換來解決方差非齊問題采用對(duì)數(shù)變換來解決方差非齊問題 2021-5-12 24 對(duì)數(shù)變換后：只有第對(duì)數(shù)變換后：只有第35號(hào)是異常點(diǎn)號(hào)是異常點(diǎn) 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號(hào)和號(hào)和 35號(hào)觀測(cè)值再建立全變量回歸方程號(hào)觀測(cè)值再建立全變量回歸方程 lm.step-step(lm.salary,direction=both) #用用“一切子集回歸法一切子集回歸法”來進(jìn)行來進(jìn)行 逐步回歸逐步回歸 y.rst-rstandard(lm.step) #計(jì)算回歸模型計(jì)算回歸模型lm.step的標(biāo)準(zhǔn)化殘差的標(biāo)準(zhǔn)化殘差 y.fit-predict(lm.step) #計(jì)算回歸模型計(jì)算回歸模型lm.ste

27、p的預(yù)測(cè)值的預(yù)測(cè)值 plot(y.rst y.fit) #繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo)繪制以標(biāo)準(zhǔn)化殘差為縱坐標(biāo),預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn)預(yù)測(cè)值為橫坐標(biāo)的殘差散點(diǎn) 圖圖 去掉去掉4 4和和3535號(hào)觀測(cè)值號(hào)觀測(cè)值再回歸再回歸 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號(hào)觀測(cè)值號(hào)觀測(cè)值 2.回歸診斷：一般的方法回歸診斷：一般的方法 殘差分析無法分析模型的影響點(diǎn)殘差分析無法分析模型的影響點(diǎn),即探測(cè)哪些點(diǎn)對(duì)模型的即探測(cè)哪些點(diǎn)對(duì)模型的 推斷有重要影響推斷有重要影響,本節(jié)給出本節(jié)給

28、出的的回歸診斷方法回歸診斷方法,可以診斷模型可以診斷模型 的基本假定是否成立的基本假定是否成立,哪些值是異常點(diǎn)哪些值是異常點(diǎn),哪些點(diǎn)是強(qiáng)影響點(diǎn)哪些點(diǎn)是強(qiáng)影響點(diǎn). 在在R中中,函數(shù)函數(shù)plot( )和和influence.measures( )可以用來繪制診斷可以用來繪制診斷 圖和計(jì)算診斷統(tǒng)計(jì)量圖和計(jì)算診斷統(tǒng)計(jì)量 2.3 模型診斷模型診斷 2021-5-12 28 例例2.4 對(duì)例對(duì)例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) #在一個(gè)在一個(gè)22網(wǎng)格中創(chuàng)建網(wǎng)格中創(chuàng)建4個(gè)繪圖區(qū)個(gè)繪圖區(qū) plot(lm.step_new) #繪制模型診斷圖繪制模型診斷圖 influence.measures(lm.step_new) #計(jì)算各個(gè)觀測(cè)值的診計(jì)算各個(gè)觀測(cè)值的診 斷統(tǒng)計(jì)量斷統(tǒng)計(jì)量 運(yùn)行上述程運(yùn)行上述程進(jìn)行進(jìn)行序可得回歸診斷圖序可得回歸診斷圖(圖圖2.4)和如下和如下 36個(gè)觀測(cè)值對(duì)應(yīng)的診斷統(tǒng)計(jì)量的值個(gè)觀測(cè)值對(duì)應(yīng)的診斷統(tǒng)計(jì)量的值. 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號(hào)觀測(cè)值可能是異常點(diǎn)和強(qiáng)影響點(diǎn)號(hào)觀測(cè)值可能是異常點(diǎn)和強(qiáng)影響點(diǎn). 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