容许度互为倒数。
下一个结果为残差对预测值的散点图,用它可以检验残差中有无异常情况,比如非线性关系、异方差、模型辨识错误、异常值、序列相关等等。此例中各散点较随机地散布在0线的上下,没有明显的模式,可认为结果是合适的(可舍弃的不显著的变量AGE不反映在残差图中)。
图4-5-图4-8是典型的非线性、异方差、异常值、序列相
关的情况的残差图。检查序列相关的残差图以观测时间(次序)图4-4 残差对预测值散点图 为横坐标。
用Tables菜单可以加入一些其它的统计量,如做共线诊断(Colinearity Diagnostics)的条件数(Conditional Index)。用Graphs菜单可以加入残差的正态概率图(Residual Normal QQ)和偏杠杆图(Partial Leverage)。
在Vars菜单中可以指定一些变量,这些变量可以加入到数据窗口中。数据窗口的内容保存在内存中,不自动改写磁盘中的数据集,所以要保存数据窗口的修改结果的话需要用“File - Save - Data”命令指定一个用来保存的数据集名。为了了解加入的变量的具体意义,选数据窗口菜单中的“Data Options”,选中“Show Variable Labels”选项。各变量中,Hat Diag为帽子矩阵的对角线元素(帽子矩阵H恰好是n?n的),即杠杆率,反映了每个观测的影响大小。Predicted为拟合值(预报值),Linear Predictor为使用线性模型拟合的结果,在线性回归时与Predicted相同。Residual为残差。Residual Normal Quantile是残差由小到大排序后对应的标准正态的分位数,第i个残差的正态分位数用??1??i?0.375??计算,其中?为标准正态分布
?n?0.25?函数,参见1.3.7关于QQ图的解释。Standardized Residual(标准化误差)为残差除以其标准误差。Studentized Residual(学生化残差)为与标准化残差类似,但计算第i个学生化残差时预测值和方差估计都是在删除第i个观测后得到的。当学生化残差的值超过2时这个观测有可能是强影响点或异常点。
关于其它的一些诊断统计量请参考帮助菜单的“SAS System Help - Help on SAS Software Products - SAS/INSIGHT Software - Multiple Regression”,或《SAS系统:SAS/STAT软件使用手册》第一章和第九章。
在SAS/INSIGHT中,为了保存结果表格,在进行分析之前选中菜单“File - Save - Initial Tables”,这是一个状态开关,选中时输出表格画在分析窗口内的同时显示在输出(Output)窗口。如果要保存某一个表格,也可以选定此表格(单击表格外框线),然后用菜单“File - Save - Tables”。为了保存分析窗口的图形,先选定此图形,然后选“File - Save - Graphics File”,输入一个文件名,选择一种文件类型如BMP即可。为了打印某一表格或图形,先选定它,然后用菜单“File - Print”。选中“File - Save - Statements”可以开始保存SAS/INSIGHT程序
103
语句,但意义不大。
4.2.3 用SAS/INSIGHT拟合广义线性模型
经典线性回归理论的估计与假设检验要求自变量X非随机,随机误差项满足
?~N(0,?2I)。广义线性模型放宽了这些假设,其模型为
Y????
??g(?)??0?X? 其中因变量Y(n?1向量)的元素为服从指数族分布(如正态、逆高斯、伽马、泊松、二项分布)的随机变量,随机误差项?(n?1向量)的元素与Y的元素分布类型相同,元素之间相互独立,单调函数g(.)叫做联系函数(Link function),它把因变量的均值?与自变量X(n?p阵)列的线性组合联系起来。?(p?1向量)为回归系数。模型中每个自变量对应于设计阵X中的一列或几列,X的第一列一般元素全为1,对应于截距项。?0(n?1向量)是表示偏移量的变量。
注:随机变量Y称为服从指数族分布,如果其分布密度(概率函数)有如下形式:
??y?b(?)?f(y;?,?)?exp??c(y,?)??a(?)?
?? 其中?为自然参数或称经典参数,?为分散度参数(与尺度参数相关),a, b, c为确定性函
数。这样的随机变量Y的均值和方差与参数的关系如下:
E(Y)?b?(?)
Var(Y)?a(?)b??(?) 为了使用SAS/INSIGHT拟合广义线性模型,在选“Analyze - Fit (Y X)”之后,在出现的对话框中先选定因变量和自变量,然后按“Method”按钮,出现选择模型的对话框(图4-9),在这里可以选因变量的分布类型(Response Dist.),选联系函数,选估计尺度参数的方法。
图4-9 模型选择对话框
各联系函数定义如下: Identity 恒等变换 Log 自然对数 104
Logit
g(?)?log?1??0???1
Probit
g(?)???1(?),其中?为标准正态分布函数
重对数变换g(?) = log(-log(1-?)),0<1
Comp. Log-Log Power
?????0,?在对话框的Power输入框指定。
g(?)???log(?)??0 指数族中每一个分布有一个特定的联系函数,使得? = g(?) = ?,即用分布的期望值表示经典参数,这样的联系函数叫经典(canonical)联系函数。正态分布的经典联系函数为恒等变换,逆高斯分布为-2次方变换,伽玛分布为-1次方变换,泊松分布为对数变换,二项分布为逻辑变换(Logit)。注意Logit、probit、重对数变换都只适用于二项分布。
例如,SASUSER.INGOTS中存放了一个铸造厂的数据,它记录了各批铸件在一定的加热、浸泡时间条件下出现的不合要求的铸件数目。HEAT为加热时间,SOAK为浸泡时间,N
R为加热浸泡后N件铸件中还不合要求的铸件数。R应该服从二项分布,为每批铸件的件数,
其分布参数(比例)可能受加热、浸泡时间的影响。因此,我们拟合以R为因变量,以HEAT和SOAK为自变量的广义线性模型,因变量分布为二项分布,使用经典联系函数(Logit函数)。模型为:
R~B(N,?)log?1????0??1?HEAT??2?SOAK
为了拟合这样的模型,选“Analyze - Fit(Y X)”,选R为因变量,选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(广义线性模型),等等。我们这里只介绍REG过程,其它过程的使用请参考《SAS系统一一SAS/STAT软件使用手册》。 REG过程的基本用法为:
PROC REG DATA=输入数据集选项;
105
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过程还在运行,并没有终止。
这些结果与SAS/INSIGHT得到的结果是一致的。同样我们发现变量AGE的作用不显著,所以我们只要再提交如下语句:
model weight=height; run;
就可以得到第二个模型结果:
事实上,RBG提供了自动选择最优自变量子集的选项。在MODBL语句中加上“SELECTION 106
=选择方法”的选项就可以自动挑选自变量,选择方法有NONE(全用,这是缺省),FORWARD(逐步引入法),BACKWARD(逐步剔除法),STEPWISE(逐步筛选法),MAXR(最大R2增量法),MINR(最小R2增量法),RSQUARE(R2选择法),ADJRSQ(修正R2选择法),CP(Mallows的Cp统计量法)。比如,我们用如下程序:
model weight=height age / selection=stepwise; run;
可得到如下结果:
可见只有变量HEIGHT进入了模型,而其它变量(AGE)则不能进入模型。
REG过程给出的缺省结果比较少。用PRINT语句和PLOT语句可以显示额外的结果。为了显示模型的预测值(拟合值)和95%预测界限,使用语句:
print cli; run;
得到如下的结果:
107