第五章 回归分析 前几章的方法都只涉及一种变量，主要是比较它的各组值之间的差异。但生物学所涉及的问题是多种多样的，对许多问题的研究需要考虑不只一个变量，例如生物的生长发育速度就与温度，营养，湿度……等许多因素有关，我们常常需要研究类似的多个变量之间的关系。这种关系可分为两大类，即相关关系与回归关系。 相关关系：两变量X，Y均为随机变量，任一变量的每一可能值都有另一变量的一个确定分布与之对应。 回归关系：X是非随机变量或随机变量，Y是随机变量，对X的每一确定值xi都有Y的一个确定分布与之对应。 从上述定义可看出相关关系中的两个变量地位是对称的，可以认为它们互为因果；而回归关系中则不是这样，我们常称回归关系中的X是自变量，而Y是因变量。即把X视为原因，而把Y视为结果。 这两种关系尽管有意义上的不同，分析所用的数学概念与推导过程也有所不同，但如果我们使用共同的标准即使y的残差平方和最小（最小二乘法，详见下述），则不管是回归关系还是相关关系都可以得到相同的参数估计式。因此本章将集中讨论数学处理较简单的回归关系，且X限定为非随机变量。从这些讨论中所得到的参数估计式也可用于X为随机变量的情况，但我们不再讨论X为随机变量时的证明与推导。 另外，回归分析和相关分析的目的也有所不同。回归分析研究的重点是建立X与Y之间的数学关系式，这种关系式常常用于预测，即知道一个新的X取值，然后预测在此情况下的Y的取值；而相关分析的重点则放在研究X与Y两个随机变量之间的共同变化规律，例如当X增大时Y如何变化，以及这种共变关系的强弱。由于这种研究目的的不同，有时也会引起标准和方法上的不同，我们将在相关分析一节中作进一步介绍。 从两个变量间相关（或回归）的程度来看，可分为以下三种情况： （1）完全相关。此时一个变量的值确定后，另一个变量的值就可通过某种公式求出来；即一个变量的值可由另一个变量所完全决定。这种情况在生物学研究中是不太多见的。 （2）不相关。变量之间完全没有任何关系。此时知道一个变量的值不能提供有关另一个变量的任何信息。 （3）统计相关（不完全相关）。介于上述两种情况之间。也就是说，知道一个变量的值通过某种公式就可以提供关于另一个变量一些信息，通常情况下是提供有关另一个变量的均值的信息。此时知道一个变量的取值并不能完全决定另一个变量的取值，但可或多或少地决定它的分布。这是科研中最常遇到的情况。本章讨论主要针对这种情况进行。为简化数学推导，本章中如无特别说明，一律假设X为非随机变量，即X只是一般数字，并不包含有随机误差。但所得结果可以推广到X为随机变量的情况。 按相关中涉及公式类型可把相关关系分为线性相关和非线性相关。在多数情况下，我们提到相关关系时都是指线性相关，这是因为线性相关的理论已经很完善，数学处理也很简单；而非线性问题则需要具体问题具体分析，常常没有什么好的解决方法，理论上能得到的结果也很有限（详见§5.4）。因此在一般情况下我们常常只能解决线性相关的问题。也正是因为如此，在不加说明的情况下提到相关时常常是指线性相关；如概率论基础部分曾提到独立可以推出不相关，而逆命题不成立。讨论回归关系时也有类似现象。 下面我们就来讨论回归关系中最简单的情况：一元线性回归。 §5.1 一元线性回归 前边已经说过，回归关系就是对每一个X的取值xi，都有Y的一个分布与之对应。在这种情况下，怎么建立X与Y的关系呢？一个比较直观的想法就是建立X与Y的分布的参数间的关系，首先是与Y的均值的关系。这就是条件均值的概念，记为：。它的意思是在X=x1的条件下，求Y的均值。更一般地，我们用代表X取一切值时，Y的均值所构成的集合。所谓一元线性回归，就是假定X与之间的关系是线性关系，而且满足：  （5.1） 此时进行回归分析的目标就是给出参数α和β的估计值。 例5.1 对大白鼠从出生第6天起，每三天称一次体重，直到第18天。数据见表5.1。试计算日龄X与体重Y之间的回归方程。 表5.1 大白鼠6-18日龄的体重 序号 1 2 3 4 5  日龄xi 6 9 12 15 18  体重yi 11 16.5 22 26 29  首先，我们可以把数对（xi, yi）标在X-Y坐标系中，这种图称为散点图。它的优点是可以使我们对X、Y之间的关系有一个直观的、整体上的印象，如它们是否有某种规律性，是接近一条直线还是一条曲线，等等。我们还可以画很多条接近这些点的直线或曲线，但这些线中的哪一条可以最好地代表X, Y之间的关系，就不是凭直观印象可以做出判断的了。例如对例5.1，我们可画出如下的散点图：  图5.1 大白鼠日龄—体重关系图 图中的点看来是呈直线关系，但那条直线是否最好地反映了这种关系呢？或者换一种说法：该如何找到最好地反映这种关系的直线呢？这就是我们以下要讨论的问题。 一元正态线性回归统计模型： 线性回归意味着条件平均数与X之间的关系是线性函数：  （5.1） 对于每个Y的观察值yi来说，由于条件均值由(5.1)式决定，观察值就应该是在条件均值的基础上再加上一个随机误差，即：  （5.2） 其中。正态线性回归中“正态”的意思是随机误差服从正态分布。(5.2)式就是一元正态线性回归的统计模型。 参数α和β的估计 统计模型中的α和β是总体参数，一般是不知道的。由于只能得到有限的观察数据，我们无法算出准确的α与β的值，只能求出它们的估计值a和b，并得到yi的估计值为：  （5.3） 那么，什么样的a和b是α和β最好的估计呢？换句话说，选取什么样的a和b可以最好地反映X和Y之间的关系呢？一个合理的想法是使残差最小。为了避免使正负ei互相抵消，同时又便于数学处理，我们定义使残差平方和达到最小的直线为回归线，即令： ，且  得：  整理后，得  （5.4） 上式称为正规方程。解此方程，得：  这种方法称为最小二乘法，它也适用于曲线回归，只要将线性模型（5.3）式换为非线性模型即可。但要注意非线性模型的正规方程一般比较复杂，有些情况下甚至没有解析解。另一方面，不管X与Y间的真实关系是什么样的，使用线性模型的最小二乘法的解总是存在的。因此正确选择模型很重要，而且用最小二乘法得出的结果一般应经过检验。 记 ，称为X的校正平方和； ，称为Y的总校正平方和； ，称为校正交叉乘积和， 则：  （5.7） 在实际计算时，可采用以下公式：  现在回到例5.1。 例5.1 对大白鼠从出生第6天起，每三天称一次体重，直到第18天。数据见表5.1。试计算日龄X与体重Y之间的回归方程。 表5.1 大白鼠6-18日龄的体重 序号 1 2 3 4 5  日龄xi 6 9 12 15 18  体重yi 11 16.5 22 26 29  解：把数据代入上述公式,得：      即：所求的回归方程为：y = 2.6996 + 1.5167 x 带有统计功能的计算器常常也可以做一元线性回归，对于这样的计算器，只需把数据依次输入，然后按一下键就可得到上述结果。 b与a的期望与方差 在介绍最小二乘法时我们曾提到，不管实际上X与Y之间有没有线性关系，用这种方法总是可以得到解的。因此我们必须有一种方法可以检验得到的结果是不是反映了X和Y之间的真实关系。为此，我们需要研究b与a的期望与方差。  注意  ∴ 原式=  ∵ 各yi互相独立，且D(yi)=σ2；各xi为常数； ∴     为估计σ2，令：，称为残差或剩余。则残差平方和为：   由于  （∵交叉项期望为0）  且 D(Sxy) = Sxx(2, E(Sxy) = (Sxx , （已证）   用MSe（剩余均方）代替(2，可得b与a的样本方差：   由于MSe的自由度为n-2，因此上述两方差的自由度也均为n-2。有了a和b的方差与均值，我们就可构造统计量对它们进行检验： H0 : ( = 0 HA: ( ( 0 （双侧检验） 或： HA: ( > 0 （或(< 0） （单侧检验） 统计量： （5.8） 当H0成立时，tb ~ t(n-2)，可查相应分位数表进行检验。 H0:( = 0 HA:( ( 0 （双侧检验） 或： HA:( > 0 （或( < 0） （单侧检验） 统计量： （5.9） 当H0成立时，ta ~ t(n-2)，可查相应分位数表进行检验。 在对一个回归方程的统计检验中，我们更关心的是(是否为0，而不是(是否为0。这是因为若( = 0，则线性模型变为Y = ( + (，与X无关；这意味着X与Y间根本没有线性关系。反之，(是否为0并不影响X与Y的线性关系。因此我们常常只对(作统计检验。 例5.2 对例5.1中的(作检验：H0: (=0 解：    查表，t0.995(3) = 5.841 < t，( 差异极显著，应拒绝H0，即( ( 0，或X与Y有着极显著的线性关系。 上述统计量还有一个用途：进行两个回归方程间的比较。即检验H0: (1 = (2和H0: (1 = (2。如果两H0均被接受，则可认为两组数据是抽自同一总体，从而可将两回归方程合并，得到一个更精确的方程。我们通过例5.3对具体方法加以说明。 例5.3 两组实验数据如下： x1 91 93 94 96 98 102 105 108  y1 66 68 69 71 73 78 82 85   x2 80 82 85 87 89 91 95  y2 55 57 60 62 64 67 71  是否可从它们得到统一的回归方程？ 解：从原始数据计算可得： 组别 n   Sxx Syy Sxy MSe b a  1 8 98.375 74.0 257.875 336.0 294.0 0.1357 1.140 -38.15  2 7 87.0 62.286 162.0 187.429 174.0 0.1080 1.074 -31.15  (1). 首先检验总体方差是否相等：   查表，F0.975(6, 5) = 6.978 > F, (接受H0，可认为两总体方差相等。 计算公共的总体方差：  (2). 检验回归系数(1与(2是否相等：H0: (1 = (2； HA: (1 ( (2  查表，得t0.975(11) = 2.201 > t, (接受H0，可认为两回归系数相等。 共同总体回归系数的估计值为：  (3). 再检验(1，(2是否相等：H0: (1 = (2； HA: (1 ( (2  查表，t0.975(11) = 2.201,  接受H0，可认为: (1 = (2。 若检验结果为(1 ( (2，此题即可结束；但若检验结果为(1 = (2，则需把全部原始数据放在一起，重新进行回归： Sxx = 902.9333, Sxy = 965.4667, Syy = 1035.7333, = 93.067, = 68.533, b == 1.0693, a == ?30.9787 从而得到合并的回归方程。 一元回归的方差分析 对回归方程的统计检验除可用上述t检验外，还有一些其他方法。这里我们再介绍一种方差分析的方法，它的基本思想仍是对平方和的分解。 无重复的情况。 y的总校正平方和可进行如下的分解：    即： Syy = SSe + SSR y的总校正平方和 残差平方和 回归平方和 自由度： n-1 n-2 1 这样就把y的总校正平方和分解成了残差平方和与回归平方和。前已证明，MSe可作为总体方差(2的估计量，而MSR可作为回归效果好坏的评价。如果MSR仅由随机误差造成的话，说明回归失败，X和Y没有线性关系；否则它应显著偏大。因此可用统计量  （5.10） 对H0: ( = 0进行检验。若F < F((1, n-2)，则接受H0，否则拒绝。 现在我们来证明这里的F检验与前述的t检验是一致的： 前已证明：SSe = Syy ? b ( Sxy, (SSR = Syy ? SSe = b ( Sxy,   例5.4 对例5.1作方差分析 解：由以前计算结果： Syy = 210.2，df = 4; SSe = 3.1704, df = 3, ( SSR = 210.2 ?3.1704 = 207.03, df = 1  查表得F0.95(1, 3) = 10.13, F0.99(1, 3) = 34.12 F > F0.99(1, 3)，拒绝H0，差异极显著。即应认为回归方程有效。 有重复的情况： 设在每一个xi取值上对Y作了m次观察，结果记为yi1, yi2, ……yim, 则线性统计模型变为： , i = 1, 2, … n, j = 1, 2, … m 估计值仍为： 现在y的总校正平方和可分解为： Syy = SSR + SSLOF + SSpe 其中SSLOF称为失拟平方和，SSpe为纯误差平方和，它们的表达式和自由度分别为：    同学们可试证明上述分解中的三个交叉项均为0。 统计检验步骤为： I. 令，它服从F(n-2, mn-n) （5.11） 若F检验差异显著，则可能的原因有： （1）除X以外还有其他变量影响Y的取值，而统计时没有加以考虑； （2）模型不当，即X与Y之间不是线性关系； 此时无必要再进一步对MSR作检验，而应想办法找出原因，并把它消除后重作回归。 若差异不显著，则把MSLOF和MSpe合并，再对MSR作检验： II. ，它服从F(1, mn-2) （5.12） 若差异显著，说明回归是成功的，X, Y间确有线性关系；若差异仍不显著，则回归失败，其可能的原因为： （1）X，Y无线性关系； （2）误差过大，掩盖了X, Y间的线性关系。 如有必要，可设法减小实验误差，或增加重复数重做实验后再重新回归。 点估计与区间估计 前边已经证明a和b是α和β的点估计，a+bx是y的点估计；但作为预测值仅给出点估计是不够的，一般要求给出区间估计，即给出置信区间。本节的重点就是讨论α，β，及y的置信区间。 α和β的区间估计 我们已经证明a和b是α和β的点估计，并求出了它们的方差。因此给出置信区间就很容易了：  ∴β的95%置信区间为：  （5.13） 同理  ((的95%置信区间为：  (5.14) 这与以前假设检验中的置信区间求法完全一样。若置信水平为99%，把分位数相应换为t0.995(n-2)即可。 例5.5 对例5.1中的(和(给出95%置信区间。 解：从前边的计算可知： a = 2.6996, b = 1.5167, Sxx = 90, MSe = 1.0568, n = 5,  查表，得t0.975(3) = 3.182  ( (的95%置信区间为： 2.6996 ( 4.3887, 即（-1.6891, 7.0883） (的95%置信区间为： 1.5167 ( 0.3448, 即（1.1719, 1.8615） 2. 对条件均值(Y? X的估计。 的点估计： 证明： 区间估计：首先需求出的方差。  (  用MSe代替(2，可得的1 ? (置信区间为：  （5.15） 注意上述置信区间的宽度与有关，当时，其宽度最小，偏离后，逐渐加大。 3. 对一次观察值y0的估计 y0的点估计： 证明： 区间估计： 一般情况下置信区间是以随机变量的期望为中点，此时只要求方差就可以了，因为方差就是衡量随机变量以数学期望为中心的离散程度的统计量。而现在是以条件均值的估计值，即另一个随机变量为中点，因此应求这两个随机变量差值的方差。由于下一次观察值y0和以前所有的观察值yi都是互相独立的，而估计值是从以前的观察值yi计算出来的，因此与y0独立，从而有：  由于y0和均为正态分布，它们的差也为正态分布。用代替后，为t分布，即：  ( 在x = x0处y0的1-(置信区间为：  (5.16) 显然y0的置信区间宽度也与x0有关，时最小，偏离时增大。y0的置信区间比的大一点，这是因为y0自己也有一个随机误差(。 例5.6 江苏武进县测定1959-1964年间3月下旬至4月中旬平均温度累积值x和一代三化螟蛾盛发期y的关系如下表(盛发期以5月10日为起算日)：试作回归分析。 表5.2 平均温度累积值与一代三化螟盛发期 年代 1956 1957 1958 1959 1960 1961 1962 1963 1964  累积温x 35.5 34.1 31.7 40.3 36.8 40.2 31.7 39.2 44.2  盛发期y 12 16 9 2 7 3 13 9 –1  解：由原始数据算得：  Sxx = 144.6356, Syy = 249.5556, Sxy= –159.0444, ( b ≈ –1.0996,   SSR = bSxy = 174.8886  查表，得：F0.95(1, 7) = 5.591, F0.99(1, 7) = 12.25, F > F0.99(1, 7), ( 拒绝H0，差异极显著。即X，Y有极显著线性关系。 为把上述回归结果用于预报，可给出观察值y0的95%置信区间：  查表，得t0.975(7) = 2.365, 把数据代入上式，得：  条件均值的95%置信区间公式为：  代入数据，得：  把不同的x0取值代入上述公式，可得置信区间的数据及图形如下： 表5.3 一代三化螟盛发期置信区间 x0  y0 的95%置信区间 y0的95%置信区间    下限 上限  下限 上限  30 15.6 10.3 20.8 6.2 24.9  32 13.4 9.2 17.5 4.6 22.1  34 11.2 7.9 14.4 2.8 19.5  36 9.0 6.3 11.6 0.8 17.1  38 6.8 4.1 9.4 -1.4 14.9  40 4.6 1.4 7.8 -3.8 12.9  42 2.4 -1.7 6.4 -6.4 11.1  44 0.2 -5.0 5.3 -9.1 9.4  46 -2.0 -8.3 4.2 -12.0 7.9   图5.2 一代三化螟盛发期置信区间 回归分析的目的常常是为了预报，也就是说下一次我们知道了x0的取值后，在观察前就对y0的取值作出估计。例如表5.3中的数据就是为了预报用的，下一年度如果我们知道了3月下旬至4月中旬的平均温度累积值，就可以估计出一代三化螟蛾盛发期是5月的什么时候。要特别注意的一点是预报范围只能是我们研究过的自变量变化范围，例如在上例中，当积温值是在32到44的范围内时，使用这一预报公式比较有把握，30和46使用已有点勉强，再大或小就不能用了。这是因为一般来说直线关系只是局部的近似，在更大的范围内，变量间常常呈现一种非线性的关系。因此若贸然把局部研究中发现的线性关系推广到更大的范围，常常是要犯严重错误的。同时从置信区间的宽度也可看出，即使是在研究的范围内，也是越接近所研究区间的中点（）预报越准确。 §5.2 相关分析 相关分析主要包括两方面的内容，即两个随机变量间回归方程的建立以及相关系数的概念及用途。前者由于X也是随机变量，所用的数学模型及分析方法与前一节中的方法相比都有所不同，而且需要更多的数学知识。但最终得到的a与b的公式则与前一节中的完全相同，因此我们不再作严格的数学推导，而是直接使用前一节中的公式。换句话说，使用那些公式时不必区分X是否是随机变量。但应该注意，a, b的公式推导都是建立在使Y的残差平方和最小这一原则上的。如果X，Y真是处在一种对等的地位，它们互为因果，那么我们就可以把Y当作自变量，X当作因变量，重新回归得到另一组回归方程的a’, b’。它们是建立在使X的残差平方和最小这一原则上的。一般来说，这两组a, b和a’, b’所代表的两条回归线在X-Y平面中不可能重合。这样一来就产生了一个问题：如果我们的目的就是研究两个处于对称位置的随机变量之间的关系，那究竟应选取哪一条回归线呢？一种简单而直观的解决办法就是认为真正代表X-Y关系的直线是通过这两条回归线交点并平分它们夹角的那一条。其他方法还有专为解决此问题而发展的主轴法，约化主轴法等，由于使用不多，我们不再详细介绍。如需要使用请参考陶澍先生编著的《应用数理统计方法》，354~360页。该书由中国环境科学出版社于1994年出版。本节的内容将主要集中在介绍相关系数的概念及它的应用上。 相关系数。 在第二章介绍随机变量的数字特征时，曾介绍了协方差及相关系数的概念： COV(X, Y) = E[(X-E(X))(Y-E(Y))] (5.17)  (5.18) 式中分别为X，Y的标准差。 对X，Y的一组观察值（xi, yi），i = 1, 2, ……n，我们可以有相应的样本协方差和样本相关系数：样本协方差定义为交叉乘积和除以它的自由度，即Sxy/(n-1)，然后用样本协方差和样本方差代替（5.18）式中的总体协方差和总体方差，可得样本相关系数：  (5.19) r的下标xy是指明为x和y的相关系数，在不会引起混淆的情况下常常把它省略。类似地，由于最常用的是样本相关系数而不是总体相关系数，因此常把样本相关系数简称为相关系数。 严格地说，只有当X，Y均为随机变量时才能定义相关系数，这一点在（5.17）, (5.18)式中可看得很清楚。这样一来，在本章的大多数情况下，由于我们假设X为非随机变量，相关系数根本就无法定义。但一方面不管X是不是随机变量，根据（5.19）式样本相关系数总是可以计算的；另一方面后边关于对样本相关系数进行统计检验的推导中，也并没有受到X必须为随机变量的限制，因此在回归分析中我们就借用了相关系数的名称和公式，而不再去区分X是否为随机变量。这一点在使用中是很方便的。 根据以前的推导结果，有：  因此, 。 当时，从上式可看出SSe = 0，即用可以准确预测y值。此时若X不是随机变量，则Y也不是随机变量了。这种情况在生物学研究中是不多见的。 当r = 0时，SSe = Syy，回归一点作用也没有，即用X的线性函数完全不能预测Y的变化。但这时X与Y间还可能存在着非线性的关系。 当时，情况介于上述二者之间隔。X的线性函数对预测Y的变化有一定作用，但不能准确预测，这说明Y还受其他一些因素，包括随机误差的影响。 综上所述，r可以作为X，Y间线性关系强弱的一种指标。它的优点是非常直观，接近于1就是线性关系强，接近于0就是线性关系弱；而其他统计量都需要查表后才知检验结果。 由于r是线性关系强弱的指标，我们当然希望能用它来进行统计检验。在一般情况下r不是正态分布，直接检验有困难。但当总体相关系数ρ= 0时，r的分布近似于正态分布，此时用MSe代替，就可以对作t检验。这种检验与对回归系数b的检验：是等价的。可证明如下： b的t检验统计量为：t = b/Sb。 b=Sxy/Sxx,  代入t的表达式，得： 。 因此我们可用上述统计量对作统计检验。 为使用方便，已根据上述公式编制专门的相关系数检验表，可根据剩余自由度及自变量个数直接查出r的临界值。 若必须对ρ≠0的情况作统计检验，可采用反双曲正切变换：  (5.20) 当n充分大时，可证明Z渐近正态分布N, 其中。利用统计量Z可对等进行检验。但这一检验方法用得很少。 例5.7 求出例5.1回归系数r，并作统计检验。 解：利用以前的计算结果，可得：  这里求得的Z值与例5.2中求得的t值是相同的，它们本来就是同一个统计量。 查表，t0.995(3) = 5.841 < t, ∴差异极显著，即X与Y有极显著的线性关系。 若直接查相关系数检验表，可得：剩余自由度为3，独立自变量为1，α=0.05的r临界值为0.878, α=0.01的临界值为0.959, ∴差异仍为极显著。 相关系数与回归系数间的关系 在X和Y均为随机变量的情况下，我们通常可以X为自变量，Y为因变量建立方程，也可反过来，以Y为自变量，X为因变量建立方程。此时它们的地位是对称的。 取X为自变量，Y为因变量，回归系b为： 取Y为自变量，X为因变量，回归系数b’为：  即：相关系数实际是两个回归系数的几何平均值。这正反映了相关与回归的不同：相关是双向的关系，而回归是单向的。 现在我们已介绍了三种对回归方程作统计检验的方法：对回归系数b作t检验，方差方析，对相关系数r作检验。对一元线性回归来说，它们的基本公式其实是等价的，因此结果也是一致的。但它们也各有自己的优缺点：对b的t检验可给出置信区间；方差分析在有重复的情况下可分解出纯误差平方和，从而可得到进一步的信息；相关系数则既直观，又方便（有专门表格可查），因此使用广泛。 最后要提请注意的一点是，不论采用什么检验方法，数据都应满足以下三个条件：独立，抽自正态总体，方差齐性。 §5.3 多元线性回归 前几节我们介绍了只有一个自变量时回归分析的方法。但在实际问题中，因变量常受不只一个自变量的影响。例如植物生长速度就可能受温度，光照，水分，营养等许多因素的影响。在这种情况下抛开其他因素不管只考虑一个因素显然是不适当的。因此我们有必要研究多个自变量的回归分析。 多元线性回归方程 k个自变量的情况下，线性回归模型变为：  （5.21） 其中，即它们为独立同分布的正态随机变量。 为求出各回归系数和, j=1, 2,……k的值，同样采用最小二乘法，即用使残差平方和  达到最小的和，j=1, 2, ……k作为和的估计值。其中  p = 1,2,……n。 令关于a和bj各的偏导数为0，可得：  整理，得正规方程如下：  由上述方程组中第一个方程解得：  （5.22） 代入其余方程，得：  （5.23） 其中  1≤i≤k, 1≤j≤k 从上述方程组中可解得b1,b2,……bk,从而求得a。可证明它们分别为β1，β2……βk和α的无偏估计量。bj称为Y对Xj的偏回归系数，它表示其他自变量固定时，Xj改变一单位所引起的Y的平均改变量。 从上述公式可见，多元回归的计算是相当麻烦的，现在通常用计算机完成。在确有多个因素影响因变量的情况下，应使用多元回归，否则会造成回归分析的失败。 矩阵解法 由于上述公式过于繁杂，应用及推导都很不方便。为简化这些表达式，可引入矩阵表示法。矩阵就是矩形的数表，一般用黑体字母表示。人们在它上面定义了一些特殊的运算规则，如加法、乘法、转置、求逆、微分等。它在数学中有许多应用，涉及多元问题时一般都要使用它。有关矩阵的初步知识可参见书后的附录。较详细的资料请参考“线性代数”教科书。 多元回归可用矩阵表示如下：令       其中β0 =α，b0 = a。 使用以上矩阵符号，线性回归模型可表示为： Y = Xβ+ε （5.24） 估计值为：  （5.25） 残差为：  （5.26） 残差平方和为： SSe = e(e = (Y - XB)((Y - XB) = (Y( - B(X()(Y - XB) = Y(Y - B(X(Y - Y(XB + B(X(XB = Y(Y - 2Y(XB + B(X(XB (5.27) （注意上式中每一项均为一个数字，而不是一个矩阵。） 对B求偏导，得：  （5.28） （根据矩阵微分法则，） 令(5.28)式等于0，得正规方程为： X(XB = X(Y （5.29） ( B = (X(X)-1X(Y （5.30） B的期望和方差为： E(B) = E[(X(X)-1 X(Y] = (X(X)-1 X(( E(Y) = (X(X)-1 X(( E(X( + () = (X(X)-1 X(((X( + E(()) = (X(X)-1 ( X(Xβ = ( 即：B为(的无偏估计。 D(B) = D[(X(X)-1 X(Y] = (X(X)-1X( ( D(Y) ( X(X(X)-1 = (X(X)-1 X( ( I ( (2 ( X(X(X)-1 = (2 (X(X)-1 （∵ Y的各分量独立，且方差均为(2） 上述矩阵主对角线上的元素是b0, b1, ……bk的方差，其他元素是各回归系数bj两两之间的协方差，因此可写为：  （5.31） 从上述推导过程可见，采用矩阵表示法后，多元回归的过程确实显得简单了不少。 多元回归的统计检验 回归方程的显著性检验 回归方程的显著性检验实际是检验所有的xj, j=1, 2, ……k作为一个整体与Y的线性关系是否显著。其假设为： H0: (1 = (2 = ……(k = 0, HA: 至少一个(j≠0 1≤j≤k 检验方法仍为方差分析。可以证明，在多元回归的情况下y的校正平方和仍可分解为回归平方和与残差平方和两部分：   它们的自由度分别为n-1, n-k-1, 和k。采用（5.23）式中的记号，可得：  其中  因此，我们可用统计量  （5.32） 作检验。当H0成立时，F～F（k,n-k-1）； H0不成立时SSR有增大的趋势，所以应使用上单尾检验。 若上述检验拒绝H0:β1 =β2 = …… =βk = 0，则应进一步对各βj，j=1,2,……k作t检验，以剔除不重要的因素。由于这里只需对各βj = 0作检验，因此可分别作t检验。 前已证明，(2（X(X）-1主对角线上的元素是各bj的方差。记C =（X(X）-1，则有：  用MSe代替总体参数(2, 得 ， 在H0:βj = 0下，统计量  （5.33） 也可采用对偏回归平方和作检验来代替上述t检验。偏回归平方和即取消一个自变量后所引起的回归平方和的减少量。即：  （5.34） 其中为去掉自变量xi后，用剩下的k-1个自变量作回归所得到的计算结果。SSPi称为Y对Xi的偏回归平方和。可以证明，SSPi = ,自由度为1。因此，可用统计量 F = SSPi / MSe （5.35） 作上单尾检验。当H0成立时，F～F(1, n-k-1)。由于  因此这一F检验与前述t检验等价。 若对某一(j的检验不显著，则接受H0: (j = 0,即说明相应的自变量xj对因变量Y没有明显影响，可将它从变量组中剔除。每剔除一个自变量后，都应对方程重新进行回归。 在剔除不重要的自变量时，应注意： 1( 每次只能剔除一个自变量。这是因为剔除掉一个自变量后，它对Y的影响很可能会转加到别的与它相关的自变量上，这样那些原先不重要的自变量也许会变得重要。 2( 由于前述原因，在一次检验中，偏回归平方和大到显著的一定应该保留；偏回归平方和最小的若不显著则可剔除，其他的不管显著与否都应待重作回归后再作检验。 复相关系数和偏相关系数 复相关系数定义为：  （5.36） 它实际上是y与的相关系数，或y与所有xj构成的整体的相关系数。对它的检验相当于对整个回归方程作方差分析。检验可通过查表进行。复相关系数与普通相关系数的不同点是它不取负值。 偏相关系数是保持其他变量不变的条件下计算的两个变量间的相关系数。它的计算公式为：  （5.37） 其中C =(X(X)-1, Cij为矩阵C的元素。对偏相关系数的检验也可通过查表进行。 偏相关系数和复相关系数查表时均使用MSe的自由度：n-k-1。对它们的检验与对回归平方和及偏回归平方和的检验是等价的。 逐步回归介绍 最优的回归方程应该是既没有包含多余的（即不显著的）自变量，也没有遗漏任何必要的（即显著的）自变量。要做到这一点可使用许多方法，而逐步回归是其中较好的一种。它的基本思想是采用偏回归平方和为检验标准，每次从未进入方程的自变量中选取偏回归平方和最大的一个进行检验。若显著，则引入回归方程，重作回归；再选已进入方程的自变量中偏回归平方和最小的一个进行检验，若不显著则剔除，并重作回归；……反复重复这一步骤，直到不能引入也不能剔除为止，这样就得到了最优的回归方程。 逐步回归的主要步骤为： 1°首先建立数据的样本相关矩阵R(0) 2°利用第n-1步的相关矩阵R(n-1)，求出未引入方程的各自变量的偏回归平方和。取其最大的作F检验，与给定的Fα作比较。若大于Fα则把对应的自变量引入回归方程，即对R(n-1)作变换，得R(n)，并建立Y与所有已引入的自变量的回归方程。 3°利用R(n)，计算所有已引入的自变量的偏回归平方和（刚引入的不必算）。选最小的作F检验。若小于给定的Fα，则把它剔除。方法仍是对R(n)作变换，得到R(n+1)，它给出了新的回归方程及其他一些信息。 4°重复步骤3°，直到没有自变量可以剔除为止。 5°重复步骤2°，直到没有自变量可以引入为止。 6°计算出最优回归方程，给出复相关系数。 关于逐步回归有几点说明如下： （1）从介绍中可看出，它的计算工作量是相当大的，不用计算机很难完成。但比起其它方法，逐步回归的计算量还是比较小的。 （2）Fα的值不象以前的检验是查表得到的，而是由使用者指定的。这是因为一方面运算过程中自由度一直在变化，因此得反复查表，会增加计算量；另一方面显著性水平α本来就是人为指定的，取值非常准确并无统计学上的意义，因此也是不必要的。一般来说可以试几个不同的Fα值，Fα越大，回归方程中包含的自变量个数越少。应以自变量个数多少为标准选一个你满意的。即在能包括主要有影响的自变量，不明显降低复相关系数的情况下，尽量选取少一些的自变量个数，一般不超过3～5个。当然自变量个数主要依赖于你的问题的复杂程度。有时也可对引入和剔除设置不同的Fα，但这样有时会形成一种循环：几个自变量走马灯一样引入又剔除，总也停不下来。此时应重新设置Fα。 （3）逐步回归是一种很有用的方法，它允许我们尽量多地收集数据，然后由计算机来选择。在问题的机理不十分清楚，无法确定哪些是真正有影响的因素时，这种方法的优越性是十分明显的。 （4）哪个自变量会进入方程与所选择的自变量变化范围有关。本来不能进入的，扩大一下范围，或换一个范围，就可能进入了。 （5）一般来说，逐步回归方法所允许考虑的自变量数应小于n-1，其中n为总的数据组数。否则正规方程系数矩阵的逆不存在，计算无法进行。 （6）由于在通常情况下我们都是利用现成的程序进行逐步回归，在本书中略去了具体的计算公式。如有需要的同学可参考其他有关多元回归的教科书。 §5.4 非线性回归 线性回归虽然比较简单，但应用非常广泛。这主要是因为如果我们缩小研究范围，则任意非线性关系最后都可以用线性关系来近似。但是范围缩得太小了使用上会很不方便，一来不能对变量间的关系有一个整体上的把握，二来在不同取值范围内还要换用不同的方程，因此在许多情况下考虑两变量间的非线性关系还是很有用的。 非线性回归可分为两种情况，即已知曲线（公式）类型和未知曲线（公式）类型。这两种情况需要用不同的方法来解决。一般来说，如果已知曲线类型，回归效果会比较有保证；同时在多数情况下我们对所研究的对象都有一定了解，可以根据理论或经验给出可能的曲线类型，因此常用的还是已知曲线类型的回归。 已知曲线类型的回归。 确定曲线类型的方法主要有： a)从专业知识判断。例如单细胞生物生长初期数量常按指数函数增长，但若考虑的生长时间相当长，后期其生长受到抑制，则会变为“S”形曲线。生态学上种群增长的情况也类似。此时常用逻辑斯蒂（Logistic）曲线进行拟合；反映药物剂量与死亡率之间关系的曲线也呈“S”形，但常用概率对数曲线描述；酶促反应动力学中的米氏方程是一种双曲线；植物叶层中的光强度分布常用指数函数描述；等等。这些公式或者来源于某种理论推导，或者是一种经验公式。 b) 如果没有足够的专业知识可判断变量间的关系是哪种类型，则可用直观的方法，即散点图的方法来判断。方法是把（x,y）数据对标在座标纸上，然后根据经验判断它们之间是什么类型。如果看来有几种类型可用，但不知哪种较好，也可多做几次回归，然后用后边介绍的方法对结果进行比较，选一种最好的。 确定曲线类型之后，回归的任务就变成确定曲线公式中的参数，因此也可称为曲线拟合。常用的回归或拟合方法有： 1. 线性化的方法。即先对数据进行适当变换，使其关系变为线性之后再按线性回归做。这种线性化的方法虽然常用，但它的缺点也是十分明显的。例如它只能保证使变换后数据的线性方程残差最小，而得到的非线性方程对原始数据没有任何最优性可谈。有时甚至会出现变换后的数据与线性回归方程吻合很好，而原始数据与非线性回归方程的差别大得不可接受的情况。因此采用线性化的方法进行曲线回归后必须用相关指数进行直观检验（见后边曲线回归的检验）。另外，也不是所有的非线性方程都能用数据变换的方法线性化。实际上，只有少数几种简单的非线性方程可用这种方法线性化，对绝大多数非线性方程来说都不行。 下面我们介绍几种生物统计中常用的变换方法。 （1）采用指数，对数，倒数等函数对自变量和因变量进行适当变换，使它们的关系变为线性。如： 指数函数： 令y( = lny, a( = lna, 可得：y( = a( + bx 幂函数：y = axb 令y( = lny, a( = lna, x( = lnx, 可得：y( = a( + bx( 对数函数：y = a + blnx 令x( = lnx, 可得：y = a + bx( 米氏方程： 令  可得:  但逻辑斯蒂方程: 是无法用变量代换的方法线性化的. (2)概率对数变换。主要用于毒理学研究中求半致死剂量。剂量与死亡率之间的关系一般呈如下曲线：  图5.3 剂量与死亡率关系 该曲线呈“S”形，但两端不对称。对于这种曲线可先把剂量取对数，使曲线对称化；然后对死亡率按标准正态分布作变换，即把死亡率作为累积概率值P，查正态分布表求出对应的单侧分位数up。它们的数学关系为 P(X<up) = p 其中X～N(0，1)。一般来说，up与剂量对数之间可呈现较好的线性关系。 综上所述，线性化方法的优缺点主要有： 优点：变量代换后可按线性回归做，简单方便。 缺点：1( 不是所有非线性方程都能用变量代换线性化； 2( 即使方程类型不对，变量代换与线性回归都可照样进行，但结果没有任何用处； 3( 线性回归效果好并不意味着变换前的非线性回归效果也好，因此必须用下面的 方法对所得的非线性方程进行检验。 4( 理论上所得回归方程是对线性化后数据最优，而不是对原始数据最优，因此影响 回归效果。 2. 曲线拟合。这种方法不需要对方程进行线性化，其基本思想是在所有参数所组成的高维空间中进行搜索，直到找到使目标函数（常为误差平方和或误差绝对值之和）达到极小值的点。具体算法有许多种，如Newton（牛顿）法，Marquardt（麦夸特）法，Powell（包维尔）法等。Newton法除了要给出曲线的公式外，还要给出一阶、二阶导数；Marguardt法也需要公式和一阶导数；而Powell法只需给出公式，不需要导数。这些方法都需要在计算机上实现。 由于这种曲线拟合的方法没有经过变量代换，而是直接使用原始数据，得到的参数至少是局部最优的，一般比用线性化方法得到的参数要好。如果采用不同的初值多拟合几次，更有可能得到接近最优的结果。在各种曲线回归的方法中，曲线拟合所得结果误差之小常是其他方法无法企及的。这种方法的缺点主要是计算量大，如果参数数量较多，甚至现代计算机解起来也有一定困难。另外，有时使用曲线拟合也会碰到迭代不收敛的问题，从而得不到参数的估计值。总起来看，随着计算机技术的发展，计算量大逐渐不成为重要限制条件，而回归误差小的优点则越来越被人们重视，因此曲线拟合的方法使用越来越多。 综上所述，曲线拟合方法的优缺点主要有： 优点：1( 不需变量代换，可保证所得参数至少局部最优，回归误差小于其他方法； 2( 常有现成软件可用； 缺点：1( 需要反复搜索，计算量大，必须用计算机； 2( 由于结果只是局部最优，一般需要试用多个初值；有时会出现不收敛的情况； 3( 参数数量多时，计算量迅速增加； 4( 有些拟合方法需要有目标函数的一、二阶导数。 未知曲线类型的回归：多项式回归。 以上我们介绍了已知曲线类型的情况。有些时候所研究的问题过于复杂，不可能进行理论上的推导；又没有前人的经验公式可利用，从散点图上也看不出明显的规律，此时就只能试用多项式回归的方法。最常用的方法为： 设自变量与因变量的关系为：  （5.38） 令 x1 = x, x2 = x2, x3 = x3, … xk = xk, 上式可化为： y = a + b1x1 + b2x2 + … + bkxk ， （5.39） 成为多元一次的线性方程，从而可用多元回归的方法求出各参数估计值。 这种方法的优点是不须对曲线类型有任何了解，如有必要，也可加上一些其他超越函数项，如对数、指数、三角函数等。它的缺点是其理论基础是任何曲线都可以在某一邻域中用多项式逼近，而这一邻域可能很小。另一个缺点是多项式的项数受到数据组数的限制：一般来说，当项数等于数据组数时，求回归系数就变成了解方程组，而不是一个优化过程。其结果是最后得到一条曲曲弯弯，但通过各数据点的曲线。此时误差为0（因为曲线通过了所有数据点），但该曲线没有任何用处，因为它根本不能反映自变量和因变量之间的关系。反之，如果项数太少又可能使回归误差变得很大，影响回归效果。因此必须保证数据组数和项数之间有适当的差距，即保证有一定剩余自由度。另外，这种方法得到的系数也很难有什么生物学意义，这也限制了它的应用。总之，这种方法现在用得较少。 优点：不需给出曲线类型； 缺点：得到的曲线只能用于一个小邻域；回归误差一般较大；参数没有生物学意义；增加了变量数，使方程复杂化。 曲线回归的检验 对曲线回归好坏的评价一定要根据变换前的原始数据作出，不能只对变换后的直线方程作检验。这种检验主要针对线性化的方法和多项式逼近的方法，曲线拟合一般就是用剩余平方和作为目标函数，方法本身就保证使它最小，因此就不用检验了。 多项式逼近方法和线性化方法在做回归前都要经过非线性变换，所得结果再变换回去后，各项方差都会有变化，所以线性回归效果好并不意味着原始数据间的非线性方程也效果好。由于非线性方程的方差等都有变化，原来的一些关系如变差的分解，均方期望等也都不再成立，建立在这些基础上的各种回归分析的统计检验方法，如t检验，F检验等也就都不能用了，只能对不同的结果作直观比较。比较的标准常用的有两个： 1. 剩余平方和：SS剩余= （5.40） 剩余平方和必须用变换前的原始数据计算。显然剩余平方和越小回归效果越好。但由于随机误差的影响，它不可能无限减小，又无法确定统计检验的阈值，因此它比较适合用于比较几种不同变换方法的优劣。 相关指数，它的定义为：  （5.41） 它的计算仍需用变换前的原始数据。由于不知道它的理论分布，因此不能用查表的方法统计检验。也不能使用SS剩余 = Syy - SSR的公式计算SS剩余，因为此时不能证明交叉项为0。与直接使用SS剩余相比，它能给人一个比较直观的印象，如越近1越好；如果接近0甚至变成负值，则说明变换公式使用不当，或X与Y没有关系。它变成负值说明用估计值预测的效果还不如用预测，这一般是由于使用了错误的变换公式。 作 业 5.1 某县儿童年龄与平均身高数据如下： 年龄X（岁） 4.5 5.5 6.5 7.5 8.5 9.5 10.5  身高Y（cm） 101.1 106.6 112.1 116.1 121.0 125.5 129.2  试作回归分析。 5.2 两个品系小麦的穗长与穗重如下： 品系 I 穗长(cm) 47 38 35 41 36 46 46 38 44 44   穗重(g) 1.9 1.5 1.1 1.4 1.2 1.8 1.7 1.3 1.7 1.8  品系 II 穗长(cm) 35 35 40 50 20 25 44 48 43 44   穗重(g) 1.2 1.4 1.6 2.0 0.6 0.7 1.7 1.9 1.6 1.8  作回归分析，并检验两回归方程能否合并。 5.3 植株生长周数与高度数据如下： 周数X 1 2 3 4 5 6 7  高度Y(cm) 5 13 16 23 33 38 40  作回归分析，并画出散点图，标出回归线，条件均值及下次观察值的95%置信区间。 5.4 根据Sb的公式考虑回归分析中，自变量取值范围大些好还是小些好？ 5.5 证明有重复一元线性回归中有： Syy + SSR + SSLOF + SSpe 5.6 棉田调查资料数据如下,其中x1为千株/亩, x2为铃数/株,y为皮棉产量(斤/亩): x1 6.21 6.29 6.38 6.50 6.52 6.55 6.61 6.77 6.82 6.96  x2 10.2 11.8 9.9 11.7 11.1 9.3 10.3 9.8 8.8 9.6  Y 190 221 190 214 219 189 183 199 182 201  试进行回归分析并检验。 5.7 烟草经X射线照射不同时间后的病斑数如下: 照射时间(分) 0 3 7.5 15 30 45 60  病斑数(个) 271 226 209 108 59 29 12  标出散点图，并进行曲线回归。(提示：曲线类型为y = aebx) 5.8 将蝽蟓暴露在-12.2℃低温下的时间x与死亡率关系如下： 时间x(小时) 0.25 0.5 1 4 12 24 48 72  死亡率y(%) 59 63 65 68 70 73 74 75  试进行曲线回归。 5.9 田间生长的大麦有如下生长曲线，计算它的回归方程（参考习题5.10中之（6）） 日数(天) 15 20 25 30 35 40 50 60 70 80 90 100 110 120  高度(cm) 4 5 6 7.5 8 10 15 20 30 48 60 65 67 69  5.10以下为几种常见曲线，试将它们直线化。 双曲线：  幂函数：  指数函数：  指数函数：  对数函数：  （6）S形曲线：