1.            相关系数计算

如果用菜单操作方法，可选择Globals/SAS/Assist/Data Analysis/Elementary/Correlation命令，再选择Active data set为study.bclass，Columns to be correlated为weight和height，然后提交运行。直接编写调用相关系数计算的程序为：

 

运行后，在Output窗口得到如表32.6所示的结果。

 

表32.6 身高与体重(weight与height)的相关系数

 

Correlation Analysis 2 'VAR' Variables:  WEIGHT   HEIGHT Simple Statistics Variable           N        Mean     Std Dev         Sum     Minimum     Maximum WEIGHT            40    47.66250    10.07415        1907    29.10000    78.10000 HEIGHT            40   153.25000    10.47525        6130   125.00000   172.00000   Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0 / N = 40                                       WEIGHT            HEIGHT WEIGHT          1.00000           0.70844                0.0               0.0001 HEIGHT          0.70844           1.00000             0.0001            0.0

    从输出表32.6可以看出，身高与体重之间的相关系数为0.70844。

 

2.            回归分析

如果用菜单操作方法，可选择Globals/SAS/Assist/Data Analysis/Regression/Linear regression命令，再选择Active data set为study.bclass，Dependent为weight，Independent为height，然后提交运行。编程实现回归方法为：

 

其中，模型参数r表示要输出残差分析，包括因变量的观察值、由输入数据和估计模型来计算的预测值、残差值、标准误差、学生化残差、COOKD统计量。模型参数clm表示对每个观察输出因变量期望值的95%置信上界和下界，仅考虑到参数估计的偏差，没有考虑误差项的偏差。模型参数cli表示对因变量的各个预测值输出95%置信上界和下界，这个置信界反映了误差的偏差以及参数估计的偏差。模型参数dw表示要进行误差项的独立性检验，计算Durbin-Watson统计量。运行后，在Output窗口得到如表32.7所示的结果。

 

表32.7    回归分析结果

 

Model: MODEL1   Dependent Variable: WEIGHT                                                                        Analysis of Variance(方差分析)                                           Sum of        Mean                 Source           DF      Squares       Square       F Value      Prob>F                  Model            1   1986.48457   1986.48457       38.287       0.0001                  Error           38   1971.56918     51.88340                  C Total         39   3958.05375                       Root MSE       7.20301     R-square       0.5019                     Dep Mean      47.66250     Adj R-sq       0.4888                     C.V.          15.11254                                       Parameter Estimates(参数估计)                                  Parameter      Standard     T for H0:                             Variable   DF      Estimate         Error    Parameter=0    Prob > |T|                INTERCEP   1     -56.748575   16.91239600        -3.355        0.0018                HEIGHT     1       0.681312    0.11010770         6.188        0.0001 误差项的独立性检验 Durbin-Watson D             1.471 (For Number of Obs.)           40 1st Order Autocorrelation   0.185 置信区间      Dep Var   Predict   Std Err  Lower95%  Upper95%  Lower95%  Upper95%             Std Err  Obs  WEIGHT     Value   Predict      Mean      Mean   Predict   Predict  Residual  Residual  1   43.1000   42.0417     1.457   39.0925   44.9908   27.1647   56.9187    1.0583     7.054  2   55.8000   44.7669     1.231   42.2743   47.2595   29.9737   59.5602   11.0331     7.097  3   33.6000   35.2286     2.310   30.5527   39.9044   19.9155   50.5417   -1.6286     6.823 …… 38   52.7000   46.8109     1.147   44.4885   49.1332   32.0453   61.5764    5.8891     7.111 39   60.8000   57.0305     1.895   53.1953   60.8658   41.9529   72.1082    3.7695     6.949 40   78.1000   60.4371     2.358   55.6639   65.2103   45.0940   75.7802   17.6629     6.806 残差分析         Student                     Cook's   Obs  Residual    -2-1-0 1 2            D     1     0.150  |      |      |     0.000     2     1.555  |      |***   |     0.036     3    -0.239  |      |      |     0.003     4     1.728  |      |***   |     0.067     5    -0.104  |      |      |     0.001     6    -0.749  |     *|      |     0.010 7     1.879  |      |***   |     0.053 ……    35    -0.110  |      |      |     0.000    36     1.242  |      |**    |     0.027    37     0.154  |      |      |     0.001    38     0.828  |      |*     |     0.009    39     0.542  |      |*     |     0.011    40     2.595  |      |***** |     0.404   Sum of Residuals                      0 Sum of Squared Residuals      1971.5692 Predicted Resid SS (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，即H₀： ，F(1,38)=MSR/MSE=1986.48457/51.88340=38.287，查F分布表，p值为0.0001小于显著水平0.05，表明可拒绝原假设，并有足够的证据断定回归线的斜率不为零。所以，这一模型拟合数据比基线模型好。

无偏的误差估计标准值Root MSE= = 7.20301，因变量weight平均值Dep Mean=47.66250，变异系数（或称方差系数）CV=Root MSE/Dep Mean×100=7.20301/47.66250×100=15.11254，它表示与单位无关的方差。R‑Square是0～1之间的值，它表示贡献给模型而不是贡献给拟合残差的总方差的那部分，它也称为决定系数或拟合优度，用于判断回归模型拟合好坏。R²=1－ESS/TSS=RSS/TSS=1986.48457/3958.05375= 0.5019，调整 R²=1－ESS/TSS×(N-i)/（N-k-i）=1-1971.56918/3958.05375×39/38=0.4888，R²越是接近1说明模型拟合得越好，等于1则说明完全拟合，没有任何信息丢失，本例的R²值表明有一半信息丢失没有被回归模型表示出来，通常R²应该超过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%，有充足的理由否决截距为零的原假设。估计斜率的标准误差计算公式见式(31.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因变量（Dep Var）的40条原始观察值和回归模型的预测均值（Predict Value），及预测均值的标准差（Std Err Predict）、预测均值的置信区间下界（Lower95% Mean）和上界（Upper95% Mean）、预测值的置信区间下界（Lower95% Predict）和上界（Upper95% Predict）、残差（Residual）、残差的标准差（Std Err Residual）。我们以第一条观察（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，学生化残差（Student Residual）=残差/标准残差=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个残差值，即 。

3.            输出带有回归线的散点图

如果我们需要输出带有回归线的散点图，菜单操作方法是通过在additional options选项菜单中选择Regression Plots/Plots of dependent by independent columns命令，重新再提交一次。注意，此时还可以同时选择输出残差图。程序的方法是在proc reg过程里增加plot语句，要注意SAS的关键字使用在plot语句中时要加小圆点，这里是预测值p关键字，增加的plot语句如下：

 

如果我们需要输出高分辨率的回归线图形，可以先在reg过程中将拟合的预测值p输出到一个SAS数据集如bclassg中，再调用gplot过程绘制图形。增加的output语句如下：

 

绘制高分辨率的带有回归线的散点图程序如下：

 

注意，我们也可以用图形自带i=rlcli95选项，直接绘制预测值的置信区间上下界。运行后，在Graph窗口得到如图32.2所示的结果。

 

图32.2  带有回归线、95%置信线的体重与身高(weight与height)散点图

从绘制出的带有回归线的图形可形象地看出模型拟合数据比只用均值预测的基线模型好。仔细观察图形，两条95%的上下预测值置信带呈现两头喇叭口。此外，还可用性别来分组，分别对男生和女生进行回归分析，分别建立男生和女生的回归模型。

 

例32.2  研究耗氧量模型。这是有关身体适应性测试的例子，肺活量与一些简单的锻炼测试数据的拟合，目的是为了在锻炼测试的基础上而不是在昂贵笨重的氧气消耗测试的基础上得到方程来预测适应性。由于回归是相关的，因此，理论上还应该请求共线性诊断。该数据名为fitness，这是一个对31位成年人心肺功能的调查结果，它包含的变量见表32.8，测试的各项数据见表32.9。

 

表32.8  fitness数据集的变量名

 

表32.9  fitness数据集中的测试数据

 

在这个锻炼测试数据里，我们感兴趣的是耗氧量是如何依赖于其他变量的。