广义线性模型()

广义线性模型广义线性模型*(Nelder和Wedderburn,1972)除了正态分布，也允许反应分布，以及模型结构中的一定程度的非线性。GLM具有基本结构g(Pi)=XP，其中山三E(Yi),g是光滑单调'链接函数'，Xi是模型矩阵的第i行，X和卩是未知参数的向量。此外，GLM通常会做出Yi是独立的和Yi服从一些指数族分布的假设。指数族分布包括许多对实际建模有用的分布，如泊松分布，二项分布，伽马分布和正态分布。GLM的综合参考文献是McCullagh和Nelder(1989)，而Dobson(2001)提供了一个全面的介绍。因为广义线性模型是以“线性预测器”X®的形式详细说明的，所以线性模型的许多一般想法和概念通过一些修改而继续存在到广义线性模型中。除了必须选择的链接函数和分布之外，基本模型公式与线性模型公式基本相同。当然，如果恒等函数被选择作为链接以及正态分布，那么普通线性模型将作为特例被恢复。然而，泛化是以某种成本为代价的：现在的模型拟合必须要迭代完成，而且用于推理的分布结果是近似的，并且由大样本限制结果证明是正确的而不是精确的。但在深入探讨这些问题之前，请考虑几个简单的例子。Pi=cexp(bti),例1：在疾病流行的早期阶段，新病例的发生率通常会随着时间以指数方式增加。因此，如果m是第ti天的新病例的预期数量，则该形式的模型为请注意，“广义”和“一般”线性模型之间存在区别-后一个术语有时用于指除简单直线以外的所有线性模型。可能是合适的，其中c和b是未知参数。通过使用对数链路，这样的模型可以变成GLM形式log(^i)=log(c)+bti=P0+tiP1(根据®0=logc和®1=b的定义)。请注意，模型的右侧现在在参数中是线性的。反应变量是每天新病例的数量，因为这是一个计数，所以泊松分布可能是一个合理的可以尝试的分布。因此，针对这种情况的GLM使用泊松反应分布，对数链路和线性预测器®0+ti®1。axi例2：狩猎动物捕获猎物的速度yi往往随着猎物密度xi的增加而增加，但最终会趋于平衡，当捕食者捕获尽可能多的猎物时。对于这种情况一个合适的模型可能是其中a是未知参数，表示最大捕获率，h是未知参数，表示捕获率为最大速率一半时的猎物密度。很显然，这个模型在其参数中是非线性的，但是通过使用倒数链路，右边的参数可以是线性的：（这里卩0三1/a和卩1三h/a）。在这种情况下，猎物捕获率的标准差可能与平均速率大致成比例，建议使用Gamma分布作为反应，并完成模型设定。们不限于示例的简单直线形式，但可以有对于线性模型可能具有的线性预测器的任何结构。2.1GLMs的理论未完成。fe（y）=exp［｛ye?b（e）｝/a©）+c（yw）］,GLM中的反应变量可以来自指数族的任何分布。如果一个分布的概率密度函数或概率质量函数可以写成，则该分布属于指数族分布其中b,a和c是任意函数，（p是任意的“尺度”参数，并且0被称为分布的“典范参数”（在GLM上下文中，0将完全依赖于模型参数卩，但是目前没有必要做这个明确）。例如，很容易看出，正态分布是指数族的一员，因为加"）=鳩卿］-◎諾T=旳「2+芻_"2_1。咖陌）一=exp"一豊一bg（（7厉）］c（0,y）=一沪/（2。）一log（V027r）三-y2/（2a2）-log（oV27r）.表2.1给出了在R中为GLM实施的指数族成员的类似分解。用a，b和p可以得到指数族分布的均值和方差的一般表达式。给定一个特定的y0的对数似然性仅仅只是将log［f0（y）］视为0的一个函数。那是E（鲁）=［E（Y）_y（0）］/a（0）将1作为一个随机变量来处理，通过用随机变量Y替换特定的观测值y，可以评估?1/?0的期望值：使用E（?1/?0）=0这个一般结果，（在0取真值时，参见2.4节中的（2.14））和重新排列意味着E（Y）=b0（0）。（2.1）即任何指数族随机变量的均值由的一阶导数给出。0，其中b的形式取决于特定的分布。该等式是将GLM的模型参数卩与指数族的典范参数联系起来的关键。在GLM中，参数卩决定了反应变量的均值，并且通过（2.1）,它们决定了每个反应观测值的典范参数。需=_b〃（&）/a（0）再次对似然性微分处理b〃(0)/a(0)=E[(Y—b«))2]/a(0)：并将其插入到一般结果中，E(?21/?G2)=-E[(?1/?G)2](衍生物在真实0值下计算，参见结果(2.16)，第2.4节)，重新安排第二个有用的一般结果：var(Y)=b00(O)a(¥).a原则上可以是申的任何函数，并且当与GLM—起工作时，如果申是已知的，处理任何形式的a都是没有困难的。然而，当申未知，事情就会变得很尴尬，除非我们可以写出(申)=9/«，其中®是一个已知常数。事实上，这种限制形式涵盖了所有有实际意义的案例(见表2.1)。a(申)=9/®允许基于正态分布的模型中的不等方差的可能性，但是在大多数情况下，®仅为1•因此，我们现在有var(Y)=b00(O)申/①.(2.2)在随后的章节中，我们通常会将var(Y)视为□三E(Y)的函数，这会很方便，并且由于卩和0通过(2.1)链接，我们总能定义一个函数V(卩)=b00(0)/①，使得var(Y)=V(卩)申。表2.1列出了几个这样的功能。2.1.2拟合广义线性模型回想一下，GLM模拟独立反应变量的n阶向量Y,其中卩三E(Y)，通过gg)=xp和Yi〜f0i(yi),其中f0i(yi)表示指数族分布，典范参数0i由pi(通过方程2.1)决定，因此最终由卩决定。给定Y的一个观测值向量y，卩的最大似然估计是可能的。由于Yi是相互独立的，卩的似然函数是nL(P)=Yfoi(yi),i=1因此卩的对数似然函数是n"0)=^log由血)]»=1n=工皿％一址(弘)]/逐(0)+6(©yi)右边对卩的依赖是通过0i对卩的依赖。请注意，函数a，b和c可能随着i而变化-例如，允许不同的二项分母ni,对于每个二项式反应的观测值，或对于正常反应的不同方差(但在常数内已知)。另一方面，对于所有的i,假设申是相同的。正如前一节所讨论的那样，对于实际工作来说，只考虑可以写出ai（9）=^/®i的情况就足够了，其中曲是已知常数（通常为1），在这种情况下通过偏导令结果表达式为零并求出卩。通过这个链式法则所以微分（2.1）,我们可以得到然后推出dl丄dl丄寸［s—耳（仇）］0他徭妙他）/0福将（2.1）和（2.2）代入最后一个方程，意味着求解卩的方程是(2.3)(2.4)然而，如果权重v（M）事先已知且与卩独立，那么这些方程就正好是为了通过非线性加权最小二乘寻找卩而必须解出的方程。在这种情况下，最小二乘的目标是其中m非线性地依赖于卩，但权重v（山）被视为固定的。要找到最小平方估计值，包括求解?S/?pj=Oj但当V（M）项被视为固定时，这个方程组很容易被看作是Q.3）。这种对应立即提出了一种求解（2.3）的迭代方法。令卩［k］表示在第k次迭代处的估计参数向量，并且令n［k］和Mk］为具有元素ni［k］=Xi卩［k］和y［ik］=g-1（币［町）的向量，其中g-1（•）是链路的反函数。从参数估计开始，卩［0］，迭代以下步骤，直到P［k］的序列收敛：计算当前P［k］隐含的V（p［ik］）项。为了得到卩［k+1］（V（p［ik］）被视为固定的而不是作为卩的函数）设置k到k+1实际上，这种方法比需要的要慢。第2步本身涉及迭代，但在「丿」已经收敛之前实际上将非线性最小二乘方法迭代到收敛没有多大意义。因此，第2步通常被替换为：2•使用P［k］作为初始来获得卩［k+1］。应用这种方法会产生一个相当紧凑和简洁的方案。为了看到这个，让我们以矩阵形式写出非线性最小二乘问题。定义对角矩阵V［k］,其中V［k］ii=V（卩阴）（2.4）附近的一阶泰勒展开替代，所以因此，没有进一步的近似Jij=?^i/?Pj|P?[k].Now由“伪数据”的定义对角线权重矩阵畔=帀说呼因此，下面的步骤可迭代到收敛使用当前的Mk］和n［k］计算伪数据z［k］和迭代权重W［k］。将对于卩的平方和最小化，以便获得卩［k+1］,因此n［k+l］=X卩［k+1］和p［k+l］。将k增加1。收敛的卩解决了（2.3）,因此是即的最大似然估计。该算法在大多数实际情况下趋于一致，但也有例外（例如，二项数据的不良或过于灵活的模型）。请注意，要开始迭代，我们只需要卩［0］和n［0］的值，不需要卩［0］。因此，迭代通常通过设置卩［0］i=yi和ni［0］=g（p［0］i）开始，并根据需要轻微调整卩［0］i,以避免无限的币［0］（例如，如果yi=0且有对数链路）。该方法被称为IterativelyRe-weightedLeastSquares（IRLS），原因很明显，在此背景下，归因于Nelder和Wedderburn（1972）。2.1.3IRLS的目标是对数似然的二次近似IRLS迭代中的工作线性模型不仅仅是寻找参数的最大似然估计的手段。在一个加性常数内（在收敛时）也是模型在卩附近的对数似然的二次近似。显然，第一个衍生工具是关于对数似然和S之间的Pj匹配：实际上它们都是零。S的二阶导数矩阵为-xwx®,并且这被证明与对数似然的期望的二阶导数矩阵匹配，并且因此在大样本极限中由大数定律来匹配二阶导数矩阵本身。为了证明这一点，首先将u定义为关于模型参数的对数似然的导数向量，因此ui=?l/?pi，然后将（2.3）中的导数以矩阵向量形式重写为u=XTG?1V?1（y?^）/^.?注意，如果V（m）被视为P的函数，算法不会最小化（2.4）,因为在这种情况下，令导数为零不会产生（2.3）。换句话说，最大似然与具有平均方差关系的最小二乘基本上不同。E（uuT）=XTG?1V?1E［（Y?卩）（Y?卩）t］V?1G?1X/3=XTG?1V?1VV?1G?1X/9=XtWX/9导数的这种对应关系足以证明S是卩附近的对数似然的二次近似，并且由于MLE的一致性，它们在真实参数值附近。2.1.4AICforGLMs通过对可能性进行直接比较的模型选择存在如下问题：如果将冗余的参数添加到正确的模型中，可能性几乎总是增加（并且从不减小），因为额外的参数让模型更接近数据，即使这只意味着对数据的组成部分进行“噪声建模”。正如在线性模型的情况下，如果我们能够根据它们拟合数据均值卩而不是数据y的能力来选择模型，这个问题会得到缓解。在GLM背景下，合理的方法是根据模型最大化1（卩卩）而不是1（卩y）的能力来选择模型，但是为了做到这一点，我们必须能够估计出1（卩卩）。实际上这个估计很明确并且因为y=y时，这也必须成立“—制滴（耳—XB）『然后，这个论据的（1.15）（只修改权重）产生了估计量'1（卩？y）?tr（A）+n/2A=X（XTWX）?iXTWandhencetr（A）=p,（可识别的）模型参数的数量。因此，在模型之间进行选择时，我们会选择哪个模型具有最高的1（卩）-p值，这相当于选择了Akaike信息准则的最低值（Akaike,1973）,AIC=2[-1（卩）+p]。前面的论点假设申已知。如果不是，那么就需要一个申的估计来计算AIC，结果AIC中的惩罚项p将变为p+1。这种归纳的大样本分布GLM的分布结果并不精确，反而是基于大样本近似，利用包括一致性在内的最大似然估计的一般性质（见2.4节）。根据最大似然估计量的一般性质，我们认为，在大样本限制下，即〜N（卩,1?】），其中I=E（uuT）是模型参数的信息矩阵，u是关于模型参数的对数似然函数的导数向量，因此即〜N(卩,(XtWX)?®).对于具有已知尺度参数申的分布，可以直接使用此结果来查找参数的置信区间，但是如果尺度参数未知（例如对于正态分布），则必须估计它，并且区间必须基于一个合适的t分布。尺度考虑检验Ho：g@）=Xo%相反Hi：g（p）=Xi卩1，其中卩是反应向量Y的期望，Y的元素是来自指数族分布的相同成员的独立随机变量，其中X0?X1。如果我们有反应向量的观测值y，则可以执行广义似然比检验。令1（卩0）和1（卩1）为两个模型的最大似然率。如果H0为真，则在大样本限制中，2[1（卩?丿?1（卩?°）〕〜X2p1?P0，（2.5）sothattwicethe?whichshou1dstrict1ybeamaximum1ike1ihoodestimate原假设是假的，则模型1倾向于具有比模型0高得多的可能性，从而两倍于严格应该是最大似然估计的?，或者在大样本限制下倾向于MLE的估计值。对数似然对于相关的X2分布的一致性差异太大。如果可以计算相关模型的对数似然性，那么近似结果（2.5）仅仅是有用的。在由IRLS估算的GLM情况下，只有当尺度参数申已知时才是这种情况。因此，结果可以与泊松和二项模型一起直接使用，但不能与正态§，伽马或反高斯分布（其中尺度参数未知）一起使用。稍后将简短地讨论在后面这些情况下做什么。偏差在实践中使用GLM时,在普通的线性建模D=2[1（'?max）?l（p?）]e（2.中，用与残差平方和类似的方式来=n6）解释数量是有用的。2Xi=yi（e??e?）?b（e?）+b（（e?）,iiiiii（2.7）这个数量是模型的偏差，被定义为thetermsinsidethesummation其中1（pmax）表示饱和模型的最大可能性：每个数据点具有一个参数的模型。在给定数据的情况下，1（卩max）是可能性具有的最大值，并且可以通过简单地设置y=y并计算出可能性来计算。G?和0分别表示关于饱和模型和兴趣模型的典范参数的最大似然估计。注意如何将偏差定义为与申独立。表2.1列出了单个数据对偏差的贡献，对于若干分布-这些是偏差定义中求和中的项（术语）。scaleddeviance,与偏差相关的是调整偏差,D?=D®,

这取决于尺度参数。对于二项分布和泊松分布，其中9=1，偏差和比例偏差是相同的，但通常情况并非如此。通过广义似然比检验结果（2.5），我们可以预计，如果模型是正确的，那么近似D?〜X2n?P，（2.8）tojustify（2.8）asalargesampleapproximationundermanycircumstances在大样本限制中。实际上，这样的论点是捏造的，因为证明（2.5）的限制论证依赖于模型中参数的数量保持不变，而样本容量趋于无穷大，但饱和模型具有与数据一样多的参数。渐近结果可用于表2.1中的一些分布，以在许多情况下将（2.8）证明为大样本近似值（参见当然，对于正态分布和一致性链接，我们使用第1章的结果。这对于正态分布情形下是精确的。但是，请注意，它完全打破了二进制数据的二项式。鉴于偏差的定义，很容易看出，本节开始的似然比检验可以通过重新表达两次对数似然比统计量作为。然后在H0下（2.9）（在大样本极限内），其中Di*是具有pi个可识别参数的模型i的偏差。但是，这只有在尺度参数已知时才有用，以便可以计算D*。与未知9的模型比较在H0下我们有近似的结果D；~E“andD；~xJ,并且，f/L门：汕」/人如果被看作是渐近独立的，这就意味着（久-（久-巧）/血-內）

卩/仇一pi）在大样本限制下（当然，在普通线性模型特例中结果确实如此）。F的有用性质是它可以在不知道9的情况下进行计算，这可以从比率收益的顶部和底部取消，在H0下，近似结果（Do—DJ/（pi—po）..（2.10）这个结果的优点是，当9未知时，它可用于基于模型比较的假设检验。缺点是对于Di的可疑的分布假设以及它所基于的独立性近似。当然，一个明显的替代方法是使用估计值9来获得每个模型的估计值Di*=Di9，然后使用（2.9）进行假设检验。然而，如果我们使用估计（2.11）来达到这个目的，则很容易看出它只是（n-pl）xF,所以我们的测试将完全等同于使用F比率结果（2.10），但是使用Fp1-pO,^作为参考分布。显然直接使用（2.10）是一种更保守的方法，因此通常是首选：它至少在估计尺度参数时考虑到了不确定性。正如我们已经看到的那样，参数卩的MLE可以在不知道尺度参数申的情况下获得，但是在这些参数未知的情况下，通常必须进行估计。近似结果（2.8）提供了一个明显的估计量。一个x2n-p随机变量的期望是n-p，所以等同于观测值达到我们的近似预期值申?D=D/?（n?p）.（2.11）第二个估计量基于Pearson统计量，定义为显然，X2®是一组零均值，单位方差，随机变量的平方和，具有n-p个自由度，这表明如果模型适当，则近似X2®?x2n-p:这个近似值是有根据的。将观测到的Pearson统计量设置为我们得到的预期值申?=X?2/（n?p）.请注意，它很明确的表明这一点\、、\\「xj,其中W和z是IRLS权重和伪数据，在收敛时进行评估。2.1.8典范链接函数典范链接gc用于分布，是链接函数，使得gc（m）=ei,其中Oi是分布的典范参数。例如，对于泊松分布，典范链接是对数函数（其他示例见表2.1）。典范链接的使用意味着Oi=Xip（其中Xi是X的第i行）。典范链接倾向于具有一些很好的性质，例如确保卩保持在反应变量的范围内，但它们也具有更多微妙精细的优点，其中之一在此处被推导出。回想一下，似然最大化涉及对每个Pj的对数似然性微分，并令结果为零，以获得方程组但是，如果正在使用典范链接，贝y?Oi/?Pj=Xij,如果像通常情况，wi=1?i，这个方程组可简化为XTy?XTp?=O,i.e.toXTy=XT^?.即，简化为现在考虑X包含1列的情况：任意其他加权求和，其中权重由模型矩阵P给出，这意味着该系统中的一个方程简单地为iyi=ipi。?Recallthatif｛Zii2〜X2n・回想一下，如果｛Zi：i=1...n｝是一组独立同分布，N（0,1）r.v然后是PZi2-%2n。列（或这些的线性组合）在原始数据和拟合值之间保存。其中一个实际结果是，对于任何具有截距项和典型联系的GLM,残差总和为零：这种“观察到的无偏性”是一个令人放心的特性。结果的另一个实际用途是在使用对数线性模型的分类数据分析中，其中它通过规定模型提供了一种方法，可以保存任何模型中保存的内置于研究设计中的总数。2.1.9残差模型检查可能是应用统计建模中最重要的部分。在普通线性模型情况下，这是基于对模型残差的检验，模型残差包含数据中的所有信息，而不是由模型的系统部分来解释。残差检验也是GLM案例模型检验的主要手段，但在这种情况下，残差标准化是必要的，并且更加困难。对于GLM而言，不仅仅只检查原始残差的主要原因是难以检查原始残差检验假设的均值方差关系的有效性。例如，如果采用泊松模型，则残差的方差应与拟合值（山）的大小成正比增加。然而，如果将原始残差与拟合值作图，则需要一种杰出的能力来判断残差的可变性是否与均值成比例增加，而不是均值的平方根或平方。出于这个原因，通常将GLM残差标准化，以这种方式，如果模型假设正确，则标准化残差应该具有大致相等的方差，并且尽可能与来自普通线性模型的残差相似（尽管见第6.5节中图6.9的替代绘图方法）。皮尔森残差根据拟合模型，将残差标准化的最明显的方法是将它们除以与其标准差成比例的量。这导致了皮尔森残差如果模型正确，它应该具有近似为零的均值和方差申。当与拟合值或任何协变量（无论是否包含在模型中）进行绘图时，这些残差不应显示均值或方差的任何趋势。“皮尔森残差”这个名字与相关的皮尔森残差的平方和给出了请注意，皮尔森残差是来自收敛IRLS方法的工作线性模型的残差除以收敛IRLS权重的平方根。异常残差在实践中，皮尔森残差的分布在零附近可以是非常不对称的，因此它们的行为不像可能希望的那样接近普通线性模型残差。异常残差在这方面通常是优选的。异常残差是通过注意到偏差对GLM起到与残差平方和对普通线性模型起相同作用来得到：确实，对于普通线性模型来说，偏差是残差平方和。在普通线性模型情况下，偏差由平方残差和组成。这就是具有适当标示的偏差组成部分的平方根的残差。因此，把di写成第i个数据所贡献的偏差的一部分（即（2.7）中求和的第i项），我们有并且通过与普通线性模型类比，我们可以定义程=sign（s：-切皿.

根据需要，这些“异常残差”的平方和给出了偏差本身。现在如果计算了所有参数已知的模型的偏差，则(2.8)将变成D*?x2n,这可能表明对于单个数据di?xl2，可推出'⑴」。当然，(2.8)不能合理地应用于单个数据，但对于良好的拟合模型来说，这表明我们可能预期异常残差具有类似于N(0,1)随机变量的行为，特别是在对(2.8)预计是一个合理的近似值的情况下。拟极大似然迄今为止，GLM的处理假定反应变量的分布是指数族的已知成员。如果有充分的理由假设该反应变量遵循特定的分布，则它对该分布上的基本模型很有吸引力，但是在许多情况下，反应分布的性质并不十分清楚，并且可能它只能详细说明反应变量的方差与其均值之间的关系。也就是说，函数V(卩)可以被指定，但是好不了多少。然后出现的问题是，是否有可能开发用于拟合和推理GLM的理论，从仅指定平均方差关系的位置开始。事实证明，基于拟极大似然的概念，开发出令人满意的方法是可能的。考虑一个具有均值pi和已知方差函数V(pi)的随机变量的观测值yi。然后对于给定yi的pi的对数拟极大似然定义为QiM=X班QiM=X班_名*TOyi(2.12)正如我们将会看到的，这个函数的关键特征是它具有li的许多有用特性，即对应于单个观测的对数似然，但只需要V的知识而不是Yi的全部分布。假设数据是独立随机变量的观测值，我们可以定义所有反应数据的平均向量p的对数拟极大似然或定义p的任意参数向量为ng(“)=工©(血)q的关键特征是，为了推导GLM,它表现得与对数似然函数非常相似，但只需要方差函数就可定义它。例如，考虑获得GLM参数卩的最大拟似然参数估计。对q关于Pj产量微分&q_s—冷。他鯛-召0叫)羽，所以参数估计是方程的解但这正是方程组(2.3),必须解决这个问题才能找到GLM的。因此，最大拟似然参数估计可以通过通常的GLMIRLS方法找到，在任何情况下只需要知道V(p)。此外，对数拟似然性与参数估计值对数似然性具有足够的性质，并且对于参数的最大拟似然估计量值也适用。类似地，当对数似然度l被对数拟似然函数q代替时，请注意，饱和模型的对数拟似然始终为零，因此GLM的拟偏差很简单Dq=?2q(p?)卩很显然，关于残差和尺度参数估计的讨论也从似然性中延续到拟似然的情况，但不超过用q代替l。拟似然性方法的实际应用要求对（2.12）中的积分进行评估，但对于大多数实际有用的V的形式，这是可能的：McCullagh和Nelder（1989）给出了例子，或者在R中你可以输入例如以获得所执行的任何特定均值方差关系的qi形式。对于表2.1中对应于指数族分布的均值方差关系，拟偏差的形式精确对应于该族的偏差形式。拟似然性的一个主要实际用途是提供一种来模拟比泊松或二项分布（具有其固定尺度参数）更可变的计数数据的方法：拟似然性方法假定申未知。这种“过度分散”的数据在实践中很常见。另一个实际用途是提供一种对具有平均方差关系的数据建模的方法，其中没有明显的指数族分布：例如，预期方差与均值成比例的连续数据。2.2GLM的几何图形GLM和GLM拟合的几何图形比普通线性模型的几何图形更难，因为用于判断模型拟合的可能性通常不意味着拟合可以通过模型和数据之间的欧几里得距离来判断。图2.1用一个含有2个参数的GLM的3个数据拟合Gamma分布和对数链接的例子说明了GLMs的几何情况。1.4节的平坦模板子空间现在被一个弯曲的“模型流形”所取代，该模型流形由模型可预测的所有可能的拟合值向量组成。由于模型流形与数据之间的欧几里德距离不再被用来进行拟合度的测量，因此必须采用不同的方法来说明估计的几何结构。图2.1右端板的黑线显示了所有反应变量的组合，这些变量产生了相同的估计模型。注意这些线条通常并不平行，并且通常不与模型流形正交。要充分理解图2.1，它可能有助于考虑一些不同的含有2个参数模型的图形。1.对于一个普通的线性模型，模型流形将是一个平面，所有的相等线拟合将是正交的（因此相互平行）。2•对于假定正态分布（但是非一致性链接）的GLM,相等线拟合将与它们遇到的模型流形的（切线空间）正交。对于拟合4个数据的含有2个参数的模型，相等线拟合将变成相等平面拟合。一般来说，图2.1所示的几何图形适用于任何GLM。有了更多的数据，相等线拟合就变成了n-p维相等面的拟合，其中n和p分别是数据和参数的数量：对于任何固定的卩，方程（2.3）给出了定义这样一个平面的y的限制。请注意这些平面x图2.1GLM的几何图形。左端板展示了所示的三个x,y数据的广义线性模型E（y）三y=exp（卩0+卩収）的最佳拟合，假设每个yi是由模型给出的带有均值的Gamma分布随机变量的一个观测值。右端板展示了使用此模型作为示例的GLM拟合的几何图形。所示的单位立方体表示向量（yl,y2,y3）T定义单点空间。灰色表面根据模型显示所有可能的预测值（在单位立方体内），即它表示所有可能的（卩1,卩2,卩3）T值。由于参数卩0和卩1允许在所有可能的值上变化，这就是相应模型“拟合值”追踪的表面：“模型流形”。连续线从立方体的一个面开始,在另一个面上离开，是等效拟合线：位于这样一条线上的反应数据（yl，y2，y3）T的值各自导致相同的卩0，卩1的最大似然估计，即相同的（卩1，卩2，卩3）T。注意等效拟合线既不相互平行，也不与模型流形正交。可以相交-稍后将返回的点。对于离散反应数据，图片没有什么不同，尽管只有在可能性连续概括下（通常可以通过在概率函数中用合适的伽玛函数代替阶乘）才能严格得到相等线拟合。只有正态分布才能与模型流形有正交的相等线/平面拟合。对于其他分布，相等线/平面拟合有时可能彼此平行，但决不会全部与模型流形正交。2.2.lIRLE的几何图形通过考虑一个参数模型对2个反应数据的拟合，IRLS估计算法的几何图形是最容易理解的。图2.2说明了这种模型的几何结构：在这种情况下，一个具有对数链接和Gamma错误的GLM，图2.2GLM的几何图形。E（yi）三yi=20exp（-pxi）其中yi?Gamma和i=1,2。左端板展示了拟合2个x,y数据的模型（连续线）的最大似然估计值，显示为?。右端板显示了拟合几何图形。15x15的正方形是空间＜2中的一部分，其中（yl,y2）定义为单点。粗体曲线是“模型流形”：它根据模型包含所有可能的点（卩1,卩2）（即，当卩变化时，（卩1,卩2）可以追踪该曲线）。细线是相等线拟合的例子。位于这些线之一上的所有点（yl,y2）共享相同的卩的卩1,卩2）：这个MLE是相等拟合线切割模型流形的地方。对于卩=.1,2,3,4,5,6,7,8,9,1,1.2,1.5,2,3,4绘制相等线拟合。（卩=.1,.7和2用实线表示，当卩=2时候的线接近图形的底部。当卩=.1时候的线在该图的绘图区外，但出现在随后的图中。）但是可以通过链接和分布假设的任意组合来为GLM构建类似的图片。现在拟合GLM的关键问题是模型流形不平坦，并且相等线拟合与它们遇到的模型流形不正交。IRLS方法对拟合问题进行线性转换和重新调整,因此在当前的卩的估计值时，模型流形和相等线拟合是正交的，并且在重新调整的空间中，当前卩的估计值的位置由X乘以当前的卩估计值给出。这种重新调整会产生一个拟合问题，可以将其视为局部线性，从而可以用最小二乘法更新卩估计值。图2.3说明了IRLS步骤如何涉及形成伪数据并对其进行加权，从而将拟合问题有效地转化为可以通过线性最小二乘法近似解决的问题。该图说明了IRLS步骤中涉及的转换，这些转换重复进行，直到IRLS方法被迭代到收敛。（a）（b）11（c）（d）11图2.3基于图2.2所示的GLM的IRLS估计的几何图形。（a）显示拟合问题的几何图形-模型流形是粗黑色曲线，等拟合线是细线（如图2.2），数据是?和拟合值的当前估计值卩［k］。（b）问题以当前拟合值为新中心（yi被yi-y［i］代替）。（c）问题被重新线性调整，使得X的列现在跨越模型流形的切线空间。切线空间用灰线表示（这个步骤用g0（y［ik］）（yi-4ik］））代替yi-|i［ik］。（d）该问题是线性转换的，因此现在由邓［k］给出位置。对于大多数GLM，现在必须通过将相对于每个轴的分量乘以〈Wi来再次调整问题，其中Wi是迭代权重：这将确保通过的相等估计线与切线空间正交。在当前的例子中，这些权重都是1，因此所需的正交性已经成立。现在对于转换后的问题，在模型流形的附近，可以用切线空间来近似，其中相等线拟合近似正交：因此，可以通过找到转换数据的最小二乘投影来获得卩和卩的更新估计，到切线空间（灰线）。图2.4拟合和收敛问题的几何图形。具有对数链接和正常错误的含有1个参数的GLM的几何图形被显示。厚曲线是模型流形-根据模型，在单位平方内包含所有可能的数据拟合值。细线是相等线拟合（水平如图2.2）。注意在图的左上角，相等线拟合是如何相互交叉的。该重叠区域中的数据将产生具有多于一个参数值的局部最小值的模型可能性。考虑IRLS拟合方法的运行，表明在这种情况下，根据用于启动拟合过程的初始值，它可以收敛到不同的估计值。?说明问题反应向量的位置，用于说明文本中的非独特收敛。2.2.2IRLS收敛的几何图形图2.4说明拟合模型的几何图形，E（yi）三卩i=exp（-®xi），其中yi是正态分布的，有两个数据yi拟合，其中xl=0.6和x2=1.5。如前两节所示，相等线拟合显示在一个图形上，其中反应向量（yl，y2）T将定义一个单点，并根据该模型，所有可能的拟合值（卩1，卩2）T的集合显示为一条厚曲线。在这个例子中，相等线拟合在图的左上角相交并交叉（对应于非常差的模型拟合）。这个交叉点是有问题的：特别是，IRLS拟合位于左上角的数据的结果将取决于IRLS处理开始的初始参数估计，因为每个这样的数据点位于两条相等拟合线的交点处。如果IRLS迭代从图形右上角的拟合值开始，则接近右上角的拟合值将被估计，然而用图形左下角的拟合值开始迭代将导致估计的拟合值不同，并更接近图形的左下角。在实践中这确实发生，可以很容易的在R中展示出来,通过拟合数据yl=.O2,y2=.9，如图2.4所示，可以很容易地证明。请注意，这里的第二个拟合实际上具有较高的可能性（较低的偏差）-拟合在可能性方面不是等同的。引起这些模糊的拟合几何图形的类型并不总是会发生：例如，一些模型具有相等拟合的平行线/平面，但对于任何具有相等线拟合的交叉线/平面的模型而言，存在一些模糊性的范围。幸运的是，如果模型是一个好的模型，那么位于模糊区域的数据往往是不太可能的。在图2.4的例子中，问题区域完全由模型只能很差地拟合的数据组成。由此可见，非常差的数据模型可能会产生这种类型的估计问题：但对于非常差的模型而言，将模型作为任何复杂数据集的早期尝试的特征并不罕见。如果遇到这样的问题，那么通过对转换后的反应数据进行线性建模可能会更好，直到已经确定足够好的候选模型转换回GLM为止。当然，如果选择合理的起始值，那么在拟合GLM时拟合过程中的模糊性不太可能引起重大问题：毕竟，该算法将收敛到可能性的局部最小值之一。然而，当它在可选最小值之间循环而不会收敛时，模糊性可以通过“性能迭代”在GAM估计中引起更严重的收敛问题，。2.3GLMswithRglm函数提供了在R中使用GLM的手段。它的使用与lm函数的使用类似，但有两点区别。模型公式的右侧指定线性预测变量的形式，现在给出了反应均值的链接函数，而不是直接使用反应均值。此外，glm需要一个族论点，用于指定要使用的指数族的分布以及与之一起使用的链接函数。在本节中，将介绍glm函数与各种简单GLM的用法，以说明GLM涵盖的各种模型结构。分布模型和心脏病如果对于心脏病患者能够提供最好的照顾，那么能够尽早发现心脏病是再好不过的。关于检测心脏病的一种方法是检测血液中肌酐激酶（CK）的水平。CKvaluePatientswithHeartattackPatientswithoutheartattack202886013261003081403051802102201912601813001313401913801504207046080表2.2作为CK水平函数的心脏病发作概率的数据（来自Hand等人，1994）。血流。一项研究被进行（Smith，1967），其中测量了怀疑患有心脏病的360名患者的CK水平。后来经过更长时间的医学调查后，确定每名患者是否确实有心脏病发作。数据见表2.2。原始报告根据CK水平的范围对患者进行分类，但在表中仅给出了该范围的中点。能够基于这些数据的诊断标准是很好的，这样CK水平可以用来估计患者发生心脏病的概率。我们可以通过构建一个模型，试图从CK水平解释心脏病发作的患者比例，从而朝着这样一个目标迈进。在接下来的数据被读入一个名为heart的data.frame中。它包含变量ha，ok和ck，给出了随后在每个CK水平发生或未发生过心脏病发作的患者数量。将观察到的比例与CK水平进行作图是有意义的。结果图是图2.5。描述这些比例的特别方便的模型是E（M=i+辭+内,其中pi是在CK水平xi下心脏病发作的比例。这条曲线是S形的

肌酸激酶水平

图2.5随后被诊断为心脏病发作的观察到的患者比例，与进入时的CK水平比较。形状，边界从0到1。（显然，心脏数据并未显示此拟建曲线的较低尾部）。这意味着心脏病患者的预期数目由下式给出：£00+01如血三^（以弘）=1+e0o+0i血",g（他）=log其中Ni是已知的每个CK水平下的患者总数。这个模型在它的参数上有些是非线性的，但是如果'logitg（他）=logNi-血丿，isappliedtoitweobtaing（^i）=P0+P1xi,在模型参数中是线性的。logit链接是二项模型的典范链接，因此是R中的默认链接。在R中有两种用glm指定二项模型的方法。反应变量可以是观察到的成功二项试验的比例，在这种情况下，提供试验次数的阵列必须作为glm的权重参数提供。对于二进制数据，不需要提供权重向量，因为默认权重为1就足够了。反应变量可以作为两列阵列提供，其中第一列给出了二项“成功”的数目，第二列是二项“失败”的数目。对于当前的例子，将使用第二种方法。提供2个阵列PredictedvaluesTheoreticalQuantilesPredictedvaluesObs.numberPredictedvaluesTheoreticalQuantilesPredictedvaluesObs.number图2.6首次尝试拟合CK数据的模型检查图。的模型公式涉及使用cbind。这是一个拟合心脏病发作模型的glm呼叫因为logit链接对于二项式是标准的，因此是R的默认值。以下是有关该模型的默认信息：零偏差是只有一个常数项模型的偏差，而偏差残差是拟合模型的偏差（以及在二项模型情况下的调整偏差）。这些可以结合起来给出比例偏差的解释，即r2的普遍化，如下所示：AIC是该模型的Akaike信息准则注意到X210随机变量的偏差非常高，如果模型拟合良好，它应该近似。事实表明X210随机变量的概率非常小，为36.93。残差图（如图2.6所示）也表明拟合效果差。这些图与普通线性模型的，预测值是在线性预测因子的比例上而不是反应变量上，有些在正常qq图的直线关系偏离往往是可以预料的。如果数据太少，这些图很难解释，但似乎有一个趋势，即对拟合值进行残差绘制的方法，这会引起关注。此外，第一点有相当大的影响力。请注意，二进制数据对残差的解释将更加困难：练习2探索了可以在二进制情况下采用的简单方法。注意这些问题在原始估计概率上叠加的拟合值的图中没有那么明显的显现出来（见图2.7）：还要注意，由glm提供的二项式模型的拟合值是估计的pi值，而不是估计的pi值。残差图建议尝试一个立体线性预测器，而不是最初的直线。肌酸激酶水平图2.7心脏病发作对CK水平的预测和观察概率。很明显，4.252与x28分布的一致性并不太大（实际上它低于预期值），而且AIC已经大幅改善。残差图（图2.8）现在显示的模式不如以前的模型清晰，但如果我们有更多的数据，则偏离常数方差就会引起关注。此外，现在的拟合更接近数据（见图2.9）：我们还可以用R来测试mod.0对于需要mod.2的替代方法是否正确的原假设。虽然它是对正在执行的偏差（即广义似然比检验）的分析，而不是对方差的分析，但是有点令人困惑的是，使用了方差分析函数来做到这一点。PredictedvaluesTheoreticalQuantilesPredictedvaluesObs.number图2.8第二次尝试拟合CK数据的模型检查图。图2.9心脏病发作对CK水平的预测和观察概率。图2.10比利时每年的艾滋病病例这个低的p值表示反对原假设的非常有力的证据-我们确实需要模型2。回想一下，模型的比较比个别偏差的检验有更稳固的理论基础。2.3.2泊松回归流行病模型本章的介绍包括了一个流行病早期阶段的简单模型。Venables和Ripley(2003)提供了1981年以后比利时每年新发艾滋病病例数的一些数据。数据可以输入到R中并绘制如下。图2.10显示了结果图。与这些数据有关的科学趣味问题是，它们是否提供任何证据来表明新病例生成的基础比率增加正在放缓。介绍中的简单模型可能为开始调查此问题提供了一个合理的模型。该模型假设每年潜在预期案例数量m根据以下公式增加：pi=cexp(bti)其中c和b是未知参数，ti是自开始以来的年数PredictedvaluesTheoreticalQuantilesPredictedvaluesObs.number图2.11AIDS数据的mo拟合残差图。数据。对数链接将其转换为GLM，log(Pi)=log(c)+bti=P0+tiP1,我们假设yi?Poi(pi)，其中yi是第ti年观察到的新病例数量。yi被假定为独立的。这本质上是一种无检查的疾病传播模型。以下拟合该模型(对数链接对于泊松分布是典范的，因此也是R的默认值)并检查它。这个偏差对于它应该的X211随机变量的观测值是非常高的PredictedvaluesTheoreticalQuantilesPredictedvaluesObs.number图2.12AIDS数据的m1拟合残差图如果模型拟合得好，可以近似。图2.11所示的残差图也令人担忧。特别是残差均值的清晰模式(与拟合值相对作图)显示违反了独立性假设，可能是由于遗漏了某些模型中的重要因素。因为对于这个模型，拟合值随时间单调递增，如果残差相对于时间绘制，我们将得到相同类型的模式-即可以有效地将时间的二次项添加到模型中。库克距离图中显示的最后一年数据的非常大的影响力也令人担忧。注意，如果泊松均值较低，则残差图的解释可能会变得困难，因此数据大多为零和一。在这种情况下，如果适应泊松案例，练习2中所涉及的模拟方法可以证明是有用的。通过增加一个二次项来修正模型似乎是明智的：□i=exp（%+P1ti+卩2t2i）.该模型允许表示疾病无限制传播以外的情况。以下拟合并检查它：注意图2.12中显示的残差图现在有了很大的改进：均值的明显趋势消失了，残差的（垂直）扩散相当均匀，点13的影响大大减少，QQ图更直。此模型显示的更完整的模型总结也有所改善：现在的偏差非常合理，即接近x210r.v.的预期，而且AIC大幅下降。总之，这个模型看起来很合理。还要注意，glm摘要的结构如何与lm摘要的结构相似。值的大样本分布。z值列仅简单地报告参数估计值除以其估计的标准偏差。由于泊松不需要离散参数估计，如果相应参数的真值为零（至少在大样本限制下），这些z值应该是N（0,1）rvs的观测值，并且报告p值是基于这种分布近似。对于mod.1,报告的p值非常低：即对于每个参数都有清楚的证据表明它不是零。注意，Fisher评分迭代是GLM背景下IRLS迭代的另一个名称。对系数汇总表的检查表明卩2=0的假设可以被坚决拒绝，这提供了明确的证据表明mod.1比mod.0更可取。同样的问题也可以使用广义似然比检验来解决：结论与以前一样：p值很小表示mod.0应该被mod.1坚决拒绝。注意，总结中的p值和偏差表的分析是不同的，因为它们基于根本不同的近似分布结果。对方差分析的测试=“Chisq”参数是有道理的，因为这个模型的尺度参数已知，如果估计时将测试设置为“F”会更可取。在这里，模型选择的假设检验方法是适当的，因为感兴趣的主要问题是从这些数据中是否有证据表明流行病的传播尚未得到检查。如果证据不够确定，则事情就会好转，这是明智的。如果我们对简单地找到预测数据的最佳模型更感兴趣，则AIC的比较会更合适，但对这些数据的结论也一样。参数卩1可以解释为在流行病开始时疾病的传播速度：即，作为没有采取控制措施的新人群的疾病内在增长率。注意这个参数的估计在第一个和第二个模型之间实际上已经大大增加了：如果我们遇到了第一个拟合不好的模型，则摘要函数可轻松提取所需的估计值和标准误差，如下所示：使用标准正态分布的临界点是恰当的，因为这个模型的尺度参数已知。如果它已经被估计，那么我们将不得不使用来自t分布的临界点，将自由度设置为模型的残余自由度（即数据的数量少于估计的卩参数的数量）。另一个想要做的明显的事情是，使用该模型随时为案例生成的基础速率找到置信区间。以下R代码说明了如何使用predict.glm函数在整个数据周期中找到基础速率的CI,并绘制它们。Year图2.13根据模型ml显示的作为连续曲线的基础艾滋病病例率，95%的置信区间

显示为虚线。该图显示在图2.13中。注意，默认情况下，predict.glm函数根据线性预测器的比例来预测：我们必须应用链接函数的反函数来返回原始反应比例。因此，这些数据提供了非常确凿的证据表明，无拘束的指数增长模型过于悲观：在数据的最后，有充分证据表明增长率正在放缓。当然这个模型不包含任何机制内容-它没有说明减速过程如何或为什么会发生，因此它完全不适合超出数据范围的预测。该模型使我们有理由相信，新案例增长率的明显放缓是真实的，而不仅仅是机会变化的结果，但它对稍后可能发生的事情几乎没有或完全没有显示。下表根据他们对来世的信仰将女性和男性的样本分类：BelieverNon-BelieverFemale435147Male375134这些数据（在Agresti被报道，1996年）来自美国综合社会调查（1991）,而非信徒”类别包括“未决定”。在坚持这种信仰时，男性和女性之间是否存在差异？我们可以通过使用偏差分析来比较两种竞争模型对这些数据的拟合程度来解决这个问题：其中一个是信仰作为与性别无关的因素被建模，另一个是信仰与性别之间有一定的相互作用。首先考虑独立的模型。如果yi是对表格中一个单元格的计数的观测值，则我们可以将预期的计数数目建模为：如果yi是性别k的数据，则pi三E（Yi）=nYkaj，并且信仰j其中n是被调查的总人数，al是信徒的比例，a2是非信徒的比例，丫1和丫2分别是女性和男性的比例。取这个模型的对数ni三log®)=log(n)+log(Yk)+log(aj).Sodefiningn?=log(n),Y?k=log(Yk)anda?j=logg)themodelcanbewrittenas所以定义n?=log(n)，Y?k=log(yk)和a?j=log(aj)模型可写成这显然是GLM这显然是GLM结构，但显然不可识别。下降丫?1和a?l解决了可识别性问题的产生一m~一m~02隔-可4_01001110n72

鱼注意，在这个模型中，性别和信仰是两个层次的因素变量。如果列联表中的计数是随机独立发生的，则使用的明显分布将是泊松分布。事实上，即使表中的事件总数或甚至其他一些边际总数是固定的，则可以证明，正确的可能性可以写成泊松的乘积，在各种不同的固定量的条件下。因此，假设拟合模型被迫匹配固定总数和任何固定的边际总数，泊松分布仍然可以某些固定总数相匹配，只不过是坚持某些条款保留在模型中。简单的'独立'模型很容易在R中估计。首先输入数据并检查它：由于性别和信仰都是因素变量，因此模型设定非常简单。以下拟合模型并检查模型矩阵是否符合预期：现在看看拟合模型对象mod.0这个拟合看起来非常接近，如果一个信仰和性别之间相互作用的模型显着改善，那将会有些令人惊讶。从来没有这样一个模型可能是：ni三log®)=n?+Y?k+a?j+Z?kjifyiisdataforgenderkandbeliefjZkj是一个'交互参数'。这个模型允许信仰和性别的每个组合独立变化。正如所写，该模型有相当多的不可识别的项。n????Y?1????nl234????????110010100????????????a??????????????12Y*?2?n=110101000a?n101010001???z-z11n101100010z?12Z?21?Z?22但是这很容易减少到可以识别的地方：1010n100072耳31111a2耳4_1100_(22以下拟合模型，检查模型矩阵并输出拟合模型对象：为了测试是否存在性别与信仰之间的相互作用的证据，使用偏差分析，对mod.0正确的原假设与mod.1正确的更一般的替代方法进行了测试。0.69的p值表明没有证据可以拒绝模型0以及来世性别和信仰之间没有关联的假设。注意，事实上，具有交互作用的模型是饱和模型，这就是为什么它的偏差在数值上为零，并且没有真正需要拟合它并明确地将其与独立模型进行比较-在这种情况下，我们也可以检查独立模型的偏差。然而，对于这种简单的2路列联表所采取的一般方法可以很容易地推广到多路表和任意数量的组。换句话说，这里概述的方法可以扩展为使用对数线性GLM来分析分类数据的相当普遍的方法。最后，注意mod.0的拟合值具有奇特性质，尽管拟合值和原始数据不同，但男性和女性的总数在数据和拟合值之间保存，信徒和非信徒的总数也是如此。这是由于对数链