SAS的基本统计功能及应用.docx

SAS的基本统计功能及应用§4.1 一些单变量检验问题 4.1.1 正态性检验 4.1.2 两独立样本的均值检验 4.1.3 成对总体均值检验 §4.2 回归分析 4.2.1 用SAS/INSIGHT进行曲线拟合 4.2.2 用SAS/INSIGHT进行线性回归分析 4.2.3 用SAS/INSIGHT拟合广义线性模型 4.2.4 用REG过程进行回归分析 §4.3 方差分析入门 4.3.1 用ANOVA过程进行单因素方差分析 4.3.2 用NPAR1WAY进行非参数单因素方差分析 4.3.3 多重比较 4.3.4 多因素方差分析 §4.4 列联表分析 4.4.1 列联表的输入与制表 4.4.2 列联表独立性检验 4.4.3 属性变量关联度计算前面我们已经看到了SAS的编程计算、数据管理能力、数据汇总、数据探索分析能力。这一章我们讲如何用SAS进行基本的统计检验、线性回归、方差分析、列联表检验等基本统计分析。我们既使用SAS语言编程，也使用SAS/INSIGHT的菜单界面。一些单变量检验问题 对单个变量，我们可能需要作正态性检验、两独立样本均值相等的检验、成对样本均值相等的检验。正态性检验在PROC UNIVARIATE语句中加上NORMAL选项可以进行正态性检验。例如，我们要检验SASUSER.GPA 中GPA是否服从正态分布，只要用如下UNIVARIATE过程：proc univariate data=sasuser.gpa normal; var gpa;run;结果（部分）如下： Univariate ProcedureVariable=GPA College Grade Point Average Moments W:Normal 0. Pr<W 0.0001其中W:Normal为Shapiro-Wilk正态性检验统计量，Pr<W为检验的显著性概率值（p值）。当N2000时正态性检验用Shapiro-Wilk统计量，N>2000时用Kolmogorov D统计量。我们可以看到，p值很小，所以在0.05水平（或0.10水平）下应拒绝零假设，即认为GPA分布非正态。在SAS/INSIGHT中为了检验GPA的分布，先选“Analyze | Distribution”菜单打开GPA 变量的分布窗口，然后选“Curves | Test for Distribution”菜单。除了可以检验是否正态分布外还可以检验是否对数正态、指数分布、Weibull分布。两独立样本的均值检验假设我们有两组样本分别来自两个独立总体，需要检验两个总体的均值或中心位置是否一样。如果两个总体都分别服从正态分布，而且方差相等，可以使用两样本t检验过程TTEST 。比如，我们要检验SASUSER.GPA数据集中男生和女生的SATM分数是否具有相等的平均值，只要用如下程序：proc ttest data=sasuser.gpa; class sex; var satm;run;过程中用CLASS语句指定分组变量，用VAR语句指定要比较的变量。结果如下： TTEST PROCEDUREVariable: SATM Math SAT ScoreSEX N Mean Std Dev Std Error-Female 145 611. 84. 6.Male 79 565. 82. 9.Variances T DF Prob>|T|-Unequal 4.0124 162.2 0.0001Equal 3.9969 222.0 0.0001For 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分数有显著差异，女生分数要高。上面的检验中对立假设是两组的均值不等，所以检验是双边的，p值的计算公式为Pr(t 分布随机变量绝对值>计算得到的t统计量的绝对值)。如果要进行单边的检验，比如对立假设为女生分数高于男生分数（右边），则p值为Pr(t分布随机变量>计算得到的t统计量) ，当计算得到的t统计量值为正数时（现在t=4.0）此单边p值为双边p值的一半，当计算得到的t统计量为负数时肯定不能否定零假设。检验左边时恰好相反。如果我们希望检验男、女生的GPA分数则无法使用两样本t检验，因为检验女生的GPA 样本的正态性发现它非正态。这种情况下我们可以使用非参数检验。检验两独立样本的位置是否相同的非参数检验有Wilcoxon秩和检验。我们用NPAR1WAY过程加Wilcoxon选项可以进行这种检验。见下例：proc npar1way data=sasuser.gpa wilcoxon; class sex; var gpa;run;其CLASS语句和VAR与TTEST过程相同。结果如下： N P A R 1 W A Y P R O C E D U R E Wilcoxon Scores (Rank Sums) for Variable GPA Classified by Variable SEX Sum of Expected Std Dev Mean SEX N Scores Under H0 Under H0 Score Female 145 16067.5000 16312.5000 463. 110. Male 79 9132.5000 8887.5000 463. 115. Average Scores Were Used for Ties Wilcoxon 2-Sample Test (Normal Approximation) (with Continuity Correction of .5) S = 9132.50 Z = 0. Prob > |Z| = 0.5978 T-Test Approx. Significance = 0.5983 Kruskal-Wallis Test (Chi-Square Approximation) CHISQ = 0.27949 DF = 1 Prob > CHISQ = 0.5970结果分为四部分：两样本的秩和的有关统计量，Wilcoxon两样本检验的结果，t检验的近似显著性，Kruskal-wallis检验结果。我们只要看Wilcoxon检验的p值Prob > |Z| = 0.5978 ，检验结果不显著，可认为男、女生的GPA分数在0.05水平下无显著差异。SAS/INSIGHT中未提供两独立样本检验的功能。成对总体均值检验我们在现实中经常遇到两个总体是相关的测量结果的比较，比如，考察同一组人在参加一年的长跑锻炼前后的心率有无显著差异。这时，每个人一年前的心率和一年后的心率是相关的，心率本来较快的人锻炼后仍相对于其它人较快。所以，检验这样的成对总体的均值不能使用两样本t检验的方法，因为独立性条件不再满足。这时，我们可以检验两个变量间的差值的均值是否为零，这等价于检验两组测量值的平均水平有无显著差异。检验单个样本的均值是否为零只要使用UNIVARIATE过程，在UNIVARIATE过程的矩部分给出了均值为零的t检验和符号检验、符号秩检验的结果。例如，我们想知道SATM和SATV这两门考试的成绩有无显著差异（SATM平均值为595.3，SASTV平均值为504.6，我们希望知道差异是否显著）。因为这两个成绩是同一个学生的成绩，所以它们之间是相关的（学得好的学生两科一般都好，学得差的一般两科都差），不能用独立两样本的t检验，但可以计算两变量间的差DMV SATMSATV，检验差值变量的均值是否为零。如果否定，则可认为SATM和SATV的平均值有显著差异。为此，我们先用一个数据步计算差值，然后对差值变量用UNIVARIATE过程进行分析就可以得到结果。程序如下：data new; set sasuser.gpa; dmv = satm - satv; keep dmv;run;proc univariate data=new; var dmv;run;结果（部分）如下： Univariate ProcedureVariable=DMV Moments N 224 Sum Wgts 224 Mean 90.73661 Sum 20325 Std Dev 92.82931 Variance 8617.28 Skewness -0.10367 Kurtosis -0.34625 USS CSS CV 102.3063 Std Mean 6. T:Mean=0 14.62923 Pr>|T| 0.0001 Num = 0 215 Num > 0 181 M(Sign) 73.5 Pr>=|M| 0.0001 Sgn Rank 9757.5 Pr>=|S| 0.0001我们只要看其中的三个检验：T: Mean=0是假定差值变量服从正态分布时检验均值为零的t 统计量值，相应的p值Pr>|T|为0.0001在0.05水平下是显著的，所以可认为两科分数有显著差异。M(Sign)是非参数检验符号检验的统计量，其p值Pr>=|M|为0.0001在0.05水平下是显著的，结论不变。Sgn Rank是非参数检验符号秩检验的统计量，其p值Pr>=|S|为0.0001 在0.05水平下是显著的，结论不变。所以这三个检验的结论都是两科成绩有显著差异。如果t检验对立假设是单边的，其p值算法与上面讲的两样本t检验p值算法相同。在SAS/INSIGHT中比较成对样本均值的显著差异，同样是先计算两变量的差值变量（在“Edit | Variables | Other”菜单中，指定两个变量，指定两个变量间的计算为减法，则可以生成差值变量，可以用数据窗口菜单的“Define Variables”改变量名），然后对此差值变量选“Analyze | Distribution”，选“Tables | Location Tests”并选中t检验、符号检验和符号秩检验即可在分布窗口显示结果。回归分析本节先讲述如何用SAS/INSIGHT进行曲线拟合，然后进一步讲如何用SAS/INSIGHT进行线性回归，简单介绍SAS/INSIGHT的广义线性模型拟合，最后介绍如何用编程进行回归分析。用SAS/INSIGHT进行曲线拟合两个变量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。用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 EquationWEIGHT = - 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对截距项系数和各斜率项系数，给出了自由度（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 语句。用SAS/INSIGHT拟合广义线性模型经典线性回归理论的估计与假设检验要求自变量 为常数（非随机），随机误差项满足 。广义线性模型放宽了这些假设，其模型为其中因变量 （ 向量）的元素为服从指数族分布（如正态、逆高斯、伽马、泊松、二项分布）的随机变量， （ 向量）的元素为与 分布类型相同的随机误差项，元素之间相互独立，单调函数 叫做联系函数，它把因变量的均值 与自变量 （ 阵）的线性组合联系起来。 （ 向量）为回归系数。模型中每个自变量对应于设计阵 中的一列或几列， 的第一列一般元素全为1，对应于截距项。 （ 向量）是表示偏移量的变量。注：随机变量Y称为服从指数族分布，如果其分布密度（概率函数）有如下形式：其中 为自然参数或称经典参数， 为分散度参数（与尺度参数相关）， a, b, c为确定性函数。这样的自变量Y的均值和方差与参数的关系如下：为了使用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”可以把这几个点排除在外，发现结果基本不变。用REG过程进行回归分析SAS/STAT中提供了几个回归分析过程，包括REG（回归）、RSREG（二次响应面回归）、ORTHOREG（病态数据回归）、NLIN（非线性回归）、TRANSREG（变换回归）、CALIS（线性结构方程和路径分析）、GLM（一般线性模型）、GENMOD（广义线性模型），等等。我们这里只介绍REG过程，其它过程的使用请参考SAS系统SAS/STAT软件使用手册。REG过程的基本用法为：PROC REG DATA=输入数据集 选项;VAR 可参与建模的变量列表;MODEL 因变量自变量表 / 选项;PRINT 输出结果;PLOT 诊断图形;RUN;REG过程是交互式过程，在使用了RUN语句提交了若干个过程步语句后可以继续写其它的REG 过程步语句，提交运行，直到提交QUIT语句或开始其它过程步或数据步才终止。例如，我们对SASUSER.CLASS中的WEIGHT用HEIGHT和AGE建模，可以用如下的简单REG 过程调用：proc reg data=sasuser.class; var weight height age; model weight=height age;run;就可以在输出窗口产生如下结果，注意程序窗口的标题行显示“PROC REG Running”表示REG 过程还在运行，并没有终止。Model: MODEL1Dependent Variable: WEIGHT Weight in pounds Analysis of Variance Sum of Mean Source DF Squares Square F Value Prob>F Model 2 7215.63710 3607.81855 27.228 0.0001 Error 16 2120.09974 132.50623 C Total 18 9335.73684 Root MSE 11.51114 R-square 0.7729 Dep Mean 100.02632 Adj R-sq 0.7445 C.V. 11.50811 Parameter Estimates Parameter Standard T for H0: Variable DF Estimate Error Parameter=0 Prob > |T| INTERCEP 1 -141. 33. -4.230 0.0006 HEIGHT 1 3. 0. 3.973 0.0011 AGE 1 1. 3. 0.411 0.6865 Variable Variable DF Label INTERCEP 1