多元線性回歸分析

1、多元線性回歸分析直線回歸概念復(fù)習(xí)例:為了研究3歲至8歲男孩身高與年齡的規(guī)律，在某地區(qū)在3歲至8歲男孩中隨機(jī)抽樣，共分6個(gè)年齡層抽樣：3歲，4歲，8歲，每個(gè)層抽10個(gè)男孩，共抽60個(gè)男孩。資料如下：60個(gè)男孩的身高資料如下年齡3歲4歲5歲6歲7歲8歲身高92.596.5106.0115.5125.5121.597.0101.0104.0115.5117.5128.596.0105.5107.0111.5118.0124.096.5102.0109.5110.0117.0125.597.0105.0111.0114.5122.0122.592.099.5107.5112.5119.0123.596

2、.5102.0107.0116.5119.0120.591.0100.0111.5110.0125.5123.096.0106.5103.0114.5120.5124.099.0100.0109.0110.0122.0126.5平均身高95.4101.8107.6113.1120.6124.0圖1 某地男童身高與年齡的散點(diǎn)圖從散點(diǎn)圖上，我們可以發(fā)現(xiàn)樣本點(diǎn)(X,Y)隨機(jī)地出現(xiàn)在一條直線附近，并且從資料背景上考察，同一年齡的兒童身高應(yīng)近似服從一個(gè)正態(tài)分布，而兒童身高的總體均數(shù)應(yīng)隨著年齡增長(zhǎng)而增大，并由每個(gè)年齡的身高樣本均數(shù)與兒童年齡的散點(diǎn)圖可以發(fā)現(xiàn)：這些點(diǎn)非常接近一條直線以及樣本均數(shù)存在抽樣誤差，

3、因此推測(cè)兒童身高的總體均數(shù)與年齡可能呈直線關(guān)系。故假定身高Y在年齡X點(diǎn)上的總體均數(shù)與X呈直線關(guān)系。其中y表示身高，x表示年齡。由于身高的總體均數(shù)與年齡有關(guān)，所以更準(zhǔn)確地標(biāo)記應(yīng)為表示在固定年齡情況下的身高總體均數(shù)。身高的樣本均數(shù)與年齡的散點(diǎn)圖故有理由認(rèn)為身高的總體均數(shù)與年齡的關(guān)系可能是一條直線關(guān)系上述公式稱為直線回歸方程。其中b為回歸系數(shù)（regression coefficient），或稱為斜率（slope）；a稱為常數(shù)項(xiàng)（constant），或稱為截距（intercept）?；貧w系數(shù)b表示x變化一個(gè)單位y平均變化b個(gè)單位。當(dāng)x和y都是隨機(jī)的，x、y間呈正相關(guān)時(shí)b>0，x、y間呈負(fù)相關(guān)時(shí)

4、b<0，x、y間獨(dú)立時(shí)b=0。一般情況而言，參數(shù)a和b是未知的。對(duì)于本例而言，不同民族和不同地區(qū)，a和b往往是不同的，因此需要進(jìn)行估計(jì)的。由于不同年齡的身高實(shí)際觀察值應(yīng)在對(duì)應(yīng)的身高總體均數(shù)附近(即：實(shí)際觀察值與總體均數(shù)之間僅存在個(gè)體變異的差異)，故可以用年齡和實(shí)際身高觀察值的資料對(duì)未知參數(shù)a和b進(jìn)行估計(jì)，一般采用最小二乘法進(jìn)行參數(shù)估計(jì)。我們將借助Stata軟件對(duì)本例資料進(jìn)行直線回歸。數(shù)據(jù)格式xy392.5397.0396.0396.5397.0392.0396.5391.0396.0399.0496.54101.04105.54102.04105.0499.54102.04100.041

5、06.54100.05106.05104.05107.05109.55111.05107.55107.05111.55103.05109.06115.56115.56111.56110.06114.56112.56116.56110.06114.56110.07125.57117.57118.07117.07122.07119.07119.07125.57120.57122.08121.58128.58124.08125.58122.58123.58120.58123.08124.08126.5回歸命令regress y x Source | SS df MS Number of obs =

6、60-+- F( 1, 58) = 777.41 Model | 5997.71571 1 5997.71571 Prob > F = 0.0000 Residual | 447.467619 58 7.71495895 R-squared = 0.9306-+- Adj R-squared = 0.9294 Total | 6445.18333 59 109.240395 Root MSE = 2.7776- y | Coef. Std. Err. t P>|t| 95% Conf. Interval-+- x | 5.854286 .2099654 27.88 0.000 5.

7、433994 6.274577 _cons | 78.18476 1.209202 64.66 0.000 75.76428 80.60524-回歸方程 b=5.854286 ， a= 78.18476se(b)= 0.2099654 回歸系數(shù)檢驗(yàn)：H0：b=0 vs H1:b¹0回歸系數(shù)統(tǒng)計(jì)量t=b/se(b)= 5.854286/ .2099654=27.88，P值<0.001，95%CI of b 為 (5.433994,6.274577)1) 簡(jiǎn)述單因素線性回歸方程y=a+bx在實(shí)際分析中要注意的問題(a) 殘差eiyiabxi，引入回歸模型yi=a+bxi+ei(b)

8、 eiN(0,s)且ei相互獨(dú)立：說明有三個(gè)條件：i) ei服從正態(tài)分布ii) ei相同的方差s2。iii) ei相互獨(dú)立。(c) 不滿足上述3個(gè)條件時(shí)，反映在實(shí)際回歸分析時(shí)，有如下情況：i) 散點(diǎn)在直線一側(cè)較多而且靠直線很近，當(dāng)在直線的另一側(cè)，散點(diǎn)較少，而且離直線較遠(yuǎn)，反映在誤差項(xiàng)e偏態(tài)分布。ii) 散點(diǎn)隨著自變量x增大而離散程度增大或減小(喇叭口狀),反映了誤差項(xiàng)e方差隨著x變而變，即不滿足相同方差(方差齊性)。iii) 隨著xi變化而ei呈某種規(guī)律性的變化。反映e還含有x的信息未利用到，還可以繼續(xù)改進(jìn)回歸模型。問題1：在同一總體中隨機(jī)抽取2個(gè)相同樣本量的樣本，每個(gè)樣本中都含有變量x和y，

9、并以y為因變量和x為自變量，作線性回歸，請(qǐng)問：兩個(gè)樣本作出的回歸方程一樣嗎？它們之間什么關(guān)系？問題2：回歸方程所示的直線與原始數(shù)據(jù)的關(guān)系是什么？1) 不同，它們之間存在抽樣誤差2) 回歸分析統(tǒng)計(jì)背景：對(duì)于固定自變量x，對(duì)y所在的總體進(jìn)行抽樣，得到在固定x情況下，y的樣本值，因此對(duì)于每個(gè)xi，得到對(duì)應(yīng)的抽樣值yi。即：資料為：(x1，y1)，(x2，y2)，(xn，yn)。因此對(duì)于同一個(gè)x值，y所對(duì)應(yīng)的總體均數(shù)相同，不同的x值，y所對(duì)應(yīng)的總體均數(shù)可能不同。如果y的總體均數(shù)值與x的關(guān)系呈直線關(guān)系，則樣本資料(x1，y1)，(x2，y2)，(xn，yn)呈帶狀直線散點(diǎn)圖。由于抽樣資料y=總體均數(shù)抽樣

10、誤差因此如果y的總體均數(shù)值與x呈直線關(guān)系，則抽樣資料當(dāng)，則對(duì)于固定x，而用樣本資料(x1，y1)，(x2，y2)，(xn，yn)所估計(jì)得到的回歸方程是固定x情況下，y的總體均數(shù)與x的線性方程的表達(dá)式。即：b是的樣本估計(jì)值(無偏估計(jì))，a是的樣本估計(jì)值(無偏估計(jì))，是的樣本估計(jì)值。抽樣誤差(估計(jì)值)樣本資料(a+bx) （即： 的估計(jì)值：殘差）所以要求回歸分析的資料，其殘差服從正態(tài)分布，且與x無關(guān)、方差齊性。2) 引入多元線性回歸模型定義(a) 例3-1，研究女中學(xué)生的肺活量與體重和胸圍的關(guān)系，隨機(jī)抽樣了10名女中學(xué)生的體重x1(kg)，胸圍x2(cm)和肺活量y(ml)，資料如表31，試建立一

11、個(gè)因變量為y對(duì)自變量x1,x2的線性回歸方程。(b) 對(duì)于相同的體重x1和胸圍x2，考查女中學(xué)生的肺活量y總是有一定的變異的，但總對(duì)應(yīng)有一個(gè)總體均數(shù)my|X，而且總體均數(shù)my|X可能與體重x1和胸圍x2有關(guān)。x1和x2與總體均數(shù)my|X最簡(jiǎn)單的關(guān)系為線性關(guān)系：i) 同樣的x1和x2，觀察值y與總體均數(shù)my總有一定的隨機(jī)誤差e，即y-my|X=e，因此ii) 若eN(0,s2)分布且獨(dú)立，而觀察值，則稱肺活量y、體重x1和胸圍x2符合線性回歸模型(c) 對(duì)于一般的線性回歸模型定義為：i) 設(shè)有p個(gè)觀察自變量x1，x2，xp ，并用向量X=( x1，x2，xp)，因變量為y，且記y的總體均數(shù)為，隨

12、機(jī)誤差eN(0,s2)且獨(dú)立，則線性回歸模型可以表示為對(duì)于觀察值(y1,X1)，(y2,X2)，(yn,Xn)，其中Xi=(xi1,xi2,，xip)，i=1,2,n。對(duì)應(yīng)的線性回歸模型為且獨(dú)立。在本例中，作線性回歸如下：(介紹一下數(shù)據(jù)結(jié)構(gòu)). regress y x1 x2SourceSSdfMSNumber of obs = 10F( 2, 7) = 6.75回歸平方和回歸均方和Model1895106.552947553.275Prob > F = 0.0232殘差平方和殘差均方和決定系數(shù)Residual982143.457140306.207R-squared = 0.6587校

13、正和決定系數(shù)Adj R-squared = 0.5611Total2877250.009319694.444Root MSE = 374.57總平方和SS總描述樣本量為n10的因變量y總的變異?；貧w平方和SSR描述了樣本量為n時(shí)，由自變量x1,x2變化而引起的因變量y的這部分變異，SSe描述了樣本量為n時(shí)，由隨機(jī)誤差項(xiàng)e所引起的因變量y的一部分變異，因此：總變異自變量引起y的變異隨機(jī)誤差e引起變異對(duì)應(yīng)：SS總SS回歸SS誤差由于SS總，SS回歸和SS誤差均與樣本量n有關(guān)，樣本量n越大，對(duì)應(yīng)變異就越大。所以取平均變異指標(biāo)：均方差MS，回歸系數(shù)回歸系數(shù)標(biāo)準(zhǔn)誤t值P值95可信區(qū)間yCoef.Std.

14、Err.tP>|t|95% Conf. Intervalx1113.998738.311092.9760.02123.40741204.5901x245.4836828.184281.6140.151-21.16155112.1289_cons-5545.8062293.933-2.4180.046-10970.1-121.5156回歸方程 解釋回歸系數(shù)的意義簡(jiǎn)述SST總SSR回歸SSE殘差，自由度df回歸模型中的回歸系數(shù)個(gè)數(shù)(不含常數(shù)項(xiàng))，df殘差=ndf回歸1， 模型的假設(shè)檢驗(yàn)H0：b1=b2=0 vs b1，b2不全為0當(dāng)H0成立時(shí)，F(xiàn)(df回歸,df殘差) 單個(gè)回歸系數(shù)檢驗(yàn)：H0

15、：b0 vs H1：b¹0當(dāng)H0：b0成立時(shí)，簡(jiǎn)述回歸系數(shù)b的95CI 意義與t檢驗(yàn)的對(duì)應(yīng)關(guān)系。(d) 假設(shè)檢驗(yàn)一般情況敘述(e) 決定系數(shù)(f) 復(fù)相關(guān)系數(shù)R(g) H0:b1=b2=br=0 vs b1,b2,br不全為0。當(dāng)H0成立時(shí)m(x1,x2,xp)的估計(jì)及其誤差(STATA命令：predict y1) (STATA命令：predict meansd,stdp)(因?yàn)橛谐闃诱`差)95%CI ，自由度v=n-1-p個(gè)體預(yù)測(cè)值和標(biāo)準(zhǔn)誤(STATA命令：predict y1)線性回歸模型應(yīng)用的條件總結(jié)理論上且獨(dú)立。具體檢查是否復(fù)合線性回歸模型步驟1. 先做線性回歸2. 計(jì)算殘差

16、ei 3. 檢查殘差ei是否服從正態(tài)分布(引起正態(tài)分布)4. 檢查殘差ei的離散程度是否與其它自變量呈某種趨勢(shì)關(guān)系。(要求無任何趨勢(shì)關(guān)系)5. 檢查殘差ei變化是否與其它自變量呈某種對(duì)應(yīng)趨勢(shì)關(guān)系。(要求無任何趨勢(shì)關(guān)系)多元線性回歸常見的應(yīng)用以及應(yīng)用中的問題l 全回歸模型(析因分析)l 多重共線對(duì)分析的影響VIFs (variance inflation factors)l 對(duì)于自變量p個(gè)自變量x1，x2，xp中，以其中一個(gè)xi作為因變量作回歸以及其它p-1個(gè)變量為自變量，得到相應(yīng)的決定系數(shù)Ri。定義xi的膨脹因子l VIFi=1對(duì)應(yīng)說明xi與其它p-1個(gè)自變量無共線。l 當(dāng)對(duì)應(yīng)VIFi>

17、1l 當(dāng)，說明xi與其它p-1個(gè)自變量完全共線，對(duì)應(yīng)VIFi成為無窮大。l 通常認(rèn)為在p個(gè)自變量x1，x2，xp中,最大的VIF>10，則認(rèn)為嚴(yán)重共線，最小二乘估計(jì)受到較嚴(yán)重的影響。l 平均VIF>>1，則認(rèn)為l 尋找影響因變量的主要因素。l 用回歸進(jìn)行兩組或多組的均數(shù)比較并校正混雜因素的影響。全回歸分析舉例例：據(jù)兒童保健部門的考察，4至7歲兒童的身高與年齡近似呈線性關(guān)系，且男女身高也有差異。下列收集了50名男孩和50名女孩的身高，年齡均在4歲至7歲之間。請(qǐng)?jiān)嚱⒒貧w方程描述年齡與身高的關(guān)系(其中sex=1表示男，sex=0表示女)sexagey14.59016.511116

18、.210716.410716.711414.48816.410914.28616.210717.4122159514.18515.610017.51211610617.312014.89316.2105159417.712515.19614.48815.610116.811317.412115.810515.610217.512214.28416.711316.811516.711414.99314.38616.310815.49917.211614.48716.310914.48917.812514.892159514.6901711715.49915.510217.812716.31101

19、7.111904.38707.2114059505.810004.59004.99104.18604.69005.19406.510907.511605.910404.99407.711807.511607.411704.79106.510706.911206.110504.38905.59904.18507.211305.61010610405.49805.19505.610104.79007.912004.79005.19504.99406.410804.38806.210706.8110059404.89405.910406.410704.79307.411606.811005.4990

20、5.49905.19607.311507.8121考慮身高總體均數(shù)為模型為:用擬合上述模型gen sexage=sex*ageregress y age sex sexage- y | Coef. Std. Err. t P>|t| 95% Conf. Interval-+- sex | -9.513794 1.119899 -8.50 0.000 -11.73678 -7.290813 age | 9.075835 .1337354 67.86 0.000 8.810372 9.341298 sexage | 1.929241 .1883106 10.24 0.000 1.555447

21、 2.303035 _cons | 48.97983 .7869668 62.24 0.000 47.41771 50.54194回歸方程為則女孩為身高與年齡的回歸方程為(sex=0)age的回歸系數(shù)的意義為每年身高增長(zhǎng)的速度則男孩為身高與年齡的回歸方程為(sex=1)age的回歸系數(shù)的意義為每年身高增長(zhǎng)的速度因此女孩身高的增長(zhǎng)速度為b2，樣本估計(jì)值為9.075835男孩身高的增長(zhǎng)數(shù)為b2b3，樣本估計(jì)值為11.005076男孩與女孩身高的增長(zhǎng)速度差異為b3，b3>0說明男孩身高增長(zhǎng)速度快，b3<0說明女孩身高增長(zhǎng)速度快，b3說明女孩與男孩的身高增長(zhǎng)速度是一樣的。樣本估計(jì)值為1.9

22、29241>0，P值<0.001。因此男孩身高速度高于女孩，并且差別有統(tǒng)計(jì)學(xué)意義。例：治療缺鐵性貧血100人，隨機(jī)分為2組，給予不同療法治療：經(jīng)過一個(gè)月治療后，治療前后的紅細(xì)胞數(shù)(萬/ml)如下：A組B組治療前y1治療后y2組別group治療前y1治療后y2組別group32533713273480312325133435403313431347368032834113173370316330135137103673801299319035436713363570311325131733803643781305326034536013623820335348131533303293

23、441370394033634913463680293306132434503453581324346036437813623830311325131833803473601329350035036413563780295308135637603693831356378032333613403620385399132234203243381310330031232513573780322336134536503403531340361033034413303510347361135838003613741306329037438913223420327340130432503353491327

24、348036337713533740338350135537603283441346369030331613693900329342132634803173311333355033434613673890334348136338403353481337360033034313683890338353133936103533661337358033234513693900303317135838003693841357378032834313453680治療前治療后第一組335.28±20.840541348.82±21.04678第二組339.98±19.8756

25、23361.14±20.188914考慮以治療前后的改變量為評(píng)價(jià)的效應(yīng)指標(biāo)先不考慮校正基線則可以用成組t檢驗(yàn)進(jìn)行統(tǒng)計(jì)分析gen y=y2-y1ttest y,by(group)結(jié)果如下：Two-sample t test with equal variances- Group | Obs Mean Std. Err. Std. Dev. 95% Conf. Interval-+- 0 | 49 21.16327 .1524933 1.067453 20.85666 21.46987 1 | 49 13.57143 .1271081 .8897565 13.31586 13.827-+

26、-combined | 98 17.36735 .3978661 3.938674 16.57769 18.157-+- diff | 7.591837 .1985212 7.197775 7.985898-Degrees of freedom: 96 Ho: mean(0) - mean(1) = diff = 0 Ha: diff < 0 Ha: diff = 0 Ha: diff > 0 t = 38.2419 t = 38.2419 t = 38.2419 P < t = 1.0000 P > |t| = 0.0000 P > t = 0.0000現(xiàn)用線性

27、回歸完成上述分析設(shè)B組(group=0)受試者的紅細(xì)胞數(shù)改變量的總體均數(shù)為md=a，設(shè)A組(group=1)受試者的紅細(xì)胞數(shù)改變量的總體均數(shù)為md=a+b因此兩組的總體均數(shù)可以表示為md=a+bgroup用線性回歸. regress y group Source | SS df MS Number of obs = 98-+- F( 1, 96) = 1462.45 Model | 1412.08163 1 1412.08163 Prob > F = 0.0000 Residual | 92.6938776 96 .965561224 R-squared = 0.9384-+- Adj

28、R-squared = 0.9378 Total | 1504.77551 97 15.5131496 Root MSE = .98263- y | Coef. Std. Err. t P>|t| 95% Conf. Interval-+- group | -7.591837 .1985212 -38.24 0.000 -7.985898 -7.197775 _cons | 21.16327 .1403757 150.76 0.000 20.88462 21.44191-Ø a的估計(jì)值為21.16327，正是B組的樣本均數(shù)Ø b的估計(jì)值為-7.591837，ab21.

29、16327-7.59183713.571433，正是A組的樣本均數(shù)Ø b的估計(jì)值為兩組樣本均數(shù)的差值，b的檢驗(yàn)統(tǒng)計(jì)量t=-38.24，與t檢驗(yàn)結(jié)果對(duì)應(yīng)，P值也對(duì)應(yīng)。Ø 可以證明：成組t檢驗(yàn)也可以用線性回歸分析進(jìn)行。Ø 從本例中可以發(fā)現(xiàn)回歸系數(shù)b的意義就是兩組總體均數(shù)的差值，其估計(jì)值同樣為兩組樣本均數(shù)的差值。gen y=y2-y1regress y group y1 Source | SS df MS Number of obs = 98-+- F( 2, 95) = 769.69 Model | 1417.30895 2 708.654475 Prob >

30、F = 0.0000 Residual | 87.4665611 95 .920700644 R-squared = 0.9419-+- Adj R-squared = 0.9407 Total | 1504.77551 97 15.5131496 Root MSE = .95953- y | Coef. Std. Err. t P>|t| 95% Conf. Interval-+- group | -7.546723 .194777 -38.75 0.000 -7.933405 -7.160042 y1 | .0114537 .0048069 2.38 0.019 .0019108 .

31、0209966 _cons | 17.27509 1.637541 10.55 0.000 14.02416 20.52602-predict e,residual 計(jì)算殘差值eisktest e 殘差正態(tài)性檢驗(yàn) Skewness/Kurtosis tests for Normality - joint - Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2-+- e | 0.233 0.221 3.00 0.2230gen ee=abs(e) 產(chǎn)生殘差e的絕對(duì)值，放在變量ee(檢驗(yàn)方差齊性：Levens方差檢驗(yàn)) ano

32、va ee group Number of obs = 98 R-squared = 0.0042 Root MSE = .589872 Adj R-squared = -0.0061 Source | Partial SS df MS F Prob > F -+- Model | .141918237 1 .141918237 0.41 0.5246 group | .141918237 1 .141918237 0.41 0.5246 Residual | 33.4030971 96 .347948928 -+- Total | 33.5450153 97 .3458249 a0.1

33、,P值>>a，因此說明兩組殘差的平均幅度差別無統(tǒng)計(jì)意義。說明殘差方差齊性。析因分析舉例例 為了研究A藥和B藥治療患免疫球蛋白偏低的兒童的療效，采用隨機(jī)對(duì)照試驗(yàn)(RCT)和析因分析的研究設(shè)計(jì)方案：第一組：僅是加強(qiáng)營養(yǎng)(作為對(duì)照組)；第二組：加強(qiáng)營養(yǎng)并服用A藥；第三組：加強(qiáng)營養(yǎng)并服用B藥；第四組：加強(qiáng)營養(yǎng)并服用A藥且B藥。每組隨機(jī)收集了25名患者進(jìn)行治療 評(píng)價(jià)藥物療效的指標(biāo)為IgA(mg/dl 血清)并用y表示定義協(xié)變量a=1表示服用A藥，a=0表示未服用A藥；b=1表示服用A藥，b=0表示未服用B藥；yab44.530150051.83540050.726990051.3288900

34、52.397040043.831780051.507170042.270480050.694540055.666340046.708150041.497520049.043630056.360920050.274930055.001340045.985350050.810790046.856950058.116750047.039850044.075180048.836150052.943210055.672110059.512231059.011921066.25341052.143921062.2821062.844391062.608621056.343021053.540261061.

35、200071058.119051064.337591053.661061060.976091053.845451069.766571056.075211056.799821055.129851063.541231058.3591058.466561069.25211061.039271064.124751048.320730159.511840151.868120154.407550149.394110151.232240146.062440150.264680152.238680156.637790161.278080154.817610151.261580162.673860161.264

36、890160.407320150.628640156.615930158.203610155.521910148.207360153.525130146.578140159.373290153.890150170.062421168.334121167.245481168.924531165.7321179.919531165.737481167.133391166.224531171.421361162.81551170.71411172.844261166.696661166.020161169.713731171.573281165.565641175.723291172.7413311

37、68.246631168.38691167.073911174.800671179.4892611對(duì)照組(a=b=0)服A藥組(a=1,b=0)服B藥組(a=0,b=1)服A藥且B藥組(a=b=1)49.76±4.56859.97±4.770 54.17±4.868 69.72±4.309gen ab=a*b 產(chǎn)生交互作用變量. regress y a b ab Source | SS df MS Number of obs = 100-+- F( 3, 96) = 86.66 Model | 5582.20784 3 1860.73595 Prob > F = 0.0000 Residual | 2061.17815 96 21.4706057 R-squared = 0.7303-+- Adj R-squared = 0.7219 Total | 7643.38599 99 77.2059191 Root MSE = 4.6336- y | Coef. Std. Err. t P>|t| 95% Conf. Interval-+- a | 10.20887 1.310591 7.79 0.000 7.607371 12.81038 b | 4.40