第4章 SAS的基本统计分析功能
前面我们已经看到了SAS的编程计算、数据管理能力、数据汇总、数据探索分析能力。这一章我们讲如何用SAS进行基本的统计检验、线性回归、方差分析等基本统计分析。我们既使用SAS语言编程，也使用SAS/INSIGHT的菜单界面。
§4.1 一些单变量检验问题
对单个变量，我们可能需要作正态性检验、两独立样本均值相等的检验、成对样本均值相等的检验。
4.1.1 正态性检验
在PROC UNIVARIATE语句中加上NORMAL选项可以进行正态性检验。例如，我们要检验SASUSER.GPA中GPA是否服从正态分布，只要用如下UNIVARIATE过程：
proc univariate data=sasuser.gpa normal;
var gpa;
run;
结果（部分）如下：
Univariate Procedure
Variable=GPA College Grade Point Average
Moments
…………
W:Normal 0.951556 Pr<W 0.0001
…………
其中W:Normal为Shapiro-Wilk正态性检验统计量，Pr<W为检验的显著性概率值（p值）。当N≤2000时正态性检验用Shapiro-Wilk统计量，N>2000时用Kolmogorov D统计量。我们可以看到，p值很小，所以在0.05水平（或0.10水平）下应拒绝零假设，即认为GPA分布非正态。
在SAS/INSIGHT中为了检验GPA的分布，先选"Analyze | Distribution"菜单打开GPA变量的分布窗口，然后选"Curves | Test for Distribution"菜单。除了可以检验是否正态分布外还可以检验是否对数正态、指数分布、Weibull分布。
4.1.2 两独立样本的均值检验
假设我们有两组样本分别来自两个独立总体，需要检验两个总体的均值或中心位置是否一样。如果两个总体都分别服从正态分布，而且方差相等，可以使用两样本t检验过程TTEST。
比如，我们要检验SASUSER.GPA数据集中男生和女生的SATM分数是否具有相等的平均值，只要用如下程序：
proc ttest data=sasuser.gpa;
class sex;
var satm;
run;
过程中用CLASS语句指定分组变量，用VAR语句指定要比较的变量。结果如下：
TTEST PROCEDURE
Variable: SATM Math SAT Score
SEX N Mean Std Dev Std Error
-----------------------------------------------------------------------------
Female 145 611.77241379 84.02056171 6.97752786
Male 79 565.02531646 82.92937599 9.33028376
Variances T DF Prob>|T|
---------------------------------------
Unequal 4.0124 162.2 0.0001
Equal 3.9969 222.0 0.0001
For H0: Variances are equal, F' = 1.03 DF = (144,78) Prob>F' = 0.9114
结果有三个部分：两个总体的SATM简单统计量，两样本均值的检验，以及两样本方差是否相等的检验。标准的两样本t检验要求两总体方差相等，所以第三部分结果检验两样本方差是否相等。如果检验的结果为相等，则可使用精确的两样本t检验，看第二部分结果的Equal那一行。如果方差
检验的结果为不等，则只能使用近似的两样本t检验，看第二部分结果的Unequal那一行。这里我们看到方差检验的p值为0.9114不显著，所以可以认为方差相等，所以我们看Equal行，p值为0.0001在0.05水平下是显著的，所以应认为男、女生的SATM分数有显著差异，女生分数要高。
上面的检验中对立假设是羇riables"改变量名），然后对此差值变量选"Analyze | Distribution"，选"Tables | Location Tests"并选中t检验、符号检验和符号秩检验即可在分布窗口显示结果。
§4.2 回归分析
本节先讲述如何用SAS/INSIGHT进行曲线拟合，然后进一步讲如何用SAS/INSIGHT进行线性回归，简单介绍SAS/INSIGHT的广义线性模型拟合，最后介绍如何用编程进行回归分析。
4.2.1 用SAS/INSIGHT进行曲线拟合
图 1 身高对体重的散点图及回归直线
两个变量Y和X之间的相关关系经常可以用一个函数来表示，一元函数可以等同于一条曲线，实际工作中经常对两个变量拟合一条曲线来近似它们的相关关系。最基本的"曲线"是直线，还可以用多项式、样条函数、核估计和局部多项式估计。其模型可表示为
例如，我们要研究SASUSER.CLASS数据集中学生体重与身高之间的相关关系。为此，我们可以先画出两者的散点图（Analyze | Scatter plot）。从图中可以看出，身高越高的人一般体重越重。我们可以把体重作为因变量、身高作为自变量拟合一条回归直线，只要选"Analyze |
Fit (Y X)"，并选体重为Y变量，身高为X变量，即可自动拟合出一条回归直线，见图 1。窗口中还给出了拟合的模型方程、参数估计、诊断信息等，我们在下一小节再详细介绍。
在拟合了直线后，为拟合多项式曲线，只要选"Curves |
Polynomial"，然后输入阶次（Degree(Polynomial)），就可以在散点图基础上再加入一条多项式曲线。对于本例，我们看到二次多项式得到的曲线与直线差别很小，所以用二次多项式拟合没有优势。还可以试用三次、四次等多项式。为了改变阶次还可以使用拟合窗口中的多项式阶次滑块
（Parametric Regression Fit中的Degree(Polynomial)）。
样条曲线是一种非参数回归的曲线拟合方法。光滑样条为分段的三次多项式，曲线在每一段内是一个三次多项式，在两段的连接点是连续、光滑的。为拟合样条曲线，只要选"Curves |
Spline"，使用缺省的GCV准则（广义交叉核实）来选取光滑系数（光滑系数c越大，得到的曲线越光滑，但拟合同时变差，光滑系数c小的时候得到的曲线较曲折，而拟合较好），就可以在散点图的基础上画出样条曲线。可以用光滑系数c的滑块来调整曲线的光滑程度/拟合优度。对于本例，
GCV准则得到的样条曲线与回归直线几乎是重合的，说明直线拟合可以得到满意的结果。
核估计是另一种非参数回归的曲线拟合方法。它定义了一个核函数 ，例如使用标准正态分布密度曲线，然后用如下公式估计经验公式 ：
其中 为光滑系数， 越大得到的曲线越光滑。为了画核估计曲线，只要选"Curves | Kernel"，权重函数使用缺省的正态核，选取光滑系数的方法采用缺省的GCV法，就可以把核估计图附加到散点图上。本例得到的核估计曲线与回归直线、样条曲线有一定差别。可以手动调整光滑系数
的值，可以看到，当 过大时曲线不仅变光滑而且越来越变水平，因为这时的拟合值基本是一个常数，这与样条曲线的情形不同，样条曲线当 增大时曲线变光滑但不趋向与常数（水平线）。
局部多项式估计（Loess）是另一种非参数回归的曲线拟合方法。它在每一自变量值处拟合一个局部多项式，可以是零阶、一阶、二阶，零阶时与核估计相同。SAS/INSIGHT缺省使用一阶（线性）局部多项式。改变Loess的系数alpha可以改变曲线的光滑度。alpha增大时曲线变光滑，而且使
用一阶或二阶多项式时曲线不会同时变水平。
固定带宽的局部多项式是另一种局部多项式拟合方法。它有一个光滑系数c。
4.2.2 用SAS/INSIGHT进行线性回归分析
上面我们已经看到，用菜单"Analyze | Fit (Y X)"就可以拟合一条回归直线，这是对回归方程
的估计结果。这样的线性回归可以推广到一个因变量、多个自变量的情况。线性模型写成矩阵形式为
其中 为 向量， 为 矩阵，一般第一列元素全是1，代表截距项。 为 未知参数向量， 为 随机误差向量，元素独立且方差为相等的 （未知）。正常情况下，系数的估计为 ，拟合值（或称预报值）为 ，其中 是 空间内向 的列张成的线性空间
投影的投影算子矩阵，叫做"帽子"矩阵。拟合残差为 ，残差平方和为 ，误差项方差的估计为（要求设计阵 满秩）均方误差（MSE） ，在线性模型的假设下，若设计阵 满秩， 和 分别是 和 的无偏估计，系数估计的方差阵
。判断回归结果优劣的一个重要指标为复相关系数平方（决定系数） （其中 ），它代表在因变量的变差中用模型能够解释的部分的比例，所以 越大说明模型越好。
例如，我们在"Fit (Y X)"的选择变量窗口选Y变量（因变量）为体重（WEIGHT），选X变量（自变量）为身高（HEIGHT）和年龄（AGE），则可以得到体重对身高、年龄的线性回归结果。下面对基本结果进行说明。
回归基本模型：
WEIGHT = HEIGHT AGE
Response Distribution: Normal
Link Function: Identity
回归模型方程：
Model Equation
WEIGHT = - 141.2238 + 3.5970 HEIGHT + 1.2784 AGE
拟合概况：
Summary of Fit
Mean of Response 100.0263 R-Square 0.7729
Root MSE 11.5111 Adj R-Sq 0.7445
其中Mean of Response为因变量（Response）的均值，Root MSE叫做根均方误差，是均方误差的平方根，R-Square即复相关系数平方，Adj R-Sq为修正的复相关系数平方，其公式为 ，其中 当有截距项时取1，否则取0，这个公式考虑到了自变量个数 的多少对拟合的影响，原来的
随着自变量个数的增加总会增大，而修正的 则因为 对它有一个单调减的影响所以 增大时修正的 不一定增大，便于不同自变量个数的模型的比较。
方差分析表：
Analysis of Variance
Source DF Sum of Squares Mean Square F Stat Prob > F
Model 2 7215.6371 3607.8186 27.2275 0.0001
Error 16 2120.0997 132.5062 . .
C Total 18 9335.7368 . . .
这是关于模型是否成立的最重要的检验。它检验的是 ：模型中所有斜率项系数都等于零，这等价于说自变量的线性组合对因变量没有解释作用。它依据的是一个标准的方差分解，把因变量的总离差平方和（C
Total）分解为能用模型解释的部分（Model）与不能被模型解释的部分（随机误差，Error）之和，如果能解释的部分占的比例大就否定 。F统计量（F Stat）就是这个比例（用自由度修正过）。从上面结果看我们这个模型很显著（p值不超过万分之一），所以可以否定 。
第三类检验：
Type III Tests
Source DF Sum of Squares Mean Square F Stat Prob > F
HEIGHT 1 2091.1460 2091.1460 15.7815 0.0011
AGE 1 22.3880 22.3880 0.1690 0.6865
这个表格给出了对各斜率项是否为零（ ）的检验结果。检验利用的是所谓第三类平方和（Type III
SS），又叫偏平方和，它代表在只缺少了本变量的模型中加入本变量导致的模型平方和的增加量。比如，HEIGHT的第三类平方和即现在的模型平方和减去删除变量HEIGHT的模型的模型平方和得到的差。第三类平方和与模型中自变量的次序无关，一般也不构成模型平方和的平方和分解。表中
用F统计量对假设进行了检验，分子是第三类平方和的均方，分母为误差的均方。实际上，当分子自由度为1时，F统计量即通常的t检验统计量的平方。从表中可见，身高的作用是显著的，而年龄的作用则不显著，有可能去掉年龄后的模型更好一些。
参数估计及相关统计量：
Parameter Estimates
Variable DF Estimate Std Error T Stat Prob >|T|
INTERCEPT 1 -141.2238 33.3831 -4.2304 0.0006
HEIGHT 1 3.5970 0.9055 3.9726 0.0011
AGE 1 1.2784 3.1101 0.4110 0.6865
Parameter Estimates
Tolerance Var Inflation
. 0.0000
0.3416 2.9276
0.3416 2.9276
图 2 残差对预测值散点图
对截距项系数和各斜率项系数，给出了自由度（DF），估计值（Estimate），估计的标准误差（Std Error），检验系数为零的t统计量，t统计量的p值，检验共线性的容许度（Tolerance）和方差膨胀因子（Var Inflation）。其中自变量 的容许度定义为1减去
对其它自变量的复相关系数平方，因此容许度越小（接近0），说明 对其它自变量的复相关系数平方大，即 可以很好地被其它自变量的线性组合近似，这样 在模型中的作用不大。记 ，则 ， 叫做方差膨胀因子，它代表
的系数估计的方差的比例系数，显然其值越大说明估计越不准确，也说明 在模型中的作用不大。方差膨胀因子与容许度互为倒数。
下一个结果为残差对预测值的散点图，用它可以检验残差中有无异常情况，比如非线性关系、异方差、模型辨识错误、异常值、序列相关等等。此例中各散点较随机地散布在0线的上下，没有明显的模式，可认为结果是合适的（多余的不显著的变量AGE不反映在残差图中）。
用Tables菜单可以加入一些其它的统计量。用Graphs菜单可以加入残差的正态概率图（Residual Normal QQ）和偏杠杆图（Partial Leverage）。
在Vars菜单中可以指定一些变量，这些变量可以加入到数据窗口中。数据窗口的内容保存在内存中，不自动改写磁盘中的数据集，所以要保存数据窗口的修改结果的话需要用"File | Save |
Data"命令指定一个用来保存的数据集名。为了了解加入的变量的具体意义，选数据窗口菜单中的"Data Options"，选中"Show Variable Labels"选项。各变量中，Hat Diag为帽子矩阵的对角线元素（帽子矩阵 恰好是
的）,即杠杆率，反映了每个观测的影响大小。Predicted为拟合值（预报值），Linear Predictor为使用线性模型拟合的结果，在线性回归时与Predicted相同。Residual为残差。Residual Normal Quantile是残差由小到大排序后对应的标准正态的分位数，第 个残差的正态分位数用
计算，其中 为标准正态分布函数。Standardized Residual（标准化误差）为残差除以其标准误差。Studentized Residual（学生化残差）为与标准化残差类似，但计算第 个学生化残差时预测值和方差估计都是在删除第
个观测后得到的。当学生化残差的值超过2时这个观测有可能是强影响点或异常点。
关于其它的一些诊断统计量请参考帮助菜单的"Extended Help | SAS System Help: Main menu | Help for SAS Products | SAS/INSIGHT | Techniques | Multiple Regression"，或《SAS系统：SAS/STAT软件使用手册》第一章和第九章。
在SAS/INSIGHT中，为了保存结果表格，在进行分析之前选中菜单"File | Save | Initial Tables"，这是一个状态开关，选中时输出表格画在分析窗口内的同时显示在输出（Output）窗口。如果要保存某一个表格，也可以选定此表格（单击表格外框线），然后用菜单"File | Save |
Tables"。为了保存分析窗口的图形，先选定此图形，然后选"File | Save | Graphics File"，输入一个文件名，选择一种文件类型如BMP即可。为了打印某一表格或图形，先选定它，然后用菜单"File | Print"。选中"File | Save | Statments"可以开始保存SAS/INSIGHT语句。
4.2.3 用SAS/INSIGHT拟合广义线性模型
经典线性回归理论的估计与假设检验要求自变量 为常数（非随机），随机误差项满足 。广义线性模型放宽了这些假设，其模型为
其中因变量 （ 向量）的元素为服从指数族分布（如正态、逆高斯、伽马、泊松、二项分布）的随机变量， （ 向量）的元素为与 分布类型相同的随机误差项，元素之间相互独立，单调函数 叫做联系函数，它把因变量的均值 与自变量 （ 阵）的线性组合联系起来。 （
向量）为回归系数。模型中每个自变量对应于设计阵 中的一列或几列， 的第一列一般元素全为1，对应于截距项。 （ 向量）是表示偏移量的变量。
注：随机变量Y称为服从指数族分布，如果其分布密度（概率函数）有如下形式：
其中 为自然参数或称经典参数， 为分散度参数（与尺度参数相关），a, b, c为确定的函数函数。这样的自变量Y的均值和方差与参数的关系如下：
图 3 模型选择对话框
为了使用SAS/INSIGHT拟合广义线性模型，在选"Analyze | Fit (Y X)"之后，选定因变量和自变量，然后按"Method"按钮，出现选择模型的对话框，在这里可以选因变量的分布类型（Response Dist.），选联系函数，选估计尺度参数的方法。
各联系函数定义如下：
Identity 恒等变换
Log 自然对数
Logit
Probit ，其中 为标准正态分布函数
Comp. Log-Log
Power ， 在对话框的Power输入框指定。
对指数族中每一个因变量分布有一个特定的联系函数，使得
，这样的联系函数叫经典（canonical）联系函数。正态分布的经典联系函数为恒等变换，逆高斯分布为-2次方变换，伽玛分布为-1次方变换，泊松分布为对数变换，二项分布为逻辑变换（Logit）。注意Logit、probit、复合重对数变换都只适用于二项分布。
例如，SASUSER.INGOTS中存放了一个铸造厂的数据，它记录了各批铸件在一定的加热、浸泡时间条件下出现的不能开始轧制的铸件数目。HEAT为加热时间，SOAK为浸泡时间，N为每批铸件的件数，R为加热浸泡后N件铸件中还不能开始轧制的铸件数。R应该服从二项分布，其分布参数（比例）
可能受加热、浸泡时间的影响。因此，我们拟合以R为因变量，以HEAT和SOAK为自变量的广义线性模型，因变量分布为二项分布，使用经典联系函数（Logit函数）。模型为
为了拟合这样的模型，选"Analyze | Fit(Y
X)"，选R为Y变量，选HEAT和SOAK为自变量，按"Method"钮，选因变量分布为二项分布（Binomial），选变量N然后按"Binomial"钮，两次OK后即可以得到模型拟合窗口。可以看到，这个模型是显著的，但变量SOAK没有显著影响。去掉变量SOAK重新拟合模型。可以看出，HEAT的系数为0.0807
是正数，说明加热时间越长不能轧制的件数越多。考察拟合结果窗口下方的残差对预报值图可以发现在右下方有三个异常点，用刷亮方法选定它们，可以看到，这三个观测都是总共只有一个铸件的，所以对一般结果意义不大。选"Edit | Observations | Exclude in
Calculation"可以把这几个点排除在外，发现结果基本不变。
4.2.4 用REG过程进行回归分析
SAS/STAT中提供了几个回归分析过程，包括REG（回归）、RSREG（二次响应面回归）、ORTHOREG（病态数据回归）、NLIN（非线性回归）、TRANSREG（变换回归）、CALIS（线性结构方程和路径分析）、GLM（一般线性模型）、GENMOD（广义鲜匝Ｎ颐窍Ｍ勒馕逯峙谱咏汉习宓哪ニ鹆坑
形尴灾畋穑绻尴灾畋鹞颐窃谘」菏本筒槐乜悸悄囊桓龈湍ザ恍杩悸羌鄹竦纫蛩兀绻峁邢灾钜煸蛴悸鞘褂媚湍バ院玫呐谱印Ｕ饫铮蛩厥墙汉习宓呐谱樱副晡ニ鹆浚备髦峙谱咏汉习迥ニ鹆坑邢灾钜焓保得饕蛩氐娜≈刀灾副暧邢灾挠跋臁Ｋ裕讲罘治
龅慕崧凼且蛩囟灾副暧形尴灾跋臁Ｗ⒁猓涞姆讲罘治鲋慌卸弦蛩氐母魉接形尴灾钜欤还芰礁鲆蛩刂涫欠裼胁钜欤热缢滴颐堑奈甯雠谱蛹词褂兴母雠谱用挥邢灾钜欤挥幸桓雠谱拥慕汉习灞日馑母雠谱拥亩己茫崧垡彩撬狄蛩厥窍灾模蛞蛩氐母魉郊溆邢灾钜臁
?nbsp;
方差分析把指标的方差分解为由因素的不同取值能够解释的部分，和剩余的不能解释的部分，然后比较两部分，当能因素解释的部分明显大于剩余的部分时认为因素是显著的。方差分析假定观测是彼此独立的，观测为正态分布的样本，由因素各水平分成的各组的方差相等。在这些假定满足
时，就可以用ANOVA过程来进行方差分析。其一般写法为
PROC ANOVA DATA＝数据集;
CLASS 因素;
MODEL 指标＝因素;
RUN;
比如，为了分析SASUSER.VENEER中各种牌子的胶合板的耐磨性有无显著差别，首先我们假定假设检验使用的检验水平为0.05，可以使用如下程序进行方差分析：
proc anova data=sasuser.veneer;
class brand;
model wear=brand;
run;
结果如下：
Analysis of Variance Procedure
Class Level Information
Class Levels Values
BRAND 5 ACME AJAX CHAMP TUFFY XTRA
Number of observations in data set = 20
Dependent Variable: WEAR Amount of material worn away
Source DF Sum of Squares Mean Square F Value Pr > F
Model 4 0.61700000 0.15425000 7.40 0.0017
Error 15 0.31250000 0.02083333
Corrected Total 19 0.92950000
R-Square C.V. Root MSE WEAR Mean
0.663798 6.155120 0.14433757 2.34500000
Source DF Anova SS Mean Square F Value Pr > F
BRAND 4 0.61700000 0.15425000 7.40 0.0017
结果可以分为四个部分，第一部分是因素水平的信息，我们看到因素只有一个BRAND，它有5个水平，分别是ACME、AJAX、CHAMP、TUFFY、XTRA。共有20个观测。第二部分就是经典的方差分析表，表前面指明了因变量（指标）为WEAR，第一列"来源"说明方差的来源，是模型的（可以用方差分
析模型解释的），误差的（不能用模型解释的），还是总和。第三列为平方和，其大小代表了各方差来源作用的大小。第二列为自由度。第四列为均方，即平方和除以自由度。第五列F值是F统计量的值，其计算公式为模型均方除以误差均方，用来检验模型的显著性，如果不显著说明模型对
指标的变化没有解释能力。第六列是F统计量的p值。由于这里p值小于0.05（我们的检验水平），所以模型是显著的，因素对指标有显著影响。结果的第三部分是一些与模型有关的简单统计量，第一个是复相关系0.005就可以保证总错误率不超过0.05，但是，这样得到的实际总第一类错误率
可能要比预定的水平小得多。在MEANS语句中使用BON语句可以执行Bonferroni t检验，缺省总错误率控制水平为0.05。对上面的胶合板数据增加如下语句：
means brand / bon;
run;
结果如下：
Bonferroni (Dunn) T tests for variable: WEAR
NOTE: This test controls the type I experimentwise error rate, but
generally has a higher type II error rate than REGWQ.
Alpha= 0.05 df= 15 MSE= 0.020833
Critical Value of T= 3.29
Minimum Significant Difference= 0.3354
Means with the same letter are not significantly different.
Bon Grouping Mean N BRAND
A 2.6000 4 TUFFY
A
B A 2.3750 4 XTRA
B A
B A 2.3750 4 CHAMP
B A
B A 2.3250 4 ACME
B
B 2.0500 4