广义线性模型.doc

广义线性模型广义线性模型*（Nelder和Wedderburn，1972）除了正态分布，也允许反应分布，以及模型结构中的一定程度的非线性。GLM具有基本结构g(i)=Xi,其中iE（Yi），g是光滑单调链接函数，Xi是模型矩阵的第i行，X和是未知参数的向量。此外，GLM通常会做出Yi是独立的和Yi服从一些指数族分布的假设。指数族分布包括许多对实际建模有用的分布，如泊松分布，二项分布，伽马分布和正态分布。GLM的综合参考文献是McCullagh和Nelder（1989），而Dobson（2001）提供了一个全面的介绍。因为广义线性模型是以“线性预测器”X的形式详细说明的，所以线性模型的许多一般想法和概念通过一些修改而继续存在到广义线性模型中。除了必须选择的链接函数和分布之外，基本模型公式与线性模型公式基本相同。当然，如果恒等函数被选择作为链接以及正态分布，那么普通线性模型将作为特例被恢复。然而，泛化是以某种成本为代价的：现在的模型拟合必须要迭代完成，而且用于推理的分布结果是近似的，并且由大样本限制结果证明是正确的而不是精确的。但在深入探讨这些问题之前，请考虑几个简单的例子。i=cexp(bti),例1：在疾病流行的早期阶段，新病例的发生率通常会随着时间以指数方式增加。因此，如果i是第ti天的新病例的预期数量，则该形式的模型为请注意，“广义”和“一般”线性模型之间存在区别-后一个术语有时用于指除简单直线以外的所有线性模型。可能是合适的，其中c和b是未知参数。通过使用对数链路，这样的模型可以变成GLM形式log(i)=log(c)+bti=0+ti1（根据0=logc和1=b的定义）。请注意，模型的右侧现在在参数中是线性的。反应变量是每天新病例的数量，因为这是一个计数，所以泊松分布可能是一个合理的可以尝试的分布。因此，针对这种情况的GLM使用泊松反应分布，对数链路和线性预测器0+ti1。,例2：狩猎动物捕获猎物的速度yi往往随着猎物密度xi的增加而增加，但最终会趋于平衡，当捕食者捕获尽可能多的猎物时。对于这种情况一个合适的模型可能是其中a是未知参数，表示最大捕获率，h是未知参数，表示捕获率为最大速率一半时的猎物密度。很显然，这个模型在其参数中是非线性的，但是通过使用倒数链路，右边的参数可以是线性的：（这里01/a和1h/a）。在这种情况下，猎物捕获率的标准差可能与平均速率大致成比例，建议使用Gamma分布作为反应，并完成模型设定。们不限于示例的简单直线形式，但可以有对于线性模型可能具有的线性预测器的任何结构。2.1GLMs的理论GLM的估计和推理基于最大似然估计理论，尽管可能性的最大化需要迭代最小二乘法，与1.8.6节的方法相关。本部分首先介绍指数分布族，它允许发展一种通用方法来最大化GLM的可能性。然后基于似然理论的一般结果（在本章末尾推导）来讨论GLM的推论。在本节中，区分反应数据y和Y的观察值的随机变量有时很有用，所以它们在符号上具有区别性：对于估计和估计值，它尚未完成。2.1.1指数族分布f(y)=expyb()/a()+c(y,),GLM中的反应变量可以来自指数族的任何分布。如果一个分布的概率密度函数或概率质量函数可以写成，则该分布属于指数族分布其中b，a和c是任意函数，是任意的“尺度”参数，并且被称为分布的“典范参数”（在GLM上下文中，将完全依赖于模型参数，但是目前没有必要做这个明确）。例如，很容易看出，正态分布是指数族的一员，因为,.表2.1给出了在R中为GLM实施的指数族成员的类似分解。用a，b和可以得到指数族分布的均值和方差的一般表达式。给定一个特定的y，的对数似然性仅仅只是将logf（y）视为的一个函数。那是.将l作为一个随机变量来处理，通过用随机变量Y替换特定的观测值y，可以评估l/的期望值：使用E（l/）=0这个一般结果，（在取真值时，参见2.4节中的（2.14）和重新排列意味着E（Y）=b0（）。（2.1）即任何指数族随机变量的均值由bw.r.t.的一阶导数给出。，其中b的形式取决于特定的分布。该等式是将GLM的模型参数与指数族的典范参数联系起来的关键。在GLM中，参数决定了反应变量的均值，并且通过（2.1），它们决定了每个反应观测值的典范参数。,再次对似然性微分处理,并将其插入到一般结果中，E（2l/2）=-E（l/）2（衍生物在真实值下计算，参见结果（2.16），第2.4节），重新安排第二个有用的一般结果：var(Y)=b00()a().a原则上可以是的任何函数，并且当与GLM一起工作时，如果是已知的，处理任何形式的a都是没有困难的。然而，当未知，事情就会变得很尴尬，除非我们可以写出（）=/，其中是一个已知常数。事实上，这种限制形式涵盖了所有有实际意义的案例（见表2.1）。a（）=/允许基于正态分布的模型中的不等方差的可能性，但是在大多数情况下，仅为1.因此，我们现在有var(Y)=b00()/.(2.2)在随后的章节中，我们通常会将var（Y）视为E（Y）的函数，这会很方便，并且由于和通过（2.1）链接，我们总能定义一个函数V（）=b00（）/，使得var（Y）=V（）。表2.1列出了几个这样的功能。2.1.2拟合广义线性模型回想一下，GLM模拟独立反应变量的n阶向量Y，其中E（Y），通过g(i)=Xi和Yifi(yi),其中fi（yi）表示指数族分布，典范参数i由i（通过方程2.1）决定，因此最终由决定。给定Y的一个观测值向量y，的最大似然估计是可能的。由于Yi是相互独立的，的似然函数是nL()=Yfi(yi),i=1因此的对数似然函数是,右边对的依赖是通过i对的依赖。请注意，函数a，b和c可能随着i而变化-例如，允许不同的二项分母ni，对于每个二项式反应的观测值，或对于正常反应的不同方差（但在常数内已知）。另一方面，对于所有的i，假设是相同的。正如前一节所讨论的那样，对于实际工作来说，只考虑可以写出ai（）=/i的情况就足够了，其中i是已知常数（通常为1），在这种情况下.通过偏导lw.r.t来最大化收益。的每个元素，令结果表达式为零并求出。,通过这个链式法则,所以微分（2.1），我们可以得到,然后推出.将（2.1）和（2.2）代入最后一个方程，意味着求解的方程是(2.3),(2.4)然而，如果权重V（i）事先已知且与独立，那么这些方程就正好是为了通过非线性加权最小二乘寻找而必须解出的方程。在这种情况下，最小二乘的目标是其中i非线性地依赖于，但权重V（i）被视为固定的。要找到最小平方估计值，包括求解S/j=0j，但当V（i）项被视为固定时，这个方程组很容易被看作是（2.3）。这种对应立即提出了一种求解（2.3）的迭代方法。令k表示在第k次迭代处的估计参数向量，并且令k和k为具有元素ik=Xik和ik=g-1（ik）的向量，其中g-1（）是链路的反函数。从参数估计开始，0，迭代以下步骤，直到k的序列收敛：1.计算当前k隐含的V（ik）项。2.为了得到k+1（V（ik）被视为固定的而不是作为的函数），给出这些估计，使用第1.8.6节的方法来最小化（2.4）。3.设置k到k+1实际上，这种方法比需要的要慢。第2步本身涉及迭代，但在已经收敛之前实际上将非线性最小二乘方法迭代到收敛没有多大意义。因此，第2步通常被替换为：2.使用k作为初始值，仅执行第1.8.6节中给出的求解（2.4）的迭代方法的一次迭代来获得k+1。应用这种方法会产生一个相当紧凑和简洁的方案。为了看到这个，让我们以矩阵形式写出非线性最小二乘问题。定义对角矩阵Vk，其中Vkii=V（ik）(2.4)按照第1.8.6节的方法，被它在k附近的一阶泰勒展开替代，所以因此，没有进一步的近似Jij=i/j|k.Now.由“伪数据”的定义对角线权重矩阵.因此，下面的步骤可迭代到收敛1.使用当前的k和k计算伪数据zk和迭代权重Wk。2.将对于的平方和最小化，以便获得k+1，因此k+1=Xk+1和k+1。将k增加1。收敛的解决了（2.3），因此是的最大似然估计。该算法在大多数实际情况下趋于一致，但也有例外（例如，二项数据的不良或过于灵活的模型）。请注意，要开始迭代，我们只需要0和0的值，不需要0。因此，迭代通常通过设置0i=yi和i0=g（0i）开始，并根据需要轻微调整0i，以避免无限的i0（例如，如果yi=0且有对数链路）。该方法被称为IterativelyRe-weightedLeastSquares（IRLS），原因很明显，在此背景下，归因于Nelder和Wedderburn（1972）。2.1.3 IRLS的目标是对数似然的二次近似IRLS迭代中的工作线性模型不仅仅是寻找参数的最大似然估计的手段。在一个加性常数内（在收敛时）也是模型在附近的对数似然的二次近似。显然，第一个衍生工具是关于对数似然和S之间的j匹配：实际上它们都是零。S的二阶导数矩阵为-XWX/，并且这被证明与对数似然的期望的二阶导数矩阵匹配，并且因此在大样本极限中由大数定律来匹配二阶导数矩阵本身。为了证明这一点，首先将u定义为关于模型参数的对数似然的导数向量，因此ui=l/i，然后将（2.3）中的导数以矩阵向量形式重写为u=XTG1V1(y)/.注意，如果V（i）被视为的函数，算法不会最小化（2.4），因为在这种情况下，令导数为零不会产生（2.3）。换句话说，最大似然与具有平均方差关系的最小二乘基本上不同。E(uuT)=XTG1V1E(Y)(Y)TV1G1X/2=XTG1V1VV1G1X/=XTWX/因为E（Y-T）（Y-）T=V。根据一般似然结果（2.19），在2.4.2节中，-E（uu）也是对数似然的期望二阶导数矩阵。导数的这种对应关系足以证明S是附近的对数似然的二次近似，并且由于MLE的一致性，它们在真实参数值附近。2.1.4 AICforGLMs通过对可能性进行直接比较的模型选择存在如下问题：如果将冗余的参数添加到正确的模型中，可能性几乎总是增加（并且从不减小），因为额外的参数让模型更接近数据，即使这只意味着对数据的组成部分进行“噪声建模”。正如在线性模型的情况下，如果我们能够根据它们拟合数据均值而不是数据y的能力来选择模型，这个问题会得到缓解。在GLM背景下，合理的方法是根据模型最大化l（;）而不是l（;y）的能力来选择模型，但是为了做到这一点，我们必须能够估计出l（;）。实际上这个估计很明确。从2.1.3节我们已经有,并且因为y=时，这也必须成立.然后，这个论据导致在第1.8.5节的（1.15）（只修改权重）产生了估计量tr(A)l(;y)tr(A)+n/2A=X(XTWX)1XTWandhencetr(A)=p,（可识别的）模型参数的数量。因此，在模型之间进行选择时，我们会选择哪个模型具有最高的l（）-p值，这相当于选择了Akaike信息准则的最低值（Akaike，1973），AIC=2-l（）+p。前面的论点假设已知。如果不是，那么就需要一个的估计来计算AIC，结果AIC中的惩罚项p将变为p+1。这种归纳在2.4.7节中是合理的。2.1.5的大样本分布GLM的分布结果并不精确，反而是基于大样本近似，利用包括一致性在内的最大似然估计的一般性质（见2.4节）。根据最大似然估计量的一般性质，我们认为，在大样本限制下，N(,I1),其中I=E（uuT）是模型参数的信息矩阵，u是关于模型参数的对数似然函数的导数向量，因此ui=l/i（见2.4节（2.20）和（2.19）。在第2.1.3节中显示E（uT）=XTWX/，因此在大样本限制中N(,(XTWX)1).对于具有已知尺度参数的分布，可以直接使用此结果来查找参数的置信区间，但是如果尺度参数未知（例如对于正态分布），则必须估计它，并且区间必须基于一个合适的t分布。尺度参数估计在第2.1.7节中介绍。2.1.6通过假设检验比较模型考虑检验H0:g()=X00相反H1:g()=X11,其中是反应向量Y的期望，Y的元素是来自指数族分布的相同成员的独立随机变量，其中X0X1。如果我们有反应向量的观测值y，则可以执行广义似然比检验。令l（0）和l（1）为两个模型的最大似然率。如果H0为真，则在大样本限制中，2l(1)l(0)2p1p0,(2.5)sothattwicethewhichshouldstrictlybeamaximumlikelihoodestimate,其中pi是模型i中（可识别的）参数（i）的数量（参见2.4.5节和2.4.6节推导该结果）。如果原假设是假的，则模型1倾向于具有比模型0高得多的可能性，从而两倍于严格应该是最大似然估计的，或者在大样本限制下倾向于MLE的估计值。对数似然对于相关的2分布的一致性差异太大。如果可以计算相关模型的对数似然性，那么近似结果（2.5）仅仅是有用的。在由IRLS估算的GLM情况下，只有当尺度参数已知时才是这种情况。因此，结果可以与泊松和二项模型一起直接使用，但不能与正态，伽马或反高斯分布（其中尺度参数未知）一起使用。稍后将简短地讨论在后面这些情况下做什么。偏差在实践中使用GLM时，在普通的线性建模中，用与残差平方和类似的方式来解释数量是有用的。D=2l(max)l()n2iyi(ii)b(i)+b(i),(2.6)(2.7)Xi=1hi这个数量是模型的偏差，被定义为thetermsinsidethesummation其中l（max）表示饱和模型的最大可能性：每个数据点具有一个参数的模型。在给定数据的情况下，l（max）是可能性具有的最大值，并且可以通过简单地设置=y并计算出可能性来计算。和分别表示关于饱和模型和兴趣模型的典范参数的最大似然估计。注意如何将偏差定义为与独立。表2.1列出了单个数据对偏差的贡献，对于若干分布-这些是偏差定义中求和中的项（术语）。scaleddeviance,与偏差相关的是调整偏差，D=D/,这取决于尺度参数。对于二项分布和泊松分布，其中=1，偏差和比例偏差是相同的，但通常情况并非如此。通过广义似然比检验结果（2.5），我们可以预计，如果模型是正确的，那么近似D2np,(2.8)tojustify(2.8)asalargesampleapproximationundermanycircumstances在大样本限制中。实际上，这样的论点是捏造的，因为证明（2.5）的限制论证依赖于模型中参数的数量保持不变，而样本容量趋于无穷大，但饱和模型具有与数据一样多的参数。渐近结果可用于表2.1中的一些分布，以在许多情况下将（2.8）证明为大样本近似值（参见当然，对于正态分布和一致性链接，我们使用第1章的结果。这对于正态分布情形下是精确的。但是，请注意，它完全打破了二进制数据的二项式。鉴于偏差的定义，很容易看出，本节开始的似然比检验可以通过重新表达两次对数似然比统计量作为。然后在H0下(2.9)（在大样本极限内），其中Di*是具有pi个可识别参数的模型i的偏差。但是，这只有在尺度参数已知时才有用，以便可以计算D*。与未知的模型比较在H0下我们有近似的结果and,并且，f如果被看作是渐近独立的，这就意味着,在大样本限制下（当然，在普通线性模型特例中结果确实如此）。F的有用性质是它可以在不知道的情况下进行计算，这可以从比率收益的顶部和底部取消，在H0下，近似结果.(2.10)这个结果的优点是，当未知时，它可用于基于模型比较的假设检验。缺点是对于Di的可疑的分布假设以及它所基于的独立性近似。当然，一个明显的替代方法是使用估计值来获得每个模型的估计值Di*=Di，然后使用（2.9）进行假设检验。然而，如果我们使用估计（2.11）来达到这个目的，则很容易看出它只是（n-p1）F，所以我们的测试将完全等同于使用F比率结果（2.10），但是使用Fp1-p0，作为参考分布。显然直接使用（2.10）是一种更保守的方法，因此通常是首选：它至少在估计尺度参数时考虑到了不确定性。2.1.7andPearsonsstatistic正如我们已经看到的那样，参数的MLE可以在不知道尺度参数的情况下获得，但是在这些参数未知的情况下，通常必须进行估计。近似结果（2.8）提供了一个明显的估计量。一个2n-p随机变量的期望是n-p，所以等同于观测值达到我们的近似预期值D=D/(np).(2.11)第二个估计量基于Pearson统计量，定义为显然，X2/是一组零均值，单位方差，随机变量的平方和，具有n-p个自由度，这表明如果模型适当，则近似X2/x2n-p：这个近似值是有根据的。将观测到的Pearson统计量设置为我们得到的预期值=X2/(np).请注意，它很明确的表明这一点,其中W和z是IRLS权重和伪数据，在收敛时进行评估。2.1.8典范链接函数典范链接gc用于分布，是链接函数，使得gc（i）=i，其中i是分布的典范参数。例如，对于泊松分布，典范链接是对数函数（其他示例见表2.1）。典范链接的使用意味着i=Xi（其中Xi是X的第i行）。典范链接倾向于具有一些很好的性质，例如确保保持在反应变量的范围内，但它们也具有更多微妙精细的优点，其中之一在此处被推导出。回想一下，似然最大化涉及对每个j的对数似然性微分，并令结果为零，以获得方程组但是，如果正在使用典范链接，则i/j=Xij，如果像通常情况，wi=1i，这个方程组可简化为XTyXT=0,i.e.toXTy=XT.即，简化为现在考虑X包含1列的情况：任意其他加权求和，其中权重由模型矩阵P给出，这意味着该系统中的一个方程简单地为iyi=ii。RecallthatifZi:i=1.nareasetofi.i.d.N(0,1)r.v.sthenPZi22n.回想一下，如果Zi：i=1.n是一组独立同分布，N（0,1）r.v.然后是PZi2-2n。列（或这些的线性组合）在原始数据和拟合值之间保存。其中一个实际结果是，对于任何具有截距项和典型联系的GLM，残差总和为零：这种“观察到的无偏性”是一个令人放心的特性。结果的另一个实际用途是在使用对数线性模型的分类数据分析中，其中它通过规定模型提供了一种方法，可以保存任何模型中保存的内置于研究设计中的总数。2.1.9残差模型检查可能是应用统计建模中最重要的部分。在普通线性模型情况下，这是基于对模型残差的检验，模型残差包含数据中的所有信息，而不是由模型的系统部分来解释。残差检验也是GLM案例模型检验的主要手段，但在这种情况下，残差标准化是必要的，并且更加困难。对于GLM而言，不仅仅只检查原始残差的主要原因是难以检查原始残差检验假设的均值方差关系的有效性。例如，如果采用泊松模型，则残差的方差应与拟合值（i）的大小成正比增加。然而，如果将原始残差与拟合值作图，则需要一种杰出的能力来判断残差的可变性是否与均值成比例增加，而不是均值的平方根或平方。出于这个原因，通常将GLM残差标准化，以这种方式，如果模型假设正确，则标准化残差应该具有大致相等的方差，并且尽可能与来自普通线性模型的残差相似（尽管见第6.5节中图6.9的替代绘图方法）。皮尔森残差根据拟合模型，将残差标准化的最明显的方法是将它们除以与其标准差成比例的量。这导致了皮尔森残差,如果模型正确，它应该具有近似为零的均值和方差。当与拟合值或任何协变量（无论是否包含在模型中）进行绘图时，这些残差不应显示均值或方差的任何趋势。“皮尔森残差”这个名字与相关的皮尔森残差的平方和给出了2.1.7节中讨论的皮尔森统计量。请注意，皮尔森残差是来自收敛IRLS方法的工作线性模型的残差除以收敛IRLS权重的平方根。异常残差在实践中，皮尔森残差的分布在零附近可以是非常不对称的，因此它们的行为不像可能希望的那样接近普通线性模型残差。异常残差在这方面通常是优选的。异常残差是通过注意到偏差对GLM起到与残差平方和对普通线性模型起相同作用来得到：确实，对于普通线性模型来说，偏差是残差平方和。在普通线性模型情况下，偏差由平方残差和组成。这就是具有适当标示的偏差组成部分的平方根的残差。因此，把di写成第i个数据所贡献的偏差的一部分（即（2.7）中求和的第i项），我们有并且通过与普通线性模型类比，我们可以定义.根据需要，这些“异常残差”的平方和给出了偏差本身。现在如果计算了所有参数已知的模型的偏差，则（2.8）将变成D*x2n，这可能表明对于单个数据dix12，可推出。当然，（2.8）不能合理地应用于单个数据，但对于良好的拟合模型来说，这表明我们可能预期异常残差具有类似于N（0,1）随机变量的行为，特别是在对（2.8）预计是一个合理的近似值的情况下。2.1.10拟极大似然迄今为止，GLM的处理假定反应变量的分布是指数族的已知成员。如果有充分的理由假设该反应变量遵循特定的分布，则它对该分布上的基本模型很有吸引力，但是在许多情况下，反应分布的性质并不十分清楚，并且可能它只能详细说明反应变量的方差与其均值之间的关系。也就是说，函数V（）可以被指定，但是好不了多少。然后出现的问题是，是否有可能开发用于拟合和推理GLM的理论，从仅指定平均方差关系的位置开始。事实证明，基于拟极大似然的概念，开发出令人满意的方法是可能的。考虑一个具有均值i和已知方差函数V（i）的随机变量的观测值yi。然后对于给定yi的i的对数拟极大似然定义为(2.12)正如我们将会看到的，这个函数的关键特征是它具有li的许多有用特性，即对应于单个观测的对数似然，但只需要V的知识而不是Yi的全部分布。假设数据是独立随机变量的观测值，我们可以定义所有反应数据的平均向量的对数拟极大似然或定义的任意参数向量为.q的关键特征是，为了推导GLM，它表现得与对数似然函数非常相似，但只需要方差函数就可定义它。例如，考虑获得GLM参数的最大拟似然参数估计。对q关于j产量微分,所以参数估计是方程的解但这正是方程组（2.3），必须解决这个问题才能找到GLM的m.l.e.s。因此，最大拟似然参数估计可以通过通常的GLMIRLS方法找到，在任何情况下只需要知道V（）。此外，对数拟似然性与参数估计值的大样本分布结果（参见2.1.5节）的对数似然性具有足够的性质，并且对于参数的最大拟似然估计量值也适用。类似地，当对数似然度l被对数拟似然函数q代替时，2.1.6节的大样本分布结果支持用GLM进行假设检验。第2.4.8节提供了这些断言的理论基础。请注意，饱和模型的对数拟似然始终为零，因此GLM的拟偏差很简单Dq=2q().很显然，关于残差和尺度参数估计的讨论也从似然性中延续到拟似然的情况，但不超过用q代替l。15拟似然性方法的实际应用要求对（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，相等线拟合将与它们遇到的模型流形的（切线空间）正交。3.对于拟合4个数据的含有2个参数的模型，相等线拟合将变成相等平面拟合。一般来说，图2.1所示的几何图形适用于任何GLM。有了更多的数据，相等线拟合就变成了n-p维相等面的拟合，其中n和p分别是数据和参数的数量：对于任何固定的，方程（2.3）给出了定义这样一个平面的y的限制。请注意这些平面1.20.00.20.40.60.81.00.20.30.40.50.60.70.8yx图2.1GLM的几何图形。左端板展示了所示的三个x，y数据的广义线性模型E（y）=exp（0+1x）的最佳拟合，假设每个yi是由模型给出的带有均值的Gamma分布随机变量的一个观测值。右端板展示了使用此模型作为示例的GLM拟合的几何图形。所示的单位立方体表示向量（y1，y2，y3）T定义单点空间。灰色表面根据模型显示所有可能的预测值（在单位立方体内），即它表示所有可能的（1，2，3）T值。由于参数0和1允许在所有可能的值上变化，这就是相应模型“拟合值”追踪的表面：“模型流形”。连续线从立方体的一个面开始，在另一个面上离开，是等效拟合线：位于这样一条线上的反应数据（y1，y2，y3）T的值各自导致相同的0，1的最大似然估计，即相同的（1，2，3）T。注意等效拟合线既不相互平行，也不与模型流形正交。可以相交-稍后将返回的点。对于离散反应数据，图片没有什么不同，尽管只有在可能性连续概括下（通常可以通过在概率函数中用合适的伽玛函数代替阶乘）才能严格得到相等线拟合。只有正态分布才能与模型流形有正交的相等线/平面拟合。对于其他分布，相等线/平面拟合有时可能彼此平行，但决不会全部与模型流形正交。2.2.1IRLE的几何图形通过考虑一个参数模型对2个反应数据的拟合，IRLS估计算法的几何图形是最容易理解的。图2.2说明了这种模型的几何结构：在这种情况下，一个具有对数链接和Gamma错误的GLM，2.01.50.00.51.08101214161820xy15051005101512图2.2GLM的几何图形。E（yi）i=20exp（-xi）其中yiGamma和i=1,2。左端板展示了拟合2个x，y数据的模型（连续线）的最大似然估计值，显示为。右端板显示了拟合几何图形。1515的正方形是空间2中的一部分，其中（y1，y2）定义为单点。粗体曲线是“模型流形”：它根据模型包含所有可能的点（1，2）（即，当变化时，（1，2）可以追踪该曲线）。细线是相等线拟合的例子。位于这些线之一上的所有点（y1，y2）共享相同的的MLE，即（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)1051500510152150510051015211(c)(d)1.50.00.51.02.02.53.00.00.51.01.52.02.53.023.00.50.01.01.52.02.50.00.51.01.52.02.53.0211图2.3基于图2.2所示的GLM的IRLS估计的几何图形。（a）显示拟合问题的几何图形-模型流形是粗黑色曲线，等拟合线是细线（如图2.2），数据是和拟合值的当前估计值k。（b）问题以当前拟合值为新中心（yi被yi-i代替）。（c）问题被重新线性调整，使得X的列现在跨越模型流形的切线空间。切线空间用灰线表示（这个步骤用g0（ik）（yi-ik）代替yi-ik。（d）该问题是线性转换的，因此现在由Xk给出位置。对于大多数GLM，现在必须通过将相对于每个轴的分量乘以Wi来再次调整问题，其中Wi是迭代权重：这将确保通过的相等估计线与切线空间正交。在当前的例子中，这些权重都是1，因此所需的正交性已经成立。现在对于转换后的问题，在模型流形的附近，可以用切线空间来近似，其中相等线拟合近似正交：因此，可以通过找到转换数据的最小二乘投影来获得和的更新估计，到切线空间（灰线）。0.60.40.00.20.81.00.00.20.40.60.81.012图2.4拟合和收敛问题的几何图形。具有对数链接和正常错误的含有1个参数的GLM的几何图形被显示。厚曲线是模型流形-根据模型，在单位平方内包含所有可能的数据拟合值。细线是相等线拟合（水平如图2.2）。注意在图的左上角，相等线拟合是如何相互交叉的。该重叠区域中的数据将产生具有多于一个参数值的局部最小值的模型可能性。考虑IRLS拟合方法的运行，表明在这种情况下，根据用于启动拟合过程的初始值，它可以收敛到不同的估计值。说明问题反应向量的位置，用于说明文本中的非独特收敛。2.2.2IRLS收敛的几何图形图2.4说明拟合模型的几何图形，E（yi）i=exp（-xi），其中yi是正态分布的，有两个数据yi拟合，其中x1=0.6和x2=1.5。如前两节所示，相等线拟合显示在一个图形上，其中反应向量（y1，y2）T将定义一个单点，并根据该模型，所有可能的拟合值（1，2）T的集合显示为一条厚曲线。在这个例子中，相等线拟合在图的左上角相交并交叉（对应于非常差的模型拟合）。这个交叉点是有问题的：特别是，IRLS拟合位于左上角的数据的结果将取决于IRLS处理开始的初始参数估计，因为每个这样的数据点位于两条相等拟合线的交点处。如果IRLS迭代从图形右上角的拟合值开始，则接近右上角的拟合值将被估计，然而用图形左下角的拟合值开始迭代将导致估计的拟合值不同，并更接近图形的左下角。在实践中这确实发生，可以很容易的在R中展示出来，通过拟合数据y1=.02，y2=.9，如图2.4所示，可以很容易地证明。请注意，这里的第二个拟合实际上具有较高的可能性（较低的偏差）-拟合在可能性方面不是等同的。引起这些模糊的拟合几何图形的类型并不总是会发生：例如，一些模型具有相等拟合的平行线/平面，但对于任何具有相等线拟合的交叉线/平面的模型而言，存在一些模糊性的范围。幸运的是，如果模型是一个好的模型，那么位于模糊区域的数据往往是不太可能的。在图2.4的例子中，问题区域完全由模型只能很差地拟合的数据组成。由此可见，非常差的数据模型可能会产生这种类型的估计问题：但对于非常差的模型而言，将模型作为任何复杂数据集的早期尝试的特征并不罕见。如果遇到这样的问题，那么通过对转换后的反应数据进行线性建模可能会更好，直到已经确定足够好的候选模型转换回GLM为止。当然，如果选择合理的起始值，那么在拟合GLM时拟合过程中的模糊性不太可能引起重大问题：毕竟，该算法将收敛到可能性的局部最小值之一。然而，当它在可选最小值之间循环而不会收敛时，模糊性可以通过“性能迭代”在GAM估计中引起更严重的收敛问题，。2.3GLMswithRglm函数提供了在R中使用GLM的手段。它的使用与lm函数的使用类似，但有两点区别。模型公式的右侧指定线性预测变量的形式，现在给出了反应均值的链接函数，而不是直接使用反应均值。此外，glm需要一个族论点，用于指定要使用的指数族的分布以及与之一起使用的链接函数。在本节中，将介绍glm函数与各种简单GLM的用法，以说明GLM涵盖的各种模型结构。2.3.1二项分布模型和心脏病如果对于心脏病患者能够提供最好的照顾，那么能够尽早发现心脏病是再好不过的。关于检测心脏病的一种方法是检测血液中肌酐激酶（CK)的水平。35CKvaluePatientswithPatientswithoutHeartattackheartattack202886013261003081403051802102201912601813001313401913801504207046080表2.2作为CK水平函数的心脏病发作概率的数据（来自Hand等人，1994）。血流。一项研究被进行（Smith，1967），其中测量了怀疑患有心脏病的360名患者的CK水平。后来经过更长时间的医学调查后，确定每名患者是否确实有心脏病发作。数据见表2.2。原始报告根据CK水平的范围对患者进行分类，但在表中仅给出了该范围的中点。能够基于这些数据的诊断标准是很好的，这样CK水平可以用来估计患者发生心脏病的概率。我们可以通过构建一个模型，试图从CK水平解释心脏病发作的患者比例，从而朝着这样一个目标迈进。在接下来的数据被读入一个名为heart的data.frame中。它包含变量ha，ok和ck，给出了随后在每个CK水平发生或未发生过心脏病发作的患者数量。将观察到的比例与CK水平进行作图是有意义的。结果图是图2.5。描述这些比例的特别方便的模型是,其中pi是在CK水平xi下心脏病发作的比例。这条曲线是S形的1002003004000.00.20.40.60.81.0ProportionHeartAttack肌酸激酶水平图2.5随后被诊断为心脏病发作的观察到的患者比例，与进入时的CK水平比较。形状，边界从0到1。（显然，心脏数据并未显示此拟建曲线的较低尾部）。这意味着心脏病患者的预期数目由下式给出：,其中Ni是已知的每个CK水平下的患者总数。这个模型在它的参数上有些是非线性的，但是如果logit链接被应用，我们可以得到，,isappliedtoitweobtaing(i)=0+1xi,在模型参数中是线性的。logit链接是二项模型的典范链接，因此是R中的默认链接。在R中有两种用glm指定二项模型的方法。1.反应变量可以是观察到的成功二项试验的比例，在这种情况下，提供试验次数的阵列必须作为glm的权重参数提供。对于二进制数据，不需要提供权重向量，因为默认权重为1就足够了。2.反应变量可以作为两列阵列提供，其中第一列给出了二项“成功”的数目，第二列是二项“失败”的数目。对于当前的例子，将使用第二种方法。提供2个阵列261012842031123ResidualsResidualsvsFitted1931.50.50.00.51.01.52010010Std.devianceresid.NormalQQplot139PredictedvaluesTheoreticalQuantiles4202681012012345Std.devianceresid.ScaleLocationplot1391224681002468CooksdistanceCooksdistanceplot135PredictedvaluesObs.number图2.6首次尝试拟合CK数据的模型检查图。的模型公式涉及使用cbind。这是一个拟合心脏病发作模型的glm呼叫因为logit链接对于二项式是标准的，因此是R的默认值。以下是有关该模型的默认信息：零偏差是只有一个常数项模型的偏差，而偏差残差是拟合模型的偏差（以及在二项模型情况下的调整偏差）。这些可以结合起来给出比例偏差的解释，即r2的普遍化，如下所示：AIC是该模型的Akaike信息准则，在2.1.4和2.4.7节中讨论过（它也可以使用AIC（mod）提取）。注意到210随机变量的偏差非常高，如果模型拟合良好，它应该近似。事实上表明210随机变量的概率非常小，为36.93。残差图（如图2.6所示）也表明拟合效果差。这些图与普通线性模型的模型检验图的解释相同，在1.5.1节中详细讨论过，除了它现在是绘制的偏差残差，预测值是在线性预测因子的比例上而不是反应变量上，有些在正常qq图的直线关系偏离往往是可以预料的。如果数据太少，这些图很难解释，但似乎有一个趋势，即对拟合值进行残差绘制的方法，这会引起关注。此外，第一点有相当大的影响力。请注意，二进制数据对残差的解释将更加困难：练习2探索了可以在二进制情况下采用的简单方法。注意这些问题在原始估计概率上叠加的拟合值的图中没有那么明显的显现出来（见图2.7）：还要注意，由glm提供的二项式模型的拟合值是估计的pi值，而不是估计的i值。残差图建议尝试一个立体线性预测器，而不是最初的直线。1002003004000.00.20.40.60.81.0ProportionHeartAttack肌酸激酶水平图2.7心脏病发作对CK水平的预测和观察概率。很明显，4.252与28分布的一致性并不太大（实际上它低于预期值），而且AIC已经大幅改善。残差图（图2.8）现在显示的模式不如以前的模型清晰，但如果我们有更多的数据，则偏离常数方差就会引起关注。此外，现在的拟合更接近数据（见图2.9）：我们还可以用R来测试mod.0对于需要mod.2的替代方法是否正确的原假设。 虽然它是对正在执行的偏差（即广义似然比检验）的分析，而不是对方差的分析，但是有点令人困惑的是，使用了方差分析函数来做到这一点。28402461.00.01.0ResidualsResiduals vs Fitted54101.51.50.50.51.0100510Std. deviance resid.Normal QQ plot5410Predicted valuesTheoretical Quantiles48622400.01.02.03.0Std. deviance resid.ScaleLocation plot5410246810120.000.150.30Cooks distanceCooks distance plot5410Predicted valuesObs. number图2.8第二次尝试拟合CK数据的模型检查图。1002003004000.00.20.40.60.81.0Proportion Heart Attack图2.9心脏病发作对CK水平的预测和观察概率。198219841986198819901992050150250New AIDS cases图2.10比利时每年的艾滋病病例这个低的p值表示反对原假设的非常有力的证据 - 我们确实需要模型2。回想一下，模型的比较比个别偏差的检验有更稳固的理论基础。2.3.2 泊松回归流行病模型本章的介绍包括了一个流行病早期阶段的简单模型。 Venables和Ripley（2003）提供了1981年以后比利时每年新发艾滋病病例数的一些数据。 数据可以输入到R中并