SAS講義第三十二課多元線性回歸分析Word版

1/31多元線性回歸分析多元回歸模型表示法通常，回歸模型包括k個變量，即一個因變量和k個自變量（包括常數(shù)項）。由于具有N個方程來概括回歸模型(32.1)模型的相應矩陣方程表示為：嵌入Equation.3錯誤！未定義書簽。(32.2)式中(32.3)其中：Y為因變量觀察的N列向量，X為自變量觀察的N×(k+1)矩陣，symbol98\f"Symbol"\s10為末知參數(shù)的(k+1))列向量，為誤差觀察的N列向量。在矩陣X表達式中，每一個元素Xij都有兩個下標，第一個下標表示相應的列（變量），第二個下標表示相應的行（觀察）。矩陣X的每一列表示相應的給定變量的N次觀察的向量，與截矩有關的所有觀察值都等于1。經(jīng)典的線性回歸模型的假設可以闡述如下：模型形式由(32.1)給定；矩陣X的元素都是確定的，X的秩為(k+1)，且k小于觀察數(shù)N；為正態(tài)分布，E()=0和，式中I為N×N單位矩陣。根據(jù)X的秩為(k+1)的假定，可以保證不會出現(xiàn)共線性。如果出現(xiàn)完全共線性，矩陣X的一列將為其余列的線性組合，而X的秩將小于(k+1))，關于誤差的假設是最有用的假設，因為用它可以保證最小二乘法估計過程的統(tǒng)計性質(zhì)。除了正態(tài)性外，我們還假定每一個誤差項的平均值為0，方差為常數(shù)，以及協(xié)方差為0。假若我們按Y的分布來表示假設(3)，則可寫成下式：(32.4)最小二乘法估計我們的目的是求出一個參數(shù)向量使得殘差平方和最小，即(32.5)式中，(32.6)(32.7)其中表示回歸殘差的N列向量，而表示Y擬合值的N列向量，表示為估計參數(shù)的(k+1)列向量，將式(32.6)和式(32.7)代入式(32.5)，則得：(32.8)為了確定最小二乘法估計量，我們求ESS對進行微分，并使之等于0，即(32.9)所以(32.10)被稱為“交叉乘積矩陣”的矩陣能夠保證逆變換，這是因為我們假設X的秩為(k+1),該假設直接導致了的非奇異性。最小化的二階條件是，是一個正定矩陣。最小二乘法殘差有一個有益的特性，即(32.11)這個結果說明自變量和殘差的交叉乘積的總和為O，這個公式在一些推導中是非常有用的。現(xiàn)在可以考慮最小二乘估計量的性質(zhì)。首先可以證明它們是無偏估計量。因為(32.12)設式中，且是常數(shù)，這樣(32.13)根據(jù)式(32.13)，可以看到，只要遺漏變量都是隨機分布的，與X無關，并且具有0均值，則最小二乘法估計量將是無偏的。(32.14)我們看到，最小二乘法估計量為線性和無偏估計量。事實上，為的最佳線性無偏估計量，也就是說它在全部無偏估計量中方差最小，這就是著名的高斯－馬爾可夫定理。為了證明高斯－馬爾可夫定理，我們需要證明，任何其他線性估計量b的方差比的方差大。請注意=AY。為了不失去一般性，我們可寫成：(32.15)假如b是無偏的，則(32.16)式(32.16)成立的一個必要和充分的條件是，這樣就可以研究矩陣。由于，所以有(32.17)由于 因為，所以，即(32.18)我們可以看出，為一半正定矩陣。該矩陣的二次型為0，只有當（所有元素為0）時才出現(xiàn)。當時，另外的估計量就是普通最小二乘法估計量，這樣，我們的定理就得到證明。的估計和t檢驗為了計算估計參數(shù)的方差-協(xié)方差矩陣，我們需要給出的估計量，該估計量自然選為(32.19)證明為的一個無偏估計量，雖很單調(diào)冗長，但不困難。因此，是Var()的估計。當為已知時，可用正態(tài)分布假設檢驗。當用近似時，我們不得不用t假設檢驗。為此，我們利用以下的統(tǒng)計結果：若已知，則服從分布，具有N－k－1個自由度；嵌入Equation.3錯誤！未定義書簽。服從分布，具有N－k－1個自由度；嵌入Equation.3錯誤！未定義書簽。，當i=0，1，2,…，k時，服從正態(tài)分布，平均值為0，方差為，其中vi為的第i個對角線元素；嵌入Equation.3錯誤！未定義書簽。和相互獨立。由此得出：(32.20)該式為t分布，具有(N-k-1)個自由度。這就使我們能按照與前面所述相同的方式確定各個回歸參數(shù)的置信區(qū)間。假如t值的絕對值相當大，就可以在適當選定的置信水平上否定原假設，參數(shù)的置信區(qū)間可由下式得出：(32.21)其中為與顯著水平有關的t分布臨界值。R2和F檢驗我們可將Y的總變差分成兩部分，一部分代表已說明變差，另一部分代表末說明變差。為了簡化公式推導過程，首先我們假定Y變量具有0平均值，即=0，則有(32.22)由于和，所以(32.23)式中為總平方和，為回歸（已說明）平方和，為殘差（未說明）平方和，歸納成回歸方差分析表，見表32.1所示。表32.1回歸方差分析表變異來源source離差平方和SS自由度df均方MSF統(tǒng)計量FP概率值P回歸RP誤差E總變異T從而，(32.24)若因變量不具有0平均值，我們必須改進一下的定義。這樣，由此可以得出：(32.25)和(32.26)注意到一個數(shù)學上的事實：隨著模型中增添新的變量，必定會增加，從而只要給模型增添越來越多的新因素，就可能使得人為地增大。在一元回歸時已經(jīng)指出較大常指模型與數(shù)據(jù)擬合得較好，在多元回歸時很容易錯誤地去尋找一個極大化的回歸模型。我們應該知道一個好的多元回歸模型，應具有合理個數(shù)的有意義自變量的簡單模型。為了解決這個問題，提出了修正，使得只有當新增變量確實對因變量有所作用時修正才會增加。我們定義為修正的，它是校正擬合優(yōu)度對自由度的依賴關系，如下式如示：(32.27)現(xiàn)在就可以考慮對回歸系數(shù)集的統(tǒng)計檢驗。最通常利用的檢驗是,這個聯(lián)合假設的檢驗。合適的F統(tǒng)計量為：(32.28)為分布，具有k和N－k－1自由度。較大的值，可使我們否定原假設。reg回歸過程在SAS/STAT中有多個進行回歸的過程，如reg、glm等，常用于進行一般線性回歸模型分析的為reg過程。procreg過程Reg過程一般由下列語句控制：procregdata=數(shù)據(jù)集集名</選項列表>;model因變量=自變量名列</選項列表>；var變量列表；outputout=數(shù)據(jù)集名</選項列表>;plot繪圖表達式</選項列表>；print關鍵字列；weight變量；freq變量；by變量；restrict方程1，方程2，…；test方程1，方程2，…；run;其中model語句是必需要有的，其他語句都是可選的。procreg語句中的<選項列表>。outest=SAS數(shù)據(jù)集——將有關模型的參數(shù)估計和選擇的統(tǒng)計量輸出到指定的SAS數(shù)據(jù)集中。outsscp=SAS數(shù)據(jù)集——要求把平方和及叉積矩陣輸出到type=sscp的數(shù)據(jù)集中。all——屏幕輸出所有內(nèi)容。usscp——對用在該過程中的所有變量輸出平方和及叉積矩陣。noprint——不在屏幕輸出任何內(nèi)容。model語句中的<選項列表>。確定變量篩選辦法的選擇項。selection=none|forward|backward|stepwise|maxr|minr|rsquare|cp|adjrsq依次表示全部變量進入法none、前進法forward、后退法backward、逐步篩選法stepwise（前進法與后退法的結合）、最大R2增量法maxr、最小R2增量法minr、R2選擇法rsquare、Mallow'sCp選擇法cp、修正R2選擇法adjrsq。其他選擇項見表3.2所示是可在model語句中選用的其他選項。表32.2model語句中的其他選項acovxpxspecpcorr1slentry＝detailsaiccovbistbpcorr2slstay＝lackfitsbccorrbpcliscorr1start＝collinss1mserclmscorr2best＝collinointss2ssebjpadjrsqinclude＝influencevifseqbdwrmsegmsepstop＝partialtolallpcspnointsigma＝noprintbic其中一些選擇項的意義如下：acov——存在異方差時，輸出參數(shù)估計量的漸近協(xié)方差陣的估計。spec——進行關于方差異性的檢驗。slentry|sle=顯著性水平——規(guī)定入選變量進人方程的顯著性水平。slstay|sls=剔除水平——規(guī)定從方程中剔除變量的顯著性水平。include＝n——強迫前n個自變量進入模型。start＝s——以含有model語句中前3個自變量的模型開始，進行比較、選擇過程（僅用于maxr或minr方法）。stop＝s——當找到最佳的s個變量模型之后，逐步回歸便停止(僅用于maxr或minr方法)。p——要求計算各觀測點上因變量的預測值。r——作殘差分析，同時給出因變量的預測值。cli——給出各自變量x0所對應的因變量y0的95％置信上、下限。clm——給出各自變量所對應的因變量預測值（均數(shù)）Eyi＝μi的95％置信上、下限。noint——指明回歸方程不帶截距項（常數(shù)項）。stb——要求輸出標準回歸系數(shù)。covb——要求輸出回歸系數(shù)估計的協(xié)方差（陣）估計。corrb——要求輸出回歸系數(shù)估計的相關矩陣估計。mse——要求輸出隨機擾動項方差的估計。rmse——要求輸出。collin——在對截距未進行校正的情形下，診斷多重共線性，條件數(shù)越大越可能存在共線性。collinoint——在對截距進行校正的情形下，診斷多重共線性。tol——表示共線性水平的容許值。對于某個變量容許值定義為1－，其中是由這個變量和模型中所有其他回歸變量建立的回歸模型所得到的。tol越小說明其可用別的自變量解釋的部分多，自然就越可能與別的自變量存在共線性關系，tol與vif互為倒數(shù)。vif——輸出變量間相關性的方差膨脹系數(shù)，vif越大，說明由于共線性的存在，使方差變大。influence——要求對異常點進行診斷。對每一觀測點，輸出如下表32.3所示統(tǒng)計量：表32.3診斷異常點的統(tǒng)計量名稱（統(tǒng)計量）含義“異?！钡呐袆e準則Leverage(hi)杠桿率hi，第i次觀測自變量的取值在模型中作用的量度（0≤hi≤1）hi越大，則第i次觀測在模型中的作用就越大Cook’sDCOOKD統(tǒng)計量，對某一觀測點引起回歸影響大小的度量。用于診斷異常點。若D＞50％，則可認為該觀測點對模型的擬合有強的影響covratio協(xié)方差矩陣的行列式之比（去掉某一觀測點后、前對比）若|covratio|≥3(自變量個數(shù)+i),則第i個觀測點值得引起注意defits此值大于2，表明該點影響較大debetas此值大于2，表明該點影響較大i——要求打?。ㄆ渲蠿為設計矩陣）。xpx——輸出模型的叉積矩陣。ss1——要求打印第一類的模型參數(shù)估計的順序平方和。ss2——要求打印第二類的模型參數(shù)估計的偏平方和。all——要求輸出SAS所分析的以下選擇項的特性：xpx，ss1，ss2，stb，covb，corrb，seqb，p，r，cli，clm，spec，acov，tol,pcorr1，pcor，r2，scorr1，scorr2。partial——給出每一回歸變量的偏回歸殘差圖。dw——一階自相關檢驗的Durbin-Watson統(tǒng)計量。其他選擇語句output語句——用于把一些計算結果輸出到指定的數(shù)據(jù)集中。有關的關鍵字及其意義如表32.4所示。表32.4reg過程的output語旬中的關鍵字關鍵字意義關鍵字意義關鍵字意義predicted預測值l95m95％clm下限stdpclm的標準差residual殘差u95m95％clm上限stdr殘差的標準差press殘差/(1–hi)l9595％cli下限stdicli的標準差rstudent刀切殘差u9595％cli上限cookedCookD統(tǒng)計量student學生氏殘差h杠桿點統(tǒng)計量hivar語句——列出叉積矩陣中的變量，僅當具有outsscp=sasdataset這個選擇時才使用。plot語句——繪制兩變量的散點圖。語句格式為：plotx*y/選項。其中x和y變量，可以是原始數(shù)據(jù)集中的變量，也可以是統(tǒng)計量關鍵字。若變量是統(tǒng)計量關鍵字時，需要在其后加上一個小圓點“·”。restrict語句——要求計算線性等式約束的最小二乘估計，其中的方程就是關于回歸系數(shù)（用自變量表示）的等式，方程與方程間用逗號分隔。例如對于模型modely=a1a2b1b2，可以用restricta1+a2=1語句，表示參數(shù)估計是在a1+a2=1的條件下，求最小二乘估計。test語句——要求進行線性等式約束的顯著性檢驗，即Tintner檢驗，其中的方程就是關于回歸系數(shù)（用自變量表示）的等式，方程與方程間用逗號分隔；test語句一般不與restrict語句同用。例如對于模型modely=a1a2b1b2，可以用testa1+a2=1語句，表示在a1+a2=1原假設條件下作F檢驗。交互式語句下面的這部分語句可以用在procreg過程中，但常用在reg過程激活后，以交互方式運行。add變量名列表——向模型中增加變量。delete變量名列表——刪除原擬合模型中的有關變量。refit——重新擬合模型。print——輸出有關模型的相關信息。reg過程其詳細用法可參閱SAS/STAT的用戶手冊。實例分析例32.1表32.5列舉了一個班級的學生情況的調(diào)查數(shù)據(jù)，試分析身高對體重的影響。表32.5bclass記錄數(shù)據(jù)name姓名age年齡sex性別height身高(厘米)weight體重(公斤)name姓名age年齡sex性別height身高(厘米)weight體重(公斤)KATE12女14543.1FREDRICK14男15442.2LOUISE12女14955.8ALFRED14男15744.9JANE12女13533.6HENRY14男15954.0JACLYN12女16265.8LEWIS14男15741.8LILLIE12女12729.1EDWARD14男16750.8TIM12男14738.1CHRIS14男15744.9JAMES12男14958.1JEFFERY14男16951.3ROBERT12男12535.9MARY15女15241.8BARBARA13女14750.8AMY15女15750.8ALICE13女14948.6ROBERT15男16458.1SUSAN13女13730.4WILLIAM15男15950.4JOHN13男15944.5CLAY15男16247.7JOE13男15447.7MARK15男15247.2MICHAEL13男14243.1DANNY15男16248.1DAVID13男14535.9MARTHA16女15950.8JUDY14女14936.8MARIAN16女14752.2ELIZABET14女15241.3PHILLIP16男16758.1LESLIE14女15964.5LINDA17女15252.7CAROL14女15438.1KIRK17男16760.8PATTY14女15238.6LAWRENCE17男17278.1分析和操作步驟過程如下。建立數(shù)據(jù)文件。首先要將上表32.5中的數(shù)據(jù)輸入到SAS數(shù)據(jù)集中，可調(diào)用SAS的數(shù)據(jù)步data過程，建立我們所需的bclass數(shù)據(jù)集。程序如下：datastudy.bclass;inputname$agesex$heightweight;cards;KATE12F14543.1LOUISE12F14955.8………LAWRENCE17M172;run;制作變量的散點圖。建立完SAS數(shù)據(jù)集bclass后，一般需要對數(shù)據(jù)集中要分析的變量weight與height制作散點圖，以便能從圖示中反映學生的身高與體重的關系。一般的處理操作都有菜單操作方法和編程方法2種。如果用菜單操作方法，在SAS/Assist環(huán)境中，從PrimaryMenu主菜單中選擇Graphics/Highresolution/Plots/Simplex*yplot…菜單命令，再選擇Activedataset為study.bclass，Verticalaxis為weight，Horizontalaxis為height，可以在additionaloptions選項菜單中通過LineandSymbol子選項選定所需要的連線類型和點的符號等，最后選擇Locals/Run菜單命令提交運行即可顯示圖形。如果用編程方法，程序如下：goptionsreset=globalgunit=pctcback=whiteborderhtitle=6htext=3ftext=swissbcolors=(back);procgplotdata=study.bclass;plotweight*height;run;運行后，在Graph窗口得到見圖321所示的結果。圖圖51體重與身高(weight與height)的散點圖相關系數(shù)計算。如果用菜單操作方法，可選擇Globals/SAS/Assist/DataAnalysis/Elementary/Correlation命令，再選擇Activedataset為study.bclass，Columnstobecorrelated為weight和height，然后提交運行。直接編寫調(diào)用相關系數(shù)計算的程序為：proccorrdata=study.bclass;varweightheight;run;運行后，在Output窗口得到見表32.6所示的結果。表32.6身高與體重(weight與height)的相關系數(shù)CorrelationAnalysis2'VAR'Variables:WEIGHTHEIGHTCorrelationAnalysis2'VAR'Variables:WEIGHTHEIGHTSimpleStatisticsVariableNMeanStdDevSumMinimumMaximumWEIGHT4047.6625010.07415190729.1000078.10000HEIGHT40153.2500010.475256130125.00000172.00000PearsonCorrelationCoefficients/Prob>|R|underHo:Rho=0/N=40WEIGHTHEIGHTWEIGHT1.000000.708440.00.0001HEIGHT0.708441.000000.00010.0從輸出的表32.6可以看出高與體重之間相關系數(shù)為0.70844?；貧w分析。如果用菜單操作方法，可選擇Globals/SAS/Assist/DataAnalysis/Regression/Linearregression命令，再選擇Activedataset為study.bclass，Dependent為weight，Independent為height，然后提交運行。編程實現(xiàn)回歸方法為：procregdata=study.bclass;modelweight=height/rclmclidw;run;其中，模型參數(shù)r表示要輸出殘差分析，包括因變量的觀察值、由輸入數(shù)據(jù)和估計模型來計算的預測值、殘差值、標準誤差、學生化殘差、COOKD統(tǒng)計量。模型參數(shù)clm表示對每個觀察輸出因變量期望值的95%置信上界和下界，僅考慮到參數(shù)估計的偏差，沒有考慮誤差項的偏差。模型參數(shù)cli表示對因變量的各個預測值輸出95%置信上界和下界，這個置信界反映了誤差的偏差以及參數(shù)估計的偏差。模型參數(shù)dw表示要進行誤差項的對立性檢驗，計算Durbin-Watson統(tǒng)計量。運行后，在Output窗口得到見表32.7所示的結果。Model:MODEL1DependentVariable:WEIGHTAnalysisofVariance(方差分析)SumofMeanSourceModel:MODEL1DependentVariable:WEIGHTAnalysisofVariance(方差分析)SumofMeanSourceDFSquaresSquareFValueProb>FModel11986.484571986.4845738.2870.0001Error381971.5691851.88340CTotal393958.05375RootMSE7.20301R-square0.5019DepMean47.66250AdjR-sq0.4888C.V.15.11254ParameterEstimates(參數(shù)估計)ParameterStandardTforH0:VariableDFEstimateErrorParameter=0Prob>|T|INTERCEP1-56.74857516.91239600-3.3550.0018HEIGHT10.6813120.110107706.1880.0001誤差項的獨立性檢驗Durbin-WatsonD1.471(ForNumberofObs.)401stOrderAutocorrelation0.185置信區(qū)間DepVarPredictStdErrLower95%Upper95%Lower95%Upper95%StdErrObsWEIGHTValuePredictMeanMeanPredictPredictResidualResidual143.100042.04171.45739.092544.990827.164756.91871.05837.054255.800044.76691.23142.274347.259529.973759.560211.03317.097333.600035.22862.31030.552739.904419.915550.5417-1.62866.823……3852.700046.81091.14744.488549.133232.045361.57645.88917.1113960.800057.03051.89553.195360.865841.952972.10823.76956.9494078.100060.43712.35855.663965.210345.094075.780217.66296.806殘差分析StudentCook'sObsResidual-2-1-012D10.150|||0.00021.555||***|0.0363-0.239|||0.00341.728||***|0.0675-0.104|||0.0016-0.749|*||0.01071.879||***|0.053……35-0.110|||0.000361.242||**|0.027370.154|||0.001380.828||*|0.009390.542||*|0.011402.595||*****|0.404SumofResiduals0SumofSquaredResiduals1971.5692PredictedResidSS(Press)2209.7166回歸分析根據(jù)所選擇的模型參數(shù)的輸出，輸出分為若干段，下面逐個段地給以說明：方差分析表提供關于擬合模型的一般信息?？傆^察數(shù)N=40，自變量個數(shù)k=1，回歸模型帶有截距i=1?；貧w模型的離差平方和RSS=1986.48457，自變量的個數(shù)k=1，所以自由度df=k=1，計算公式見(31.29)式。因變量的樣本離差平方和TSS=3958.05375，自由度為df=N－1=40－1=39，計算公式見(31.34)式。誤差項的樣本離差平方和ESS=1971.56918，自由度df=N－k－1=40-1-1=38，計算公式見(31.32)式。注意TSS=RSS+ESS，即3958.05375=1986.48457+1971.56918?；貧w模型的離差平方和平均值MSR=RSS/df=1986.48457/1=1986.48457，誤差項的離差平方和平均值MSE=ESS/df=1971.56918/38=51.88340。在原假設所有自變量的回歸系數(shù)都為0的情況下，本例只有一個自變量height，即H0：，F(xiàn)(1,38)=MSR/MSE=1986.48457/51.88340=38.287，查F分布表，p值為0.0001小于顯著水平0.05，表明可拒絕原假設并有足夠的證據(jù)斷定回歸線的斜率不為零。所以，這一模型擬合數(shù)據(jù)比基線模型好。無偏的誤差估計標準值RootMSE==7.20301，因變量weight平均值DepMean=47.66250，變異系數(shù)（或稱方差系數(shù)）CV=RootMSE/DepMean×100=7.20301/47.66250×100=15.11254，它表示與單位無關的方差。RSquare是0～1之間的值，它表示貢獻給模型而不是貢獻給擬合殘差的總方差的那部分，它也稱為決定系數(shù)或擬合優(yōu)度，用于判斷回歸模型擬合好壞。R2=1－ESS/TSS=RSS/TSS=1986.48457/3958.05375=0.5019，調(diào)整R2=1－ESS/TSS×(N-i)/（N-k-i）=1-1971.56918/3958.05375×39/38=0.4888，R2越是接近1說明模型擬合的越好，等于1則說明完全擬合，沒有任何信息丟失，本例的R2值表明有一半信息丟失沒有被回歸模型表示出來，通常R2應該超過0.7以上才比較好。參數(shù)估計表給出截距和斜率的估計值，方程表明截距的估計值為－56.748575，斜率的估計值為0.681312，計算公式見(31.17)和(31.19式)。估計截距的標準誤差計算公式見(31.37)式，其中，自變量height的平均值=153.25，自變量height的離差平方和=4279.5，估計誤差51.88340，所以估計截距的標準誤差=16.912396，截距等于零的原假設下，計算出的t(38)=－56.748575/16.912396=－3.355，大于此臨界點絕對值出現(xiàn)的概率為0.0018，遠遠地小于5%，有充足的理由否決截距為零的原假設。估計斜率的標準誤差計算公式見(5.1.38)式，估計斜率的標準誤差=0.1101077，斜率等于零的原假設下，計算出的t(38)=0.681312/0.1101077=6.188，大于此臨界點絕對值出現(xiàn)的概率為0.0001，遠遠地小于5%，有充足的理由否決斜率為零的原假設。自由度為38的T分布，95%置信區(qū)間的雙側臨界值為2.0243941，所以截距的95%置信區(qū)間的下界=－56.74857556.748575－2.0243941×16.912396=－90.98593007，上界=－56.74857556.748575＋2.0243941×16.912396=－22.5112，斜率的95%置信區(qū)間的下界=0.681312－2.0243941×0.1101077=0.458410683，上界=0.681312＋2.0243941×0.1101077=0.9042135。置信區(qū)間分析，輸出了weight因變量（DepVar）的40條原始觀察值和回歸模型的預測均值（PredictValue），及預測均值的標準差（StdErrPredict）、預測均值的置信區(qū)間下界（Lower95%Mean）和上界（Upper95%Mean）、預測值的置信區(qū)間下界（Lower95%Predict）和上界（Upper95%Predict）、殘差（Residual）、殘差的標準差（StdErrResidual）。我們以第一條觀察（Obs=1）為例來說明計算過程，已知第一條的觀察=43.1，=145，根據(jù)回歸模型最小二乘法計算出的估計參數(shù)，可以得到預測均值為－56.748575＋0.681312×145=42.0417。第一條觀察的杠桿率計算公式見(31.42)式，=0.040904311，所以預測均值的標準差=1.457。預測均值服從自由度為38的T分布，這樣預測均值的95%置信區(qū)間下界=42.0417－2.0243941×1.457=39.0925，上界=42.0417＋2.0243941×1.457=44.9908。預測值的方差除了要考慮參數(shù)估計的偏差，還要考慮誤差項的偏差，所以要在預測均值的偏差上加上一個誤差項的偏差，計算公式見(31.44)式，預測值的標準差=7.34885394，這樣預測值的95%置信區(qū)間下界=42.0417－2.0243941×7.34885394=27.1647，上界=42.0417＋2.0243941×7.34885394=56.9187，我們從上面的置信區(qū)間計算中可以發(fā)現(xiàn)二個知識點，第一知識個點，預測值的置信區(qū)間要大于預測均值的置信區(qū)間，第二個知識點，越是接近自變量height平均值153.25的height觀察值，它的因變量weight預測均值和預測值的置信區(qū)間越是窄，而越是偏離自變量平均值153.25的height觀察值，它的因變量weight預測均值和預測值的置信區(qū)間越是寬，從圖形上直觀地看置信區(qū)間為中間窄，兩頭形成喇叭口。殘差分析，我們?nèi)匀灰缘谝粭l觀察為例來說明計算過程。殘差=43.1000－42.0417=1.0583。標準殘差的計算公式見(31.46)式，標準殘差=7.054，學生化殘差（StudentResidual）=殘差/標準殘差=1.0583/7.054=0.150。由于學生化殘差服從標準正態(tài)分布，將學生化殘差畫在殘差圖上，我們可以清楚地看到大約68%的學生化殘差值落在一個標準差－1到＋1之間，而大約95%學生化殘差值落在二個標準差－2到＋2之間?；旧险J為模型的誤差項正態(tài)分布及同方差假設，在診斷上沒有太大問題。殘差之和=0，殘差的平方和=1971.5692。COOKD統(tǒng)計量用于預測每個觀察點是否為強影響點或稱異常點，它是通過刪除這個觀察點后重新用最小二乘估計求解參數(shù)值，來分析這個觀察點。觀察點的COOKD統(tǒng)計量小于50%，我們認為不存在異常情況。PRESS統(tǒng)計量是預測殘差的平方和，第i個觀察的殘差定義為，其中為刪除第i個觀察后從余下的組數(shù)據(jù)中重新用最小二乘法求出的參數(shù)估計，計算出的第i個觀察的預測值。第i個觀察的預測殘差為。誤差的獨立性檢驗，它是回歸模型的三大假設之一。我們采用針對殘差一階自相關性進行計算的Durbin-Watson統(tǒng)計量來檢驗，計算公式見(31.48)式，相鄰殘差之差的平方和=2899.603，DW=2899.603/1971.56918=1.471，DW值靠近2說明誤差基本上是獨立的，小于2說明是正相關。殘差一階自相關系數(shù)=0.185，接近0也說明了誤差基本上是獨立的。殘差一階自相關系數(shù)的計算方法與一般的相關系數(shù)計算公式類似，殘差值的第一個序列數(shù)據(jù)為第1個殘差到第39個殘差，第二個序列數(shù)據(jù)為第2個殘差到第40個殘差，簡化公式為第一、二個序列殘差數(shù)據(jù)的平均值為0，標準化時（公式的分母）取殘差值從1個到40個，即。輸出帶有回歸線的散點圖。如果我們需要輸出帶有回歸線的散點圖，菜單操作方法是通過在additionaloptions選項菜單中選擇RegressionPlots/Plotsofdependentbyindependentcolumns命令，重新再提交一次。注意，此時還可以同時選擇輸出殘差圖。程序的方法是在procreg過程里增加plot語句，要注意SAS的關鍵字使用在plot語句中時要加小圓點，這里是預測值p關鍵字，增加的plot語句如下：plotweight*height=’+’p.*height=’*’/overlay;如果我們需要輸出高分辯率的回歸線圖形，可以先在reg過程中將擬合的預測值p輸出到一個SAS數(shù)據(jù)集如bclassg中，再調(diào)用gplot過程繪制圖形。增加的output語句如下：outputout=study.bclassgp=predictl95=clil95u95=cliu95;繪制高分辨率的帶有回歸線的散點圖程序如下：goptionsreset=globalgunit=pctcback=whiteborderhtitle=6htext=3ftext=swissbcolors=(back);procgplotdata=bclassg;plotweight*heightpredict*heightclil95*heightcliu95*height/overlay;symbol1v=plusc=redi=noneh=2.5;symbol2i=splinev=nonec=blue;symbol3i=splinev=nonec=redl=3;symbol4i=splinev=nonec=blackl=3;run;注意，我們也可以用圖形自帶i=rlcli95選項，直接繪制預測值的置信區(qū)間上下界。運行后，在Graph窗口得到見圖322所示結果。圖5圖52帶有回歸線、95%置信線的體重與身高(weight與height)散點圖此外，還可用性別來分組，分別對男生和女生進行回歸分析，分別建立男生和女生的回歸模型。例32.2研究耗氧量模型。這是有關身體適應性測試的例子，肺活量與一些簡單的鍛煉測試數(shù)據(jù)的擬合，目的是為了在鍛煉測試的基礎上而不是在昂貴笨重的氧氣消耗測試的基礎上得到方程來預測適應性。由于回歸是相關的，所以理論上還應該請求共線性診斷。該數(shù)據(jù)名為fitness，這是一個對31位成年人心肺功能的調(diào)查結果，它包含的變量見下表32.8，測試的各項數(shù)據(jù)見下表32.9。表32.8fitness數(shù)據(jù)集的變量名變量名含義age年齡weight體重oxygen耗氧量runtime跑15英哩的時間（分）rstpulse休息時每分鐘心跳次數(shù)runpulse跑步時每分鐘心跳次數(shù)maxpulse每分鐘心跳次數(shù)最大值表32.9fitness數(shù)據(jù)集中的測試數(shù)據(jù)ageweightoxygenruntimerstpulserunpulsemaxpulse4489.4744.60911.37621781824075.0745.31310.07621851854485.8454.2978.65451561684268.1559.5718.17401661723889.0249.8749.22551781804777.4544.81111.63581761764075.9845.68111.95701761804381.1949.09110.85641621704481.4239.44213.08631741763881.8760.0558.63481701864473.0350.54110.13451681684587.6637.38814.03561861924566.4544.75411.12511761764779.1547.27310.60471621645483.1251.85510.33501661704981.4249.1568.95441801855169.6340.83610.95571681725177.9146.67210.00481621684891.6346.77410.25481621644973.3750.38810.08761681685773.3739.40712.63581741765479.3846.08011.17621561655276.3245.4419.63481641665070.8754.6258.92481461555167.2545.11811.08481721725491.6339.20312.88441681725173.7145.79010.47591861885759.0850.5459.93491481554976.3248.6739.40561861884861.2447.92011.50521701765282.7847.46710.5053170172在這個鍛煉測試數(shù)據(jù)里，我們感興趣的是耗氧量是如何依賴于其他變量的。建立數(shù)據(jù)文件。程序如下：datafitness;inputageweightoxygenruntimerstpulserunpulsemaxpulse;cards;4489.4744.60911.37621781824075.0745.31310.0762185185………5282.7847.46710.5053170172;run;制作變量的散點圖。fitness數(shù)據(jù)集中的變量較多，我們需要制作每兩個不同變量oxygen、age、weight、runtime、rstpulse、runpulse和maxpulse之間的所有散點圖，即散點圖矩陣，共有7×6=42個散點圖。我們可以通過在SAS/Insight軟件中繪制散點圖矩陣，操作步驟為：在SAS命令框中鍵入insight后按Enter，在SAS/Insight：Open對話單中，選擇work.fitness數(shù)據(jù)集后單擊Open按鈕，將在屏幕的窗口中顯示當前打開的數(shù)據(jù)集work.fitness內(nèi)容，再選擇菜單上的Analyze/ScatterPlot（YX）命令，在出現(xiàn)的ScatterPlot（YX）對話單中，把fitness數(shù)據(jù)集中的7個變量依上面的次序全部加入的Y軸和X軸的列表框中，最后單擊OK。見圖323所示。圖323fitness變量間的散點圖矩陣散點圖矩陣圖中第一行的六個散點圖分別表示oxygen變量作為y軸,其他六個變量作為x軸的散點圖,第一列的六個散點圖分別表示oxygen變量作為x軸,其他六個變量作為y軸的散點圖。對角線上第一個oxygen小方格里的左下角數(shù)字37.388和右上角數(shù)字60圖323fitness變量間的散點圖矩陣繪制散點圖矩陣圖是為了研究變量間的相關性，從圖323中，我們發(fā)現(xiàn)變量runpulse與maxpulse之間存在有較強的共線性，如果在回歸模型中增加方差膨脹系數(shù)（vif），共線性水平的容許值（tol），條件數(shù)（collin）選項對回歸進行診斷，也會得到相同的結論。另外，我們從圖中還發(fā)現(xiàn)耗氧量oxygen與變量runtime有較強的負相關。從下面的相關系數(shù)計算中也能得到這些相同認識。相關系數(shù)計算。編寫的相關系數(shù)計算程序如下：proccorrdata=fitness;varoxygenageweightruntimerstpulserunpulsemaxpulse;labeloxygen='Oxygenconsumption'age='Ageinyears'weight='weightinkg'runtime='Min.torun1.5miles'rstpulse='Heartratewhileresting'runpulse='Heartratewhilerunning'maxpulse='Maximumheartrate';run;運行后，在Output窗口得到表32.10所示的結果。CorrelationAnalysis7'VAR'Variables:OXYGENAGEWEIGHTRUNTIMERSTPULSERUNPULSEMAXPULSESimpleStatisticsVariableNMeanCorrelationAnalysis7'VAR'Variables:OXYGENAGEWEIGHTRUNTIMERSTPULSERUNPULSEMAXPULSESimpleStatisticsVariableNMeanStdDevSumMinimumMaximumLabelOXYGEN3147.37585.32721468.737.388060.0550OxygenconsumptionAGE3147.67745.21141478.038.000057.0000AgeinyearsWEIGHT3177.44458.32862400.859.080091.6300WeightinkgRUNTIME3110.58611.3874328.28.170014.0300Min.torun1.5milesRSTPULSE3153.74198.29441666.040.000076.0000HeartratewhilerestingRUNPULSE31169.645210.25205259.0146.0000186.0000HeartratewhilerunningMAXPULSE31173.77429.16415387.0155.0000192.0000MaximumheartratePearsonCorrelationCoefficients/Prob>|R|underHo:Rho=0/N=31OXYGENAGEWEIGHTRUNTIMERSTPULSERUNPULSEMAXPULSEOXYGEN1.00000-0.30459-0.16275-0.86219-0.34641-0.39797-0.23674Oxygenconsumption0.00.09570.38170.00010.05630.02660.1997AGE-0.304591.00000-0.233540.18875-0.14157-0.33787-0.43292Ageinyears0.09570.00.20610.30920.44750.06300.0150WEIGHT-0.16275-0.233541.000000.143510.022700.181520.24938Weightinkg0.38170.20610.00.44120.90350.32840.1761RUNTIME-0.862190.188750.143511.000000.400540.313650.22610Min.torun1.5miles0.00010.30920.44120.00.02560.08580.2213RSTPULSE-0.34641-0.141570.022700.400541.000000.317970.25750Heartratewhileresting0.05630.44750.90350.02560.00.08130.1620RUNPULSE-0.39797-0.337870.181520.313650.317971.000000.92975Heartratewhilerunning0.02660.06300.32840.08580.08130.00.0001MAXPULSE-0.23674-0.432920.249380.226100.257500.929751.00000Maximumheartrate0.19970.01500.17610.22130.16200.00010.0回歸分析。編寫的回歸分析程序如下：procregdata=fitness;modeloxygen=agemaxpulserstpulserunpulseruntimeweight/ss1ss2;run;Model語句中增加ss1和ss2選項，回歸模型將計算每個自變量兩種不同類型的平方和，其中，ss1是按model語句中自變量的排列順序依此計算每個自變量的平方和，也稱為第一類平方和或稱順序平方和，ss2是把model語句中每個自變量排到變量列表的最后，所計算的第一類平方和，也稱為第二類平方和或稱部分平方和。通過分析每個自變量的這兩類平方和，能知道回歸模型總的平方和的構成和各個自變量所貢獻的平方和，進而能知道哪些自變量是最重要的回歸變量，哪些回歸變量可能是無關緊要的，配合參數(shù)估計的t檢驗，最終為縮減回歸變量提供依據(jù)，到達簡化模型的目的。運行后，在Output窗口得到表32.11所示的結果。Model:MODEL1DependentVariable:OXYGENAnalysisofVarianceModel:MODEL1DependentVariable:OXYGENAnalysisofVarianceSumofMeanSourceDFSquaresSquareFValueProb>FModel6721.97421120.3290422.3160.0001Error24129.407335.39197CTotal30851.38154RootMSE2.32206R-square0.8480DepMean47.37581AdjR-sq0.8100C.V.4.90137ParameterEstimatesParameterStandardTforH0:VariableDFEstimateErrorParameter=0Prob>|T|TypeISSTypeIISSINTERCEP1102.23833912.453047198.2100.000169578363.432659AGE1-0.2199160.09959154-2.2080.037078.98822726.291488MAXPULSE10.3047350.137224722.2210.0361142.35542626.590540RSTPULSE1-0.0008440.05863130-0.0140.988682.4478650.001118RUNPULSE1-0.3731640.12068038-3.0920.005098.36406551.555411RUNTIME1-2.6805160.37488355-7.1500.0001310.368687275.671437WEIGHT1-0.0723800.05467334-1.3240.19809.4499429.449942這一輸出與簡單線性回歸的輸出是完全相仿的，只是進入回歸的自變量有6個。從參數(shù)的估計值容易得到擬合的回歸為：oxygen=102.238339－0．219916age+0．304735maxpulse－0.000844rstpulse－0.373164runpulse－2.680516runtime－0.072380weight方差分析類似于例30.1，多元線性回歸時，F(xiàn)Value是除了截距外其他所有參數(shù)都為零的假設檢驗。RSquare同時還是多重相關系數(shù)的平方，換言之，它是自變量與預測值之間的相關系數(shù)平方。我們這里不再重復分析和計算。多元線性回歸模型的第一個重要問題是模型的簡化，也就是說如何正確地縮減自變量到達最優(yōu)的簡化模型。由于自變量除了對總方差貢獻大小外，還存在著自變量間的相關性，刪除哪一個自變量并不是一個簡單問題。我們在這里重點分析參數(shù)估計表中的TypeISS和TypeIISS平方和，配合分析參數(shù)估計的t檢驗p值，對縮減自變量的依據(jù)有所掌握。類型1的平方和（TypeISS）被稱為序列平方和，它是按model語句中自變量的先后序列，每加入一個自變量到這個序列中，回歸模型因此而增加的平方和，把這個增加的平方和算為這個自變量的類型1的平方和。截距INTERCEP的TypeISS為，稱為修正均值，31×47.375812=69578。age的類型1的平方和等于回歸模型只包括截距和age時的RSS值（78.988227），maxpulse的類型1的平方和等于回歸模型只包括截距、age、maxpulse時的RSS值（221.344）減回歸模型只包括截距和age時的RSS值（78.988227），即221.344－78.988227=142.355426，實質(zhì)上是指回歸模型在前面模型的基礎上因為增加了自變量maxpulse而增加的RSS，用表達式表示為TypeISS（maxpulse）=RSS（modeloxygen=agemaxpulse）－RSS（modeloxygen=age）以此類推，weight的類型1的平方和為TypeISS（weight）=RSS（modeloxygen=agemaxpulserstpulserunpulseruntimeweight）－RSS（modeloxygen=agemaxpulserstpulserunpulseruntime）所以有恆等式TypeISS（age）+TypeISS（maxpulse）+TypeISS（rstpulse）+TypeISS（runpulse）+TypeISS（runtime）+TypeISS（weight）=RSS（modeloxygen=agemaxpulserstpulserunpulseruntimeweight）即721.97421=78.988227+142.355426+82.447865+98.364065+310.368687+9.449942說明回歸模型的平方和可以分解成回歸模型中各自變量的類型1平方和。類型2的平方和（TypeIISS）通常被稱為部分平方和，對于一個給定自變量的類型2平方和等于把這個自變量放在model語句等號右邊的最后，所求得的類型1平方和。例如TypeIISS（age）=RSS（modeloxygen=maxpulserstpulserunpulseruntimeweightage）－RSS（modeloxygen=maxpulserstpulserunpulseruntimeweight）所以model語句最后一個自變量的類型1平方和與類型2平方和總是相等的，例如，表32.11中的最后一個變量weight，TypeISS（weight）=TypeIISS（weight）=9.449942。從類型2平方和的計算可知，自變量的類型2平方和不依賴于在model語句中的列表次序，而類型1平方和卻與自變量在model語句中的列表次序密切相關。每個自變量的類型2平方和總是小于等于類型1平方和，所以所有自變量的類型2平方和之和小于回歸模型的平方和，原因是一個自變量先進入模型，它的貢獻總是大于最后進入模型的自變量，因為變量間有相關性，所以先進入模型的自變量總或多或少把后進入模型相關變量的貢獻也算在自己頭上一點。除非每個自變量之間是不相關的，類型1平方和與類型2平方和才是等價的。判斷回歸模型是否還能縮減自變量，可以通過這兩類平方和構造F檢驗來比較確定。首先，注意到每個自變量的平方和只有一個自由度，所以平方和也就是均方和，除以均方誤差MSE構造了模型中包含這個自變量的F統(tǒng)計量，原假設為包含這個自變量模型與不包含這個自變量的模型是一樣好。這樣，例如，maxpulse變量的類型1平方和的F統(tǒng)計量為相應p=0.000017，拒絕原假設，接受回歸模型中包含age和maxpulse比只包含age要好。變量的類型1平方和的F統(tǒng)計量只有weight的p值是不顯著，不能拒絕回歸模型中包含age、maxpulse、rstpulse、runpulse、runtime、weight與包含age、maxpulse、rstpulse、runpulse、runtime一樣好。類似地我們可以計算maxpulse變量的類型2平方和的F統(tǒng)計量為相應p=0.03434219，拒絕原假設，接受回歸模型中包含age、maxpulse、rstpulse、runpulse、runtime、weight比包含age、rstpulse、runpulse、runtime、weight要好。類型2平方和的F統(tǒng)計量只有rstpulse和weight不顯著。兩類F統(tǒng)計量的p值差異最大的是rstpulse自變量，這是我們首先要考慮刪除的。實際上，自變量類型2平方和的F檢驗等同于這個自變量的參數(shù)t檢驗，因為F值等于t值的平方，例如，。其實，有時直觀地分析兩類平方和的數(shù)值大小就能得到基本的判斷，runtime自變量的兩類平方和都是最大的且占的比例很大，說明是回歸模型中第一重要的自變量。而rstpulse自變量在第一類平方和中有比較大的數(shù)值卻在第二類平方和中是最小的，這是rstpulse自變量應該被考慮第一個刪除的主要原因。剔除不顯著的回歸變量。從參數(shù)估計表檢驗部分也可以看出，對自變量rstpulse和weight的回歸系數(shù)的t假設檢驗，不能拒絕它們?yōu)?的原假設。當然在這里必須小心地看待這些檢驗，因為它們都是在其他自變量都加入回歸的前提下進行顯著性檢驗的，完全可能因為自變量間存在較強的相關而掩蓋他們對回歸的貢獻。所以在剔除不顯著的回歸變量時必須逐個進行。另外，從自變量rstpulse和weight的回歸系數(shù)接近于0，也提示我們應先考慮刪除哪一個變量。因為過程reg具有連續(xù)的交互功能，在執(zhí)行了提交的部分語句后，仍可繼續(xù)提交語句讓它執(zhí)行，直至提交quit語句或因執(zhí)行其他過程而終止。根據(jù)上面刪除自變量的分析，我們先從已加入回歸模型的自變量中剔除rstpulse。可直接提交如下的程序：deleterstpulse;print;run;運行后，得到表32.12所示的結果。Model:MODEL1DependentVariable:OXYGENAnalysisofVarianceModel:MODEL1DependentVariable:OXYGENAnalysisofVarianceSumofMeanSourceDFSquaresSquareFValueProb>FModel5721.97309144.3946227.8950.0001Error25129.408455.17634CTotal30851.38154RootMSE2.27516R-square0.8480DepMean47.37581AdjR-sq0.8176C.V.4.80236ParameterEstimatesParameterStandardTforH0:VariableDFEstimateErrorParameter=0Prob>|T|TypeISSTypeIISSINTERCEP1102.20427511.979289728.5320.000169578376.789349AGE1-0.2196210.09550245-2.3000.030178.98822727.374291MAXPULSE10.3049080.133936422.2770.0316142.35542626.826403RUNPULSE1-0.3734010.11714109-3.1880.0038138.17217952.596237RUNTIME1-2.682523