第十讲广义线性模型

第十讲广义线性模型 • 方差分析模型 • 广义线性模型

一方差分析模型 • 概述 • 单因素方差分析模型 • 多因素方差分析模型 • 交互效应模型 • 用编程做方差分析 • 图示分析 • 拟合方差分析模型 • 生成诊断图 • 拟合交互效应模型 • 用菜单功能做方差分析 • 图示分析 • 拟合方差分析模型 • 生成诊断图 • 拟合交互效应模型

1.1 概述 方差分析是进行多个正态总体均值比较最常用的一种统计方法。按影响因素的多少，方差分析模型可分为单因素方差分析模型和多因素方差分析模型。 • 单因素方差分析模型例如某家钢铁公司采用a种不同的工艺加工钢材，现从第i种工艺生产的钢材中随机选出ni件样品，i=1,…,a。测得j样品的质量指标为有Yij，j=1,…, ni。我们希望用得到的数据来比较这些工艺的优劣。这是一个典型的方差分析问题。统计上把影响质量指标的工艺称为因素（因子,factor)，每种工艺条件为一个水平(level)。不同水平下得到的样本看成来自不同的总体。希望由此对不同水平下的各个总体的均值进行比较。对此，观测到的Yij常用以下的模型来表示：

单因素方差分析模型概述 其中表示第i种水平下产量的均值，为各水平下均值的平均值，称为各水平的效应（效应的总和为0）；为包含所有其它因素的的随机误差。如果把工艺条件因素以外的其它因素（原材料，操作人员，生产环境等等）都尽可能控制在同一条件下，且每次试验都独立进行，则不同水平下的总体都可近似看成独立且服从正态分布。因此，在方差分析中我们通常假定相互独立且有相同的分布。比较不同水平下均值是否相同的问题就归为检验如下的假设：不全相等。在上述模型中，观测到的数据{Yij}的均值受一个因素不同水平的影响，即依赖于水平。这类问题称为单因素方差分析问题

多因素方差分析模型概述 实际中影响总体均值的重要因素不止一个。例如用a种测定方法来测定b种食品的放射性同位素的含量。假定在每种食品中都随机地抽出a×n个样品，再随机地分配于不同的测定方法中，使得每种方法都测定同一种食品的n个样品。用Yijk表示第i种测定方法测定的第j种食品的第k个样品的放射性同位素含量，i=1,…, a，j=1,…, b，k=1,…, n.。试验的目的是分析不同的食品和方法是否对含量的测定值有显著影响。在上述问题中，包含有两个因素：一是测定方法，记为因素A，有a个水平；二是食品种类，记为因素B，有b个水平。类似于单因素方差分析问题的讨论，对Yijk考虑如下的模型：其中表示各因素的平均效应，表示第一个因素A第 i 个水平的附加效应，表示第二个因素B第 j 个水平的附加效应，为随机的误差,通常假定为独立具有相同正态分布,该模型就称为两因素方差分析模型.

多因素方差分析模型概述（续） 为避免模型参数，和取值上的随意性，还附加如下限制，对上述模型，要说明因素A无显著性影响，就是检验如下的假设: 不全相等。要说明因素B无显著影响，就是检验如下的假设，不全相等。而模型无显著性效果是指上述两个原假设同时成立。

交互效应模型概述 在多个因素的问题中，例如要考察添加剂和工艺条件对产品产量的影响，不同的添加剂往往会要求不同的工艺条件，在一种工艺条件下有效的添加剂，在另一种工艺条件下可能是完全无效的。在这种情况下，两因素对产量的影响不再是两者附加效果的简单叠加，而是两因素之间有着交互影响。我们称这种现象为“两个因素存在交互作用”。对于存在交互作用的观测，我们采用下面的方差分析模型，其中表示平均的效应，与前面的模型有着相同的含义，表示因素A的第i个水平和因素B的第j个水平交互作用的附加效应。这里我们也假定是独立的且服从相同的正态分布。

交互效应模型概述（续） 为避免模型参数取值的任意性，模型中还常附加如下限制，该模型称为交互效应模型。对于上述交互效应模型，我们不仅可以检验单一因素A或B的各个水平之间是否有显著差异，也可以检验因素A与B之间有无明显的交互影响，既检验如下的假设，不全为零。在上面的内容中，我们只介绍了几种常见的方差分析模型，相应的方差分析方法和过程可参考相关的统计教科书。值得注意的是，方差分析模型事实上也是线性模型，它们之间唯一不同之处在于统计分析的目的不全一样。回归分析的目的在于建立因变量和解释变量之间的回归关系，而方差分析的目的只在于检验各因素是否对因变量有显著性影响。这种检验都是通过总的离差平方和的分解来实现的。

1.2 用编程做方差分析 数据集Catalyst 的四个变量 >catalyst Temp Conc Cat Yield 1 160 20 A 60 2 180 20 A 72 3 160 40 A 54 4 180 40 A 68 5 160 20 B 52 6 180 20 B 83 7 160 40 B 45 8 180 40 B 80

1.2.1 图示分析 • plot.design( ) >plot.design(catalyst) #作各因素影响比较图 • plot.factor( ) >par(mfrow=c(1,3)) >plot.factor(catalyst) #各因素不同水平上产量分布盒子图 • interaction.plot( ) > attach(catalyst) > par(mfrow=c(2,2)) > interaction.plot(Temp,Conc,Yield) > interaction.plot(Temp,Cat,Yield) > interaction.plot(Conc,Cat,Yield) >detach(catalyst) #作各因素之间的交互作用图注：方差分析中的数据集必须是试验设计对象，否则需要用fac.design( )重建数据框

图示分析结果图（plot.design） 各因素影响比较图分析：从图中可以看出Temp是最有影响的因素，两个不同的水平180和160所对应的平均产量的差异最大，然后是Conc和Cat

图示分析结果图(plot.factor) 分析：各因素不同水平上产量分布盒子图 temperature是比较显著的因子，在Temp的两个水平上，Yield的均值有很大差异，而且在这两个水平上，分布没有重叠的部分。而对酸性浓度Conc和催化剂Cat，它们的两个不同水平对Yield的均值影响就不那么明显了。从第三组图上我们还可以看出，Yield在Cat的两个不同水平上的均值虽然很接近，方差却是明显不同

图示分析结果图(interaction.plot) 各因素之间的交互作用图分析：第二张图表明Temp与Cat之间有一定的交互影响存在

1.2.2 拟合方差分析模型 在S-Plus中，所有方差分析模型的拟合都是使用函数aov。注：anova则只是返回方差分析表 > aovcat.all<-aov(Yield~.,data=catalyst) 我们用aovcat.all储存生成的方差分析对象。模型公式 “Yield~.”表示对全部因素的方差分析。通过简单地调用该方差分析对象可得到下面的分析报告 > aovcat.all 用summary函数调用模型对象aovcat.all，可得到下面的方差分析表 > summary(aovcat.all) > fitted(aovcat.all) #返回拟合值 > predict(aovcat.all) #返回设计点的预测值

summary(aovcat.all)分析 > summary(aovcat.all) ＃输入命令＃输出结果： Df Sum of Sq Mean Sq F Value Pr(F) Temp 1 1058.0 1058.00 20.64390 0.0104688 Conc 1 50.0 50.00 0.97561 0.3792013 Cat 1 4.5 4.50 0.08780 0.7817346 Residuals 4 205.0 51.25 结果分析：从方差分析结果可以看出Temp的p值为0.01，而Conc和Cat的p值都较大，也就是说在不同的温度水平上产量Yield是有显著差异的，所以在考虑影响时，温度（Temp）应是首选的考虑因素。

1.2.3 方差分析模型诊断图 >par(mfrow=c(3,2)) >plot(aovcat.all) #对方差分析模型对象aovcat作诊断图 1、残差对拟合值散点图； 2、残差绝对值平方根对拟合值散点图； 3、因变量对拟合值散点图； 4、残差得正态QQ图； 5、残差和拟合值的取值范围图； 6、 Cook距离图。

1.2.4 拟合交互效应模型 Temp和Cat有交互影响，我们可以把这些影响考虑到模型中去。我们可以在aov函数中修改模型公式，生成新的模型对象；也可用update函数在原有的模型拟合基础上添加交互影响成分，如 >aovcat.int<-update(aovcat.all,.~.+Temp:Conc+Temp:Cat,data=catalyst) 即生成包含Temp与Conc交互作用以及Temp和Cat交互作用的新模型对象。对此模型对象直接调用和用summary函数调用，分别得到模型的拟合报告与方差分析表。 > summary(aovcat.int)＃方差分析模型的方差分析表 Df Sum of Sq Mean Sq F Value Pr(F) Temp 1 1058.0 1058.00 4232 0.00023621 Conc 1 50.0 50.00 200 0.00496281 Cat 1 4.5 4.50 18 0.05131670 Temp:Conc 1 4.5 4.50 18 0.05131670 Temp:Cat 1 200.0 200.00 800 0.00124766 Residuals 2 0.5 0.25 分析：从方差分析表中可以看出，Temp和Cat的交互影响是显著的（p值＝0.0012），而模型中其它项的p值与原模型相比也有明显的下降。

1.3 用菜单功能作方差分析 • 图示分析 Graph 2Dgraph Box Plot • 拟合方差分析模型 Statistics ANOVA  Fixed Effects • 生成诊断图方法同上,编辑Plot页面 • 拟合交互效应模型 ANOVA对话框设定Model页面 Formula模型公式： Yield~Temp+Conc+Cat+Temp:Conc+Temp:Cat 注：具体过程和结果分析同1.2

二广义线性模型 • 概述 • 背景 • 定义 • 用编程做logistic回归 • 作图分析 • 拟合模型 • 离差分析 • 单项添加 • 模型选择 • 用菜单功能做logistic回归 • 作图分析 • 模型拟合

2.1.1 背景 前面介绍的一般线性回归模型中，我们总是假定因变量Y的均值可表示为解释变量的线性函数，因变量取值为正态的连续变量。然而在实际问题中，我们经常会遇到因变量取值并非为正态分布。一类重要的情形是因变量为取离散值的属性变量，例如用来表示产品是否合格、成功与失败的二项分布变量，公司的信用等级以及文化程度等各种分类指标。对于属性变量，它的均值通常是不能表示为解释变量的线性函数。比如考虑某项新产品是否能被市场中消费者接受，因变量取值为1/0（对应于接受/拒绝）的二项分布。当然新产品是否被接受不仅取决于产品的质量、价格、品牌的知名度等等（这里我们假定共有k个因素，设它们分别为），还受到其他（随机）因素的影响。这时Y的条件数学期望是某个特定的概率，如果再简单的将概率表示为解释变量的线性组合就不合适了，因为该线性组合的取值可以是中的任何值。

2.1.2 定义 一种可行的办法是找一个变换，使得上面的条件期望经过变换后可以在中取值。一种常用的变换是称为logit变换。此时，假设要比假设要合理得多。这种变换将模型中的系统部分（即）与解释变量的线性组合“联系” 了起来，故又称之为“联系函数”（Link function）。上例中使用的是logit联系函数，相应的统计模型称为Logistic回归模型。我们同样可以取联系函数为正态分布的逆函数，相应的模型就称为probit回归模型。它们都是的广义线性模型的重要类型。

2.1.2 定义 1972年Nelder and Wedderburn统一地给出了这些模型（包括传统的线性模型）的理论框架和计算方法，并将这类模型（包括通常的线性模型）统称为广义线性模型。下面，我们来简单了解一下广义线性模型的定义。该模型假定： 1. 是n个服从指数分布族的独立样本，i＝1，…，n； 2. 是个解释变量的线性组合 3.存在一个严增可微的联系函数（Link function）g，使得与有下面的关系，其中称为模型的系统部分，称为模型的随机部分。

常见分布及其联系函数 广义线性模型可以拟合各种满足因变量服从指数分布族的数据。指数分布族是一个较为广泛的分布族，常见的重要分布如正态分布、两点分布、泊松分布、指数分布以及伽玛分布等都属于指数族。对非正态广义线性模型，经典的最小二乘已不能用于这种模型的拟合，而是采用最大似然估计方法。其拟合算法就是解广义线性模型似然方程的最大值问题，可通过Fisher记分迭代算法来实现。广义线性模型在S-Plus中的处理方式是统一的

2.2 用编程做logistic回归 数据集为S—Plus的内部数据集kyphosis，含有四个变量，81个观测，来自81个做过脊椎矫正手术的儿童。数据分析的目的在于建立变量kyphosis和其它变量的回归关系，用于解释术后驼背的风险的大小。由于因变量是取值为1/0（分别代表present/absent）的属性变量，因此直接用线性回归模型来拟合数据是不恰当的。为此，我们考虑使用广义线性模型中的logistic回归。

广义线性模型函数glm( ) >fitted.model<-glm(formula,family=family(link=link),data=data.frame) 其中data.frame为各变量所在的数据框，formula为模型公式，family为分布族，link为联系函数（不同的分布选择不同的联系函数），fitted.model用来储存广义线性模型的拟合结果，称为广义线性模型的模型对象（其分类属性为glm）。其它选项设置细节可以通过Help系统了解。

2.2.1 作图分析 在这个例子中，我们感兴趣的是Age、Number和Start对手术后是否出现驼背的影响，可以通过盒子图来做因子分析，了解一下大致情况。 > par(mfrow=c(1,3)) > plot.factor(kyphosis) 图形反映出的信息表明，三个变量分别对于Kyphosis的两个不同取值表现出了一定的分布差异，说明三个变量对Kyphosis都有一定的影响。而变量Start表现出来的分布差异最大，因而作用也最明显。因子分析图

2.2.2 拟合模型 > kyp.fit1<-glm(Kyphosis~.,family=binomial,data=kyphosis) #做包含所有变量的logistic回归，当分布族family设定为二项分布binomial时，link函数的默认值即为logit，无须设定，拟合得到的结果存储到kyp.fit1中，为一个属性为glm的模型对象，可以直接输入此名称，得到简短输出。调用summary( )可得拟合的详细信息。 > summary(kyp.fit1) #结果见下页分析：第一部分是用来说明本次运行所调用的函数、数据来源及相应的参数设定。第二部分是残差序列的简单统计描述。第三部分是估计的系数、标准差和t统计量值。可以看到，t值很高的Start是最显著的，说明对于是否出现驼背，Start是一个有着重要影响的因素，而Age和Number的影响不是很明显。由于GLM模型估计采用极大似然估计，回归参数是通过迭代得到的，所以没有给出p值。Null Deviance给出了零模型拟合下的离差（相当于通常线性模型中的残差平方和，用以度量模型的拟合优度）以及所对应的自由度；Residual Deviance描述的是拟合模型下的离差以及所对应的自由度。对比中我们可以看到模型中的三个变量在降低离差方面所起的作用。最后一行给出了模型拟合时，Fisher记分算法的迭代次数。

输出结果summary( ) *** Generalized Linear Model *** Call: glm(formula = Kyphosis ~ Age + Number + Start, family = binomial(link = logit), data = kyphosis, na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50,trace = F)) Deviance Residuals: Min 1Q Median 3Q Max -2.312363 -0.5484308 -0.3631876 -0.1658653 2.16133 Coefficients: Value Std. Error t value (Intercept) -2.03693225 1.44918287 -1.405573 Age 0.01093048 0.00644419 1.696175 Number 0.41060098 0.22478659 1.826626 Start -0.20651000 0.06768504 -3.051043 (Dispersion Parameter for Binomial family taken to be 1 ) Null Deviance: 83.23447 on 80 degrees of freedom Residual Deviance: 61.37993 on 77 degrees of freedom Number of Fisher Scoring Iterations: 5

2.2.3 离差分析(Analysis of Deviance) 通过离差分析，我们可以了解各解释变量在模型中的相对影响和作用。 > anova(kyp.fit1,test="Chi") Analysis of Deviance Table Binomial model Response: Kyphosis Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(Chi) NULL 80 83.23447 Age 1 1.30198 79 81.93249 0.2538510 Number 1 10.30593 78 71.62656 0.0013260 Start 1 10.24663 77 61.37993 0.0013693 这是一张离差分析表，给出在零模型基础上依次（模型公式中变量的次序）加入变量Age、Number、Start对拟合效果的影响。如在零模型上增加变量Age，离差减少1.30198，对应的自由度为1；剩余离差为81.93249，对应的自由度为79。通过卡方检验的p值可以看出Number,Start都是显著的，但是如果Start变量在前拟合，则会得到Number变量的影响并不显著。如 Kyphosis~Start+Number+Age

2.2.4 单项添加 addl( ) > kyp.fit1.null <- glm(Kyphosis ~ 1, family = binomial, data = kyphosis) #先产生一个零模型拟合 > add1(kyp.fit1.null, ~ . + Age + Number + Start) #再使用add1函数对零模型来做单项添加 Single term additions Model: Kyphosis ~ 1 Df Sum of Sq RSS Cp <none> 81.00000 83.02500 Age 1 1.29546 79.70454 83.75454 Number 1 10.55222 70.44778 74.49778 Start 1 16.10805 64.89195 68.94195 Cp 值用来比较没有嵌套的模型，较小的Cp值对应的模型较好。从上面的分析可见，Start明显是首先要考虑加入到模型中的一个重要变量，而添加 Age反而增加了Cp 值。

2.2.5 模型选择 在广义线性模型中也可以使用S-Plus提供的step函数来做模型选择我们用前面生成的全模型对象kyp.fit1作为搜索的初始模型，键入下面的命令 > step(kyp.fit1) #输出结果略P167-168 这里由于全模型kyp.fit1已被设定为搜索的初始模型，所以默认的搜索过程为单项删除过程。从上面的输出信息可知，模型中包含全部解释变量是最佳选择。也可以考虑是在包含交互项的模型中搜索。 > step(glm.object,list(upper = ~.^2, lower = ~ Age)) 上面的模型搜索范围被设定为：上限为包含所有变量的直到2阶的交互作用模型，下限为只有一个变量Age的模型。搜索的结果表明还是包含全部变量的全模型为最优。

2.3 用菜单功能做logistic回归 • 作图分析 2Dplot作图形分析：Graph 2Dgraph box plot • 模型拟合 Statistics Regression  Logistic 设定各页面即可

第十讲 广义线性模型