SAS讲义第三十二课多元线性回归分析Word版

1/31多元线性回归分析多元回归模型表示法通常，回归模型包括k个变量，即一个因变量和k个自变量（包括常数项）。由于具有N个方程来概括回归模型(32.1)模型的相应矩阵方程表示为：嵌入Equation.3错误！未定义书签。(32.2)式中(32.3)其中：Y为因变量观察的N列向量，X为自变量观察的N×(k+1)矩阵，symbol98\f"Symbol"\s10为末知参数的(k+1))列向量，为误差观察的N列向量。在矩阵X表达式中，每一个元素Xij都有两个下标，第一个下标表示相应的列（变量），第二个下标表示相应的行（观察）。矩阵X的每一列表示相应的给定变量的N次观察的向量，与截矩有关的所有观察值都等于1。经典的线性回归模型的假设可以阐述如下：模型形式由(32.1)给定；矩阵X的元素都是确定的，X的秩为(k+1)，且k小于观察数N；为正态分布，E()=0和，式中I为N×N单位矩阵。根据X的秩为(k+1)的假定，可以保证不会出现共线性。如果出现完全共线性，矩阵X的一列将为其余列的线性组合，而X的秩将小于(k+1))，关于误差的假设是最有用的假设，因为用它可以保证最小二乘法估计过程的统计性质。除了正态性外，我们还假定每一个误差项的平均值为0，方差为常数，以及协方差为0。假若我们按Y的分布来表示假设(3)，则可写成下式：(32.4)最小二乘法估计我们的目的是求出一个参数向量使得残差平方和最小，即(32.5)式中，(32.6)(32.7)其中表示回归残差的N列向量，而表示Y拟合值的N列向量，表示为估计参数的(k+1)列向量，将式(32.6)和式(32.7)代入式(32.5)，则得：(32.8)为了确定最小二乘法估计量，我们求ESS对进行微分，并使之等于0，即(32.9)所以(32.10)被称为“交叉乘积矩阵”的矩阵能够保证逆变换，这是因为我们假设X的秩为(k+1),该假设直接导致了的非奇异性。最小化的二阶条件是，是一个正定矩阵。最小二乘法残差有一个有益的特性，即(32.11)这个结果说明自变量和残差的交叉乘积的总和为O，这个公式在一些推导中是非常有用的。现在可以考虑最小二乘估计量的性质。首先可以证明它们是无偏估计量。因为(32.12)设式中，且是常数，这样(32.13)根据式(32.13)，可以看到，只要遗漏变量都是随机分布的，与X无关，并且具有0均值，则最小二乘法估计量将是无偏的。(32.14)我们看到，最小二乘法估计量为线性和无偏估计量。事实上，为的最佳线性无偏估计量，也就是说它在全部无偏估计量中方差最小，这就是著名的高斯－马尔可夫定理。为了证明高斯－马尔可夫定理，我们需要证明，任何其他线性估计量b的方差比的方差大。请注意=AY。为了不失去一般性，我们可写成：(32.15)假如b是无偏的，则(32.16)式(32.16)成立的一个必要和充分的条件是，这样就可以研究矩阵。由于，所以有(32.17)由于 因为，所以，即(32.18)我们可以看出，为一半正定矩阵。该矩阵的二次型为0，只有当（所有元素为0）时才出现。当时，另外的估计量就是普通最小二乘法估计量，这样，我们的定理就得到证明。的估计和t检验为了计算估计参数的方差-协方差矩阵，我们需要给出的估计量，该估计量自然选为(32.19)证明为的一个无偏估计量，虽很单调冗长，但不困难。因此，是Var()的估计。当为已知时，可用正态分布假设检验。当用近似时，我们不得不用t假设检验。为此，我们利用以下的统计结果：若已知，则服从分布，具有N－k－1个自由度；嵌入Equation.3错误！未定义书签。服从分布，具有N－k－1个自由度；嵌入Equation.3错误！未定义书签。，当i=0，1，2,…，k时，服从正态分布，平均值为0，方差为，其中vi为的第i个对角线元素；嵌入Equation.3错误！未定义书签。和相互独立。由此得出：(32.20)该式为t分布，具有(N-k-1)个自由度。这就使我们能按照与前面所述相同的方式确定各个回归参数的置信区间。假如t值的绝对值相当大，就可以在适当选定的置信水平上否定原假设，参数的置信区间可由下式得出：(32.21)其中为与显著水平有关的t分布临界值。R2和F检验我们可将Y的总变差分成两部分，一部分代表已说明变差，另一部分代表末说明变差。为了简化公式推导过程，首先我们假定Y变量具有0平均值，即=0，则有(32.22)由于和，所以(32.23)式中为总平方和，为回归（已说明）平方和，为残差（未说明）平方和，归纳成回归方差分析表，见表32.1所示。表32.1回归方差分析表变异来源source离差平方和SS自由度df均方MSF统计量FP概率值P回归RP误差E总变异T从而，(32.24)若因变量不具有0平均值，我们必须改进一下的定义。这样，由此可以得出：(32.25)和(32.26)注意到一个数学上的事实：随着模型中增添新的变量，必定会增加，从而只要给模型增添越来越多的新因素，就可能使得人为地增大。在一元回归时已经指出较大常指模型与数据拟合得较好，在多元回归时很容易错误地去寻找一个极大化的回归模型。我们应该知道一个好的多元回归模型，应具有合理个数的有意义自变量的简单模型。为了解决这个问题，提出了修正，使得只有当新增变量确实对因变量有所作用时修正才会增加。我们定义为修正的，它是校正拟合优度对自由度的依赖关系，如下式如示：(32.27)现在就可以考虑对回归系数集的统计检验。最通常利用的检验是,这个联合假设的检验。合适的F统计量为：(32.28)为分布，具有k和N－k－1自由度。较大的值，可使我们否定原假设。reg回归过程在SAS/STAT中有多个进行回归的过程，如reg、glm等，常用于进行一般线性回归模型分析的为reg过程。procreg过程Reg过程一般由下列语句控制：procregdata=数据集集名</选项列表>;model因变量=自变量名列</选项列表>；var变量列表；outputout=数据集名</选项列表>;plot绘图表达式</选项列表>；print关键字列；weight变量；freq变量；by变量；restrict方程1，方程2，…；test方程1，方程2，…；run;其中model语句是必需要有的，其他语句都是可选的。procreg语句中的<选项列表>。outest=SAS数据集——将有关模型的参数估计和选择的统计量输出到指定的SAS数据集中。outsscp=SAS数据集——要求把平方和及叉积矩阵输出到type=sscp的数据集中。all——屏幕输出所有内容。usscp——对用在该过程中的所有变量输出平方和及叉积矩阵。noprint——不在屏幕输出任何内容。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——存在异方差时，输出参数估计量的渐近协方差阵的估计。spec——进行关于方差异性的检验。slentry|sle=显著性水平——规定入选变量进人方程的显著性水平。slstay|sls=剔除水平——规定从方程中剔除变量的显著性水平。include＝n——强迫前n个自变量进入模型。start＝s——以含有model语句中前3个自变量的模型开始，进行比较、选择过程（仅用于maxr或minr方法）。stop＝s——当找到最佳的s个变量模型之后，逐步回归便停止(仅用于maxr或minr方法)。p——要求计算各观测点上因变量的预测值。r——作残差分析，同时给出因变量的预测值。cli——给出各自变量x0所对应的因变量y0的95％置信上、下限。clm——给出各自变量所对应的因变量预测值（均数）Eyi＝μi的95％置信上、下限。noint——指明回归方程不带截距项（常数项）。stb——要求输出标准回归系数。covb——要求输出回归系数估计的协方差（阵）估计。corrb——要求输出回归系数估计的相关矩阵估计。mse——要求输出随机扰动项方差的估计。rmse——要求输出。collin——在对截距未进行校正的情形下，诊断多重共线性，条件数越大越可能存在共线性。collinoint——在对截距进行校正的情形下，诊断多重共线性。tol——表示共线性水平的容许值。对于某个变量容许值定义为1－，其中是由这个变量和模型中所有其他回归变量建立的回归模型所得到的。tol越小说明其可用别的自变量解释的部分多，自然就越可能与别的自变量存在共线性关系，tol与vif互为倒数。vif——输出变量间相关性的方差膨胀系数，vif越大，说明由于共线性的存在，使方差变大。influence——要求对异常点进行诊断。对每一观测点，输出如下表32.3所示统计量：表32.3诊断异常点的统计量名称（统计量）含义“异常”的判别准则Leverage(hi)杠杆率hi，第i次观测自变量的取值在模型中作用的量度（0≤hi≤1）hi越大，则第i次观测在模型中的作用就越大Cook’sDCOOKD统计量，对某一观测点引起回归影响大小的度量。用于诊断异常点。若D＞50％，则可认为该观测点对模型的拟合有强的影响covratio协方差矩阵的行列式之比（去掉某一观测点后、前对比）若|covratio|≥3(自变量个数+i),则第i个观测点值得引起注意defits此值大于2，表明该点影响较大debetas此值大于2，表明该点影响较大i——要求打印（其中X为设计矩阵）。xpx——输出模型的叉积矩阵。ss1——要求打印第一类的模型参数估计的顺序平方和。ss2——要求打印第二类的模型参数估计的偏平方和。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统计量。其他选择语句output语句——用于把一些计算结果输出到指定的数据集中。有关的关键字及其意义如表32.4所示。表32.4reg过程的output语旬中的关键字关键字意义关键字意义关键字意义predicted预测值l95m95％clm下限stdpclm的标准差residual残差u95m95％clm上限stdr残差的标准差press残差/(1–hi)l9595％cli下限stdicli的标准差rstudent刀切残差u9595％cli上限cookedCookD统计量student学生氏残差h杠杆点统计量hivar语句——列出叉积矩阵中的变量，仅当具有outsscp=sasdataset这个选择时才使用。plot语句——绘制两变量的散点图。语句格式为：plotx*y/选项。其中x和y变量，可以是原始数据集中的变量，也可以是统计量关键字。若变量是统计量关键字时，需要在其后加上一个小圆点“·”。restrict语句——要求计算线性等式约束的最小二乘估计，其中的方程就是关于回归系数（用自变量表示）的等式，方程与方程间用逗号分隔。例如对于模型modely=a1a2b1b2，可以用restricta1+a2=1语句，表示参数估计是在a1+a2=1的条件下，求最小二乘估计。test语句——要求进行线性等式约束的显著性检验，即Tintner检验，其中的方程就是关于回归系数（用自变量表示）的等式，方程与方程间用逗号分隔；test语句一般不与restrict语句同用。例如对于模型modely=a1a2b1b2，可以用testa1+a2=1语句，表示在a1+a2=1原假设条件下作F检验。交互式语句下面的这部分语句可以用在procreg过程中，但常用在reg过程激活后，以交互方式运行。add变量名列表——向模型中增加变量。delete变量名列表——删除原拟合模型中的有关变量。refit——重新拟合模型。print——输出有关模型的相关信息。reg过程其详细用法可参阅SAS/STAT的用户手册。实例分析例32.1表32.5列举了一个班级的学生情况的调查数据，试分析身高对体重的影响。表32.5bclass记录数据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分析和操作步骤过程如下。建立数据文件。首先要将上表32.5中的数据输入到SAS数据集中，可调用SAS的数据步data过程，建立我们所需的bclass数据集。程序如下：datastudy.bclass;inputname$agesex$heightweight;cards;KATE12F14543.1LOUISE12F14955.8………LAWRENCE17M172;run;制作变量的散点图。建立完SAS数据集bclass后，一般需要对数据集中要分析的变量weight与height制作散点图，以便能从图示中反映学生的身高与体重的关系。一般的处理操作都有菜单操作方法和编程方法2种。如果用菜单操作方法，在SAS/Assist环境中，从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)的散点图相关系数计算。如果用菜单操作方法，可选择Globals/SAS/Assist/DataAnalysis/Elementary/Correlation命令，再选择Activedataset为study.bclass，Columnstobecorrelated为weight和height，然后提交运行。直接编写调用相关系数计算的程序为：proccorrdata=study.bclass;varweightheight;run;运行后，在Output窗口得到见表32.6所示的结果。表32.6身高与体重(weight与height)的相关系数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可以看出高与体重之间相关系数为0.70844。回归分析。如果用菜单操作方法，可选择Globals/SAS/Assist/DataAnalysis/Regression/Linearregression命令，再选择Activedataset为study.bclass，Dependent为weight，Independent为height，然后提交运行。编程实现回归方法为：procregdata=study.bclass;modelweight=height/rclmclidw;run;其中，模型参数r表示要输出残差分析，包括因变量的观察值、由输入数据和估计模型来计算的预测值、残差值、标准误差、学生化残差、COOKD统计量。模型参数clm表示对每个观察输出因变量期望值的95%置信上界和下界，仅考虑到参数估计的偏差，没有考虑误差项的偏差。模型参数cli表示对因变量的各个预测值输出95%置信上界和下界，这个置信界反映了误差的偏差以及参数估计的偏差。模型参数dw表示要进行误差项的对立性检验，计算Durbin-Watson统计量。运行后，在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(参数估计)ParameterStandardTforH0:VariableDFEstimateErrorParameter=0Prob>|T|INTERCEP1-56.74857516.91239600-3.3550.0018HEIGHT10.6813120.110107706.1880.0001误差项的独立性检验Durbin-WatsonD1.471(ForNumberofObs.)401stOrderAutocorrelation0.185置信区间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回归分析根据所选择的模型参数的输出，输出分为若干段，下面逐个段地给以说明：方差分析表提供关于拟合模型的一般信息。总观察数N=40，自变量个数k=1，回归模型带有截距i=1。回归模型的离差平方和RSS=1986.48457，自变量的个数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。回归模型的离差平方和平均值MSR=RSS/df=1986.48457/1=1986.48457，误差项的离差平方和平均值MSE=ESS/df=1971.56918/38=51.88340。在原假设所有自变量的回归系数都为0的情况下，本例只有一个自变量height，即H0：，F(1,38)=MSR/MSE=1986.48457/51.88340=38.287，查F分布表，p值为0.0001小于显著水平0.05，表明可拒绝原假设并有足够的证据断定回归线的斜率不为零。所以，这一模型拟合数据比基线模型好。无偏的误差估计标准值RootMSE==7.20301，因变量weight平均值DepMean=47.66250，变异系数（或称方差系数）CV=RootMSE/DepMean×100=7.20301/47.66250×100=15.11254，它表示与单位无关的方差。RSquare是0～1之间的值，它表示贡献给模型而不是贡献给拟合残差的总方差的那部分，它也称为决定系数或拟合优度，用于判断回归模型拟合好坏。R2=1－ESS/TSS=RSS/TSS=1986.48457/3958.05375=0.5019，调整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以上才比较好。参数估计表给出截距和斜率的估计值，方程表明截距的估计值为－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，大于此临界点绝对值出现的概率为0.0018，远远地小于5%，有充足的理由否决截距为零的原假设。估计斜率的标准误差计算公式见(5.1.38)式，估计斜率的标准误差=0.1101077，斜率等于零的原假设下，计算出的t(38)=0.681312/0.1101077=6.188，大于此临界点绝对值出现的概率为0.0001，远远地小于5%，有充足的理由否决斜率为零的原假设。自由度为38的T分布，95%置信区间的双侧临界值为2.0243941，所以截距的95%置信区间的下界=－56.74857556.748575－2.0243941×16.912396=－90.98593007，上界=－56.74857556.748575＋2.0243941×16.912396=－22.5112，斜率的95%置信区间的下界=0.681312－2.0243941×0.1101077=0.458410683，上界=0.681312＋2.0243941×0.1101077=0.9042135。置信区间分析，输出了weight因变量（DepVar）的40条原始观察值和回归模型的预测均值（PredictValue），及预测均值的标准差（StdErrPredict）、预测均值的置信区间下界（Lower95%Mean）和上界（Upper95%Mean）、预测值的置信区间下界（Lower95%Predict）和上界（Upper95%Predict）、残差（Residual）、残差的标准差（StdErrResidual）。我们以第一条观察（Obs=1）为例来说明计算过程，已知第一条的观察=43.1，=145，根据回归模型最小二乘法计算出的估计参数，可以得到预测均值为－56.748575＋0.681312×145=42.0417。第一条观察的杠杆率计算公式见(31.42)式，=0.040904311，所以预测均值的标准差=1.457。预测均值服从自由度为38的T分布，这样预测均值的95%置信区间下界=42.0417－2.0243941×1.457=39.0925，上界=42.0417＋2.0243941×1.457=44.9908。预测值的方差除了要考虑参数估计的偏差，还要考虑误差项的偏差，所以要在预测均值的偏差上加上一个误差项的偏差，计算公式见(31.44)式，预测值的标准差=7.34885394，这样预测值的95%置信区间下界=42.0417－2.0243941×7.34885394=27.1647，上界=42.0417＋2.0243941×7.34885394=56.9187，我们从上面的置信区间计算中可以发现二个知识点，第一知识个点，预测值的置信区间要大于预测均值的置信区间，第二个知识点，越是接近自变量height平均值153.25的height观察值，它的因变量weight预测均值和预测值的置信区间越是窄，而越是偏离自变量平均值153.25的height观察值，它的因变量weight预测均值和预测值的置信区间越是宽，从图形上直观地看置信区间为中间窄，两头形成喇叭口。残差分析，我们仍然以第一条观察为例来说明计算过程。残差=43.1000－42.0417=1.0583。标准残差的计算公式见(31.46)式，标准残差=7.054，学生化残差（StudentResidual）=残差/标准残差=1.0583/7.054=0.150。由于学生化残差服从标准正态分布，将学生化残差画在残差图上，我们可以清楚地看到大约68%的学生化残差值落在一个标准差－1到＋1之间，而大约95%学生化残差值落在二个标准差－2到＋2之间。基本上认为模型的误差项正态分布及同方差假设，在诊断上没有太大问题。残差之和=0，残差的平方和=1971.5692。COOKD统计量用于预测每个观察点是否为强影响点或称异常点，它是通过删除这个观察点后重新用最小二乘估计求解参数值，来分析这个观察点。观察点的COOKD统计量小于50%，我们认为不存在异常情况。PRESS统计量是预测残差的平方和，第i个观察的残差定义为，其中为删除第i个观察后从余下的组数据中重新用最小二乘法求出的参数估计，计算出的第i个观察的预测值。第i个观察的预测残差为。误差的独立性检验，它是回归模型的三大假设之一。我们采用针对残差一阶自相关性进行计算的Durbin-Watson统计量来检验，计算公式见(31.48)式，相邻残差之差的平方和=2899.603，DW=2899.603/1971.56918=1.471，DW值靠近2说明误差基本上是独立的，小于2说明是正相关。残差一阶自相关系数=0.185，接近0也说明了误差基本上是独立的。残差一阶自相关系数的计算方法与一般的相关系数计算公式类似，残差值的第一个序列数据为第1个残差到第39个残差，第二个序列数据为第2个残差到第40个残差，简化公式为第一、二个序列残差数据的平均值为0，标准化时（公式的分母）取残差值从1个到40个，即。输出带有回归线的散点图。如果我们需要输出带有回归线的散点图，菜单操作方法是通过在additionaloptions选项菜单中选择RegressionPlots/Plotsofdependentbyindependentcolumns命令，重新再提交一次。注意，此时还可以同时选择输出残差图。程序的方法是在procreg过程里增加plot语句，要注意SAS的关键字使用在plot语句中时要加小圆点，这里是预测值p关键字，增加的plot语句如下：plotweight*height=’+’p.*height=’*’/overlay;如果我们需要输出高分辩率的回归线图形，可以先在reg过程中将拟合的预测值p输出到一个SAS数据集如bclassg中，再调用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选项，直接绘制预测值的置信区间上下界。运行后，在Graph窗口得到见图322所示结果。图5图52带有回归线、95%置信线的体重与身高(weight与height)散点图此外，还可用性别来分组，分别对男生和女生进行回归分析，分别建立男生和女生的回归模型。例32.2研究耗氧量模型。这是有关身体适应性测试的例子，肺活量与一些简单的锻炼测试数据的拟合，目的是为了在锻炼测试的基础上而不是在昂贵笨重的氧气消耗测试的基础上得到方程来预测适应性。由于回归是相关的，所以理论上还应该请求共线性诊断。该数据名为fitness，这是一个对31位成年人心肺功能的调查结果，它包含的变量见下表32.8，测试的各项数据见下表32.9。表32.8fitness数据集的变量名变量名含义age年龄weight体重oxygen耗氧量runtime跑15英哩的时间（分）rstpulse休息时每分钟心跳次数runpulse跑步时每分钟心跳次数maxpulse每分钟心跳次数最大值表32.9fitness数据集中的测试数据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在这个锻炼测试数据里，我们感兴趣的是耗氧量是如何依赖于其他变量的。建立数据文件。程序如下：datafitness;inputageweightoxygenruntimerstpulserunpulsemaxpulse;cards;4489.4744.60911.37621781824075.0745.31310.0762185185………5282.7847.46710.5053170172;run;制作变量的散点图。fitness数据集中的变量较多，我们需要制作每两个不同变量oxygen、age、weight、runtime、rstpulse、runpulse和maxpulse之间的所有散点图，即散点图矩阵，共有7×6=42个散点图。我们可以通过在SAS/Insight软件中绘制散点图矩阵，操作步骤为：在SAS命令框中键入insight后按Enter，在SAS/Insight：Open对话单中，选择work.fitness数据集后单击Open按钮，将在屏幕的窗口中显示当前打开的数据集work.fitness内容，再选择菜单上的Analyze/ScatterPlot（YX）命令，在出现的ScatterPlot（YX）对话单中，把fitness数据集中的7个变量依上面的次序全部加入的Y轴和X轴的列表框中，最后单击OK。见图323所示。图323fitness变量间的散点图矩阵散点图矩阵图中第一行的六个散点图分别表示oxygen变量作为y轴,其他六个变量作为x轴的散点图,第一列的六个散点图分别表示oxygen变量作为x轴,其他六个变量作为y轴的散点图。对角线上第一个oxygen小方格里的左下角数字37.388和右上角数字60图323fitness变量间的散点图矩阵绘制散点图矩阵图是为了研究变量间的相关性，从图323中，我们发现变量runpulse与maxpulse之间存在有较强的共线性，如果在回归模型中增加方差膨胀系数（vif），共线性水平的容许值（tol），条件数（collin）选项对回归进行诊断，也会得到相同的结论。另外，我们从图中还发现耗氧量oxygen与变量runtime有较强的负相关。从下面的相关系数计算中也能得到这些相同认识。相关系数计算。编写的相关系数计算程序如下：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语句中每个自变量排到变量列表的最后，所计算的第一类平方和，也称为第二类平方和或称部分平方和。通过分析每个自变量的这两类平方和，能知道回归模型总的平方和的构成和各个自变量所贡献的平方和，进而能知道哪些自变量是最重要的回归变量，哪些回归变量可能是无关紧要的，配合参数估计的t检验，最终为缩减回归变量提供依据，到达简化模型的目的。运行后，在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个。从参数的估计值容易得到拟合的回归为：oxygen=102.238339－0．219916age+0．304735maxpulse－0.000844rstpulse－0.373164runpulse－2.680516runtime－0.072380weight方差分析类似于例30.1，多元线性回归时，FValue是除了截距外其他所有参数都为零的假设检验。RSquare同时还是多重相关系数的平方，换言之，它是自变量与预测值之间的相关系数平方。我们这里不再重复分析和计算。多元线性回归模型的第一个重要问题是模型的简化，也就是说如何正确地缩减自变量到达最优的简化模型。由于自变量除了对总方差贡献大小外，还存在着自变量间的相关性，删除哪一个自变量并不是一个简单问题。我们在这里重点分析参数估计表中的TypeISS和TypeIISS平方和，配合分析参数估计的t检验p值，对缩减自变量的依据有所掌握。类型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，实质上是指回归模型在前面模型的基础上因为增加了自变量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统计量，原假设为包含这个自变量模型与不包含这个自变量的模型是一样好。这样，例如，maxpulse变量的类型1平方和的F统计量为相应p=0.000017，拒绝原假设，接受回归模型中包含age和maxpulse比只包含age要好。变量的类型1平方和的F统计量只有weight的p值是不显著，不能拒绝回归模型中包含age、maxpulse、rstpulse、runpulse、runtime、weight与包含age、maxpulse、rstpulse、runpulse、runtime一样好。类似地我们可以计算maxpulse变量的类型2平方和的F统计量为相应p=0.03434219，拒绝原假设，接受回归模型中包含age、maxpulse、rstpulse、runpulse、runtime、weight比包含age、rstpulse、runpulse、runtime、weight要好。类型2平方和的F统计量只有rstpulse和weight不显著。两类F统计量的p值差异最大的是rstpulse自变量，这是我们首先要考虑删除的。实际上，自变量类型2平方和的F检验等同于这个自变量的参数t检验，因为F值等于t值的平方，例如，。其实，有时直观地分析两类平方和的数值大小就能得到基本的判断，runtime自变量的两类平方和都是最大的且占的比例很大，说明是回归模型中第一重要的自变量。而rstpulse自变量在第一类平方和中有比较大的数值却在第二类平方和中是最小的，这是rstpulse自变量应该被考虑第一个删除的主要原因。剔除不显著的回归变量。从参数估计表检验部分也可以看出，对自变量rstpulse和weight的回归系数的t假设检验，不能拒绝它们为0的原假设。当然在这里必须小心地看待这些检验，因为它们都是在其他自变量都加入回归的前提下进行显著性检验的，完全可能因为自变量间存在较强的相关而掩盖他们对回归的贡献。所以在剔除不显著的回归变量时必须逐个进行。另外，从自变量rstpulse和weight的回归系数接近于0，也提示我们应先考虑删除哪一个变量。因为过程reg具有连续的交互功能，在执行了提交的部分语句后，仍可继续提交语句让它执行，直至提交quit语句或因执行其他过程而终止。根据上面删除自变量的分析，我们先从已加入回归模型的自变量中剔除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