第5章一元线性回归 线性回归是统计分析中最常被用到的一种技术。在其他的领域,例如机器学习理论和计量经济研究中,回归分析也是不可或缺的重要组成部分。本章将要介绍的一元线性回归是最简单的一种回归分析方法,其中所讨论的诸多基本概念在后续更为复杂的回归分析中也将被常常用到。 5.1回归分析的性质 回归一词最早由英国科学家弗朗西斯·高尔顿(Francis Galton)提出,他还是著名生物学家、进化论奠基人查尔斯·达尔文(Charles Darwin)的表弟。高尔顿深受进化论思想的影响,并把该思想引入到人类研究,从遗传的角度解释个体差异形成的原因。高尔顿发现,虽然存在一个趋势——父母高,儿女也高; 父母矮,儿女也矮。但给定父母的身高,儿女辈的平均身高却趋向于或者“回归”到全体人口的平均身高。换句话说,即使父母双方都异常高或者异常矮,儿女的身高还是会趋向于人口总体的平均身高。这也就是所谓的普遍回归规律。高尔顿的这一结论被他的朋友,数学家、数理统计学的创立者卡尔·皮尔逊(Karl Pearson)所证实。皮尔逊收集了一些家庭的1000多名成员的身高记录,发现对于一个父亲高的群体,儿辈的平均身高低于他们父辈的身高; 而对于一个父亲矮的群体,儿辈的平均身高则高于其父辈的身高。这样就把高的和矮的儿辈一同“回归”到所有男子的平均身高,用高尔顿的话说,这是“回归到中等”。 回归分析是被用来研究一个被解释变量(Explained Variable)与一个或多个解释变量(Explanatory Variable)之间关系的统计技术。被解释变量有时也被称为因变量(Dependent Variable),与之相对应地,解释变量也被称为自变量(Independent Variable)。回归分析的意义在于通过重复抽样获得的解释变量的已知或设定值来估计或者预测被解释变量的总体均值。 图51父亲身高与儿子身高的关系 在高尔顿的普遍回归规律研究中,他的主要兴趣在于发现为什么人口的身高分布存在有一种稳定性。现在关心的是,在给定父辈身高的条件下,找出儿辈平均身高的变化规律。也就是一旦知道了父辈的身高,怎样预测儿辈的平均身高。图51展示了对应于设定的父亲身高,儿子在一个假想人口总体中的身高分布情况。不难发现,对于任一给定的父亲身高,我们都能从图中确定出儿子身高的一个分布范围,同时随着父亲身高的增加,儿子的平均身高也会增加。为了更加清晰地表示这种关系,在散点图上勾画了一条描述这些数据点分布规律的直线,用来表明被解释变量与解释变量之间关系,即儿子的平均身高与父亲身高之间的关系。这条直线就是所谓的回归线,后面还会对此进行详细讨论。 在回归分析中,变量之间的关系与物理学公式中所表现的那种确定性依赖关系不同。回归分析中因变量与自变量之间所呈现出来的是一种统计性依赖关系。在变量之间的统计性依赖关系中,主要研究的是随机变量,也就是有着概率分布的变量。但是函数或确定性依赖关系中所要处理的变量并非是随机的,而是一一对应的关系。例如,粮食产量对气温、降雨和施肥的依赖关系是统计性质的。这个性质的意义在于: 这些解释变量固然重要,但并不能据此准确地预测粮食的产量。首先是因为对这些变量的测量有误差,其次是还有很多影响收成的因素,很难一一列举。事实上,无论考虑多少个解释变量都不可能完全解释粮食产量这个因变量,毕竟粮食作物的生长过程是受到许许多多随机因素影响的。 与回归分析有密切关联的另外一种技术是相关分析,但两者在概念上仍然具有很大差别。相关分析是用来测度变量之间线性关联程度的一种分析方法。例如,常常会研究吸烟与肺癌发病率、金融发展与经济增长等之间的关联程度。而在回归分析中,对变量之间的这种关系并不感兴趣,回归分析更多的是通过解释变量的设定值来估计或预测因变量的平均值。 回归与相关在对变量进行分析时是存在很大分歧的。在回归分析中,对因变量和自变量的处理方法上存在着不对称性。此时,因变量被当作统计的、随机的,也就是存在着一个概率分布,而解释变量则被看成是(在重复抽样中)取有规定值的一个变量。因此在图51中,假定父亲的身高变量是在一定范围内分布的,而儿子的身高却反映在重复抽样后的一个由回归线给出的稳定值。但在相关分析中,将对称地对待任何变量,即因变量和自变量之间不加区别。例如,同样是分析父亲身高与儿子身高之间的相关性,那么这时我们所关注的将不再是由回归线给出的那个稳定值,儿子的身高变量也是在一定范围内分布的。大部分的相关性理论都建立在变量的随机性假设上,而回归理论往往假设解释变量是固定的或非随机的。 虽然回归分析研究是一个变量对另外一个或几个变量的依赖关系,但它并不意味着因果关系。莫里斯·肯达尔(Maurice Kendall)和艾伦·斯图亚蒂(Alan Stuart)曾经指出: “一个统计关系式,不管多强也不管多么有启发性,都永远不能确立因果关系的联系; 对因果关系的理念必须来自统计学以外,最终来自这种或那种理论。”比如前面谈到的粮食产量的例子中,将粮食产量作为降雨等因素的因变量没有任何统计上的理由,而是出于非统计上的原因。而且常识还告诉我们不能将这种关系倒转,即我们不可能通过改变粮食产量的做法来控制降雨。再比如,古人将月食归因于“天狗吃月”,所以每当发生月食时,人们就会敲锣打鼓意图吓走所谓的天狗。而且这种方法屡试不爽,只要人们敲锣打鼓一会儿,被吃掉的月亮就会恢复原样。显然,敲锣打鼓与月食结束之间有一种统计上的关系。但现代科技告诉我们月食仅仅是一种自然现象,它与敲锣打鼓之间并没有因果联系,事实上即使人们不敲锣打鼓,被“吃掉”的月亮也会恢复原状。总之,统计关系本身不可能意味着任何因果关系。要谈及因果关系必须进行先验的或理论上的思考。 5.2回归的基本概念 本节将从构建最简单的回归模型开始,结合具体例子向读者介绍与回归分析相关的一些基本概念。随着学习的深入,我们渐渐会意识到,更为一般的多变量之间的回归分析,在许多方面都是最简情形的逻辑推广。 5.2.1总体的回归函数 经济学中的需求法则认为,当影响需求的其他变量保持不变时,商品的价格和需求量之间呈反向变动的关系,即价格越低,需求量越多; 价格越高,需求量越少。据此,假设总体回归直线是线性的,便可以用下面的模型来描述需求法则 Ey|xi=w0+w1xi 这是直线的数学表达式,它给出了与具体的x值相对应的(或条件的)y的均值,即y的条件期望或条件均值。下标i代表第i个子总体,读作“在x取特定值xi时,y的期望值”。该式也称为非随机的总体回归方程。 这里需要指出,E(y|xi)是xi的函数,这意味着y依赖于x,也称为y对x的回归。回归可以简单地定义为在给定x值的条件下y值分布的均值,即总体回归直线经过y的条件期望值,而上式就是总体回归函数的数学形式。其中,w0和w1为参数,也称为回归系数。w0又称为截距,w1又称为斜率。斜率度量了x每变动一个单位,y的均值的变化率。 回归分析就是条件回归分析,即在给定自变量的条件下,分析因变量的行为。所以,通常可以省略“条件”二字,表达式E(y|xi)也简写成E(y)。 5.2.2随机干扰的意义 现通过一个例子来说明随机干扰项的意义。表51给出了21种车型燃油消耗(单位: L/100km)和车重(单位: kg)。下面在R中使用下列命令读入数据文件,并绘制散点图,还可以用一条回归线拟合这些散点。 > cars <- read.csv("c:/racv.csv") > plot(lp100km ~ mass.kg, data=cars, + xlab="Mass (kg)", ylab="Fuel consumption (l/100km)") > abline(lm(lp100km ~ mass.kg, data = cars)) 从表51中不难看出,车的油耗与车重呈正向关系,即车辆越重,油耗越高。如果用数学公式来表述这种关系,很自然地会想到采用直线方程来将这种依赖关系表示成下式 yi=Ey+ei=w0+w1xi+ui 其中,ui表示误差项。上式也称为随机总体回归方程。 表51车型及相关数据 MakeL/100kmmass/kg Alpha Romeo9.51242 Audi A38.81160 BA Falcon Futura12.91692 Chrysler PT Cruiser Classic9.81412 Commodore VY Acclaim12.31558 Falcon AU II Futura11.41545 Holden Barina7.31062 Hyundai Getz6.9980 Hyundai LaVita8.91248 Kia Rio7.31064 Mazda 27.91068 Mazda Premacy10.21308 Mini Cooper8.31050 Mitsubishi Magna Advance10.91491 Mitsubishi Verada AWD12.41643 Peugeot 3079.11219 Suzuki Liana8.31140 Toyota Avalon CSX10.81520 Toyota Camry Ateva V611.51505 Toyota Corolla Ascent7.91103 Toyota Corolla Conquest7.81081 易见,某一款车型的燃油消耗量等于两部分之和: 第一部分是由相应重量决定的燃油消耗期望Ey=w0+w1xi,也就是在重量取xi时,回归直线上相对应的点,这一部分称为系统的或者非随机的部分; 第二部分ui称为非系统的或随机的部分,在本例中由除了车重以外的其他因素所决定。 误差项ui是一个随机变量,因此,其取值无法先验地知晓,通常用概率分布来描述它。随机误差项可能代表了人类行为中一些内在的随机性。即使模型中已经包含了所有的决定燃油消耗的有关变量,燃油消耗的内在随机性也会发生变化,这是做任何努力都无法解释的。即使人类行为是理性的,也不可能是完全可以预测的。所以在回归方程中引入ui是希望可以反映人类行为中的这一部分内在随机性。 此外,随机误差项可以代表测量误差。在收集、处理统计数据时,由于仪器的精度、操作人员的读取或登记误差,总是会导致有些变量的观测值并不精准地等于实际值。所以误差项ui也代表了测量误差。 图52油耗与车重的关系 随机误差项也可能代表了模型中并未包括变量的影响。有时在建立统计模型时,并非事无巨细、无所不包的模型就是最好的模型。恰恰相反,有时只要能说明问题,建立的模型可能越简单越好。即使知道其他变量可能对因变量有影响,我们也倾向于将这些次要因素归入随机误差项ui中。 5.2.3样本的回归函数 如何求得总体回归函数中的参数w0和w1呢?显然在实际中,很难获知整个总体的全部数据。更多的时候,我们仅有来自总体的一部分样本。于是任务就变成了根据样本提供的信息来估计总体回归函数。下面来看一个类别数据的例子。 一名园艺师想研究某种树木的树龄与树高之间的关系,于是他随机选定了24株树龄在2~7年的树苗,每个树龄选择4棵,并记录每棵树苗的高度,具体数据如表52所示。表中同时给出了每个树龄对应的平均树高,例如对于树龄为2的4棵树苗,它们的平均树高是5.35。但在这个树龄下,并没有哪棵树苗的树高恰好等于5.35。那么我们如何解释在某一个树龄下,具体某一棵树苗的树高呢?不难看出每个树龄对应的一棵树苗的高度等于平均树高加上或减去某一个数量,用数学公式表达即为 yij=w0+w1xi+uij 某一个树龄i下,第j棵树苗的高度可以看作两个部分的和: 第一部分为该树龄下所有树苗的平均树高,即w0+w1xi,反映在图形上,就是在此树龄水平下,回归直线上相对应的点; 另一部分是随机项uij。 表52树高与树龄 树龄/年树高/m平均树高/m 25.64.85.35.75.350 36.25.96.46.16.150 46.26.76.46.76.500 57.17.36.96.97.050 67.27.57.87.87.575 78.99.28.58.78.825 在上述例子中,并无法获知所有树苗的高度数据,而仅仅是从每个树龄中抽取了4棵树苗作为样本。而且类别数据也可以向非类别数据转换,我们也会在后面演示R中处理这类问题的方法。 样本回归函数可以用数学公式表示为 y^i=w^0+w^1xi 其中,y^i是总体条件均值Ey|xi的估计量; w^0和w^1分别表示w0和w1的估计量。并不是所有样本数据都能准确地落在各自的样本回归线上,因此,与建立随机总体回归函数一样,我们需要建立随机的样本回归函数。即 yi=w^0+w^1xi+ei 式中,ei表示ui的估计量。通常把ei称为残差(Residual)。从概念上讲,它与ui类似,样本回归函数生成ei的原因与总体回归函数中生成ui的原因是相同的。 回归分析的主要目的是根据样本回归函数 yi=w^0+w^1xi+ei 来估计总体回归函数 yi=w0+w1xi+ui 样本回归函数是总体回归函数的近似。那么能否找到一种方法,使得这种近似尽可能地接近真实值?换言之,一般情况下很难获得整个总体的数据,那么如何建立样本回归函数,使得w^0和w^1尽可能接近w0和w1呢?我们将在下一小节介绍相关技术。 5.3回归模型的估计 本小节介绍一元线性回归模型的估计技术,并结合之前给出的树龄与树高关系的例子,演示在R中进行线性回归分析的方法。 5.3.1普通最小二乘法原理 在回归分析中,最小二乘法是求解样本回归函数时最常被用到的方法。本小节就来介绍它的基本原理。一元线性总体回归方程为 yi=w0+w1xi+ui 由于总体回归方程不能进行参数估计,因此只能对样本回归函数 yi=w^0+w^1xi+ei 进行估计。因此有 ei=yi-y^i=yi-w^0-w^1xi 从上式可以看出,残差ei是yi的真实值与估计值之差。估计总体回归函数的最优方法是,选择w0、w1的估计值w^0、w^1,使得残差ei尽可能小。最小二乘法的原则是选择合适的参数w^0、w^1,使得全部观察值的残差平方和为最小。 最小二乘法用数学公式可以表述为 min∑e2i=∑yi-y^i2=∑yi-w^0-w^1xi2 总而言之,最小二乘原理就是所选择的样本回归函数使得所有y的估计值与真实值差的平方和为最小。这种确定参数w^0和w^1的方法就叫做最小二乘法。 对于二次函数y=ax2+b来说,当a>0时,函数图形的开口朝上,所以必定存在极小值。根据这一性质,因为∑e2i是w^0和w^1的二次函数,并且是非负的,所以∑e2i的极小值总是存在的。根据微积分中的极值原理,当∑e2i取得极小值时,∑e2i对w^0和w^1的一阶偏导数为零,即 ∑e2iw^0=0,∑e2iw^1=0 由于 ∑e2i=∑yi-w^0-w^1xi2=∑yi-w^1xi2+w^20-2w^0yi-w^1xi 则得 ∑e2iw^0=-2∑yi-w^0-w^1xi=0 ∑e2iw^1=-2∑yi-w^0-w^1xixi=0 即 ∑yi=nw^0+w^1∑xi ∑xiyi=w^0∑xi+w^1∑x2i 以上两式构成了以w^0和w^1为未知数的方程组,通常叫做正规方程组,或简称正规方程。解正规方程,得到 w^0=∑x2i∑yi-∑xi∑xiyin∑x2i-∑xi2 w^1=n∑xiyi-∑xi∑yin∑x2i-∑xi2 等式左边的各项数值都可以由样本观察值计算得到。由此便可求出w0、w1的估计值w^0、w^1。 若设 x-=1n∑xi,y-=1n∑yi 则可以将w^0的表达式整理为 w^0=y--w^1x- 由此便得到了总体截距w0的估计值。其中,w^1的表达式如下 w^1=∑xiyi-nx-y-∑x2i-nx-2 这也就是总体斜率w1的估计值。 为了方便起见,在实际应用中,经常采用离差的形式表示w^0和w^1。为此设 x′i=xi-x-,y′i=yi-y- 因为 ∑x′iy′i=∑xi-x-yi-y-=∑xiyi-x-yi-xiy-+x-y- =∑xiyi-x-∑yi-y-∑xi+nx-y-=∑xiyi-nx-y- ∑x′2i=∑xi-x-2=∑x2i-2x-∑xi+nx-2=∑x2i-nx-2 所以w^0、w^1的表达式可以写成 w^0=y--w^1x-,w^1=∑x′iy′i∑x′2i 5.3.2一元线性回归的应用 上一小节中已经给出了最小二乘法的基本原理,而且还给出了计算斜率的几种不同方法。现在就以树高与树龄关系的数据为例来实际计算回归函数的估计结果。 正如前面说过的那样,类别数据可以转化成非类别数据,进而完成一元线性回归分析。其方法就是通过重复类别项从而将原来以二维数据表示的因变量转化为一维数据的形式。例如,在R中可以采用下列方法组织树高与树龄关系的数据。 > plants <- data.frame(age = rep(2:7, rep(4, 6)), + height = c(5.6, 4.8, 5.3, 5.7, 6.2, 5.9, 6.4, 6.1, + 6.2, 6.7, 6.4, 6.7, 7.1, 7.3, 6.9, 6.9, + 7.2, 7.5, 7.8, 7.8, 8.9, 9.2, 8.5, 8.7)) 上述代码将会得到如表53所示的数据组织形式。根据上一小节所得出的计算公式,我们还需计算相应的x2i和xiyij,这些数据也一并在表中列出。 表53树龄与树高数据 树龄/年 xi树高/m yijx2ixiyij树龄/年 xi树高/m yijx2ixiyij 25.6411.257.12535.5 24.849.657.32536.5 25.3410.656.92534.5 25.7411.456.92534.5 36.2918.667.23643.2 35.9917.767.53645.0 36.4919.267.83646.8 36.1918.367.83646.8 46.21624.878.94962.3 46.71626.879.24964.4 46.41625.678.54959.5 46.71626.878.74960.9 基于表53中的数据进而可以算得 x-=4.5,y-=6.908 nx-2=486,nx-y-=746.1 ∑x2i=556,∑xiyi=790.5 进而可以算得模型中估计的截距和斜率如下 w^1=∑xiyi-nx-y-/∑x2i-nx-2≈0.63429 w^0=y--w^1x-≈4.05405 由此便得到最终的估计模型为 y^i=4.05405+0.63429xi 或 yi=4.05405+0.63429xi+ei 当然,在R中并不需要这样繁杂的计算过程,仅需几条简单的命令就可以完成数据的线性回归分析。示例代码如下。 > plants.lm <- lm(height ~ age, data = plants) > summary(plants.lm) 由上述代码产生的模型估计如下,其中截距的估计值由Intercept项中的Estimate条目给出,斜率的估计值由age项中的Estimate条目给出,具体数值已经用方框标出。这些数据与我们人工算得的结果是一致的。输出结果中的其他数据将在后续的篇幅中加以讨论。 Call: lm(formula = height ~ age, data = plants) Residuals: Min1QMedian3QMax -0.65976 -0.22476 -0.008330.215240.70595 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept)4.054050.1937820.92 5.19e-16 *** age0.634290.0402615.76 1.82e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3368 on 22 degrees of freedom Multiple R-squared: 0.9186, Adjusted R-squared: 0.9149 F-statistic: 248.2 on 1 and 22 DF, p-value: 1.821e-13 模型的拟合结果由图53给出,代码如下。 > plot(height ~ age, data = plants) > abline(plants.lm) 图53线性回归拟合结果 5.3.3经典模型的基本假定 为了对回归估计进行有效的解释,就必须对随机干扰项ui和解释变量Xi进行科学的假定,这些假定称为线性回归模型的基本假定。主要包括以下几个方面。 1. 零均值假定 由于随机扰动因素的存在,yi将在其期望附近上下波动,如果模型设定正确,yi相对于其期望的正偏差和负偏差都会有,因此随机项ui可正可负,而且发生的概率大致相同。平均来看,这些随机扰动项有相互抵消的趋势。 2. 同方差假定 对于每个xi,随机干扰项ui的方差等于一个常数σ2,即解释变量取不同值时,ui相对于各自均值的分散程度是相同的。同时也不难推证因变量yi与ui具有相同的方差。因此,该假定表明,因变量yi可能取值的分散程度也是相同的。 前两个假设可以用公式ui~N0,σ2来表述,通常我们都认为随机扰动(噪声)符合一个均值为0,方差为σ2的正态分布。 3. 相互独立性 随机扰动项彼此之间都是相互独立的。如果干扰的因素是全随机的,相互独立的,那么变量yi的序列值之间也是互不相关的。 4. 因变量与自变量之间满足线性关系 这是建立线性回归模型所必需的。如果因变量与自变量之间的关系是杂乱无章、全无规律可言的,那么谈论建立线性回归模型就显然是毫无意义的。 R中提供了4种基本的统计图形,用于对线性回归模型的假设基础进行检验。下面就用车重与燃油消耗的例子来说明这几种图形的意义。在R中输入下列代码,则可绘制出如图54所示的4张统计图形。 > cars.lm <- lm(lp100km ~ mass.kg, data = cars) > par(mfrow = c(2, 2)) > plot(cars.lm) 图54(a)是一幅残差对拟合值的散点图。图中的x轴是拟合值,也就是当i取不同值时,有相应的y^i值。y轴表示的是残差值,即ei值。该图用于检验回归模型是否合理,是否有异方差性以及是否存在异常值。其中实线表示的附加线是采用局部加权回归散点修匀法(LOcally WEighted Scatterplot Smoothing,LOWESS)绘制的。如果残差的分布大致围绕着x轴,或红色附加线基本贴近x轴,则模型基本是无偏的; 另外,如果残差的分布范围不随预测值的改变而大幅变化,则可以认为同方差假设成立。所以图形显示其模型基本上没有什么问题。 图54线性回归模型的诊断信息 图54(b)展示了一幅标准化残差的QQ图,即将每个残差都除以残差标准差,然后再将结果与正态分布做比较。本书前面也已经对QQ图进行过较为详细的介绍,理想的结果是QQ图中的散点排列成一条直线,当然适度的偏离也是可以接受的。毕竟我们的采样点有限,根据中央极限定理,我们可以认为当采样点的数量足够大时,其结果会更加逼近正态分布。注意到在应用线性回归分析时,随机干扰项ui应当满足正态分布这个假定,而残差相当于是对ui的估计。如果图中散点的分布较大地偏离了直线,表明残差的分布是非正态的或者不满足同方差性,那么随机干扰的正态性自然也是不满足的。在我们给出的例子中,残差的正态性得到了较好的满足。 图54(c)作用大致与第一幅图相同。图中的x轴是拟合值,y轴表示的是相应的标准化残差值绝对值的平方根。如果标准化残差的平方根大于1.5,则说明该样本点位于95%置信区间之外。中间的实线偏离于水平直线的程度较大,则意味着异方差性。尽管图中的实线表示的附加线不是一条完全水平的直线,但这种小的偏离主要是因为样本点的数量较小,所以图形显示我们的模型基本上没有什么问题。 图54(d)是标准化残差对杠杆值的散点图,其作用是检查样本点中是否有异常值。如果删除样本点中的某一条数据,由此造成的回归系数变化过大,就表明这条数据对回归系数的计算产生了明显的影响,这条数据就是异常值。需要好好考虑是否在模型中使用这条数据。设有帽子矩阵H,该矩阵的诸对角线元素记为hii,这就是杠杆值(Leverage)。杠杆值用于评估第i个观测值离其余n-1个观测值的距离有多远。对一元回归来说,其杠杆值为 hii=1n+xi-x-2∑xi-x-2 此外,图中还添加了LOWESS曲线和库克距离(Cooks Distance)曲线。库克距离用于诊断各种回归分析中是否存在异常数据。库克距离太大的样本点可能是模型的强影响点或异常值点,值得进一步检验。一个通常的判断准则是当库克距离大于1时就需要引起注意,图中显示所有点的库克距离都在0.5以内,所以没有异常点。 在本小节最后,尝试在R中自行绘制图54(d)。这个过程有助于读者更好地理解杠杆值的意义。表54给出了操作步骤计算所得的中间结果。这些计算步骤需要用到的三个值,即斜率0.008024、截距-0.817768和的残差标准差0.3891,这些值都可以从线性回归的输出结果中直接得到。 表54中间结果数据 xi-x-2杠杆值y^iei标准化残差 2308.5740.0499039.1480400.3519600.904549 16912.380.0643478.4900720.3099280.796525 161565.70.20742712.758840.1411600.362786 14872.380.06232310.51212-0.712120-1.830170 71798.480.11863611.683620.6163761.584107 65000.720.11191311.57931-0.179310-0.460840 52005.720.0990597.703720-0.403720-1.037570 96129.530.1427037.045752-0.145750-0.374590 1768.0020.0493689.196184-0.296180-0.761200 51097.530.0981617.719768-0.419770-1.078820 49305.150.0963887.7518640.1481360.380714 322.28800.0479389.6776240.5223761.342524 57622.860.1046157.6074320.6925681.779923 40381.860.08756211.14601-0.246020-0.632270 124575.40.17083912.365660.0343360.088245 5047.7640.0526128.9634880.1365120.350840 22514.290.0698888.329592-0.029590-0.076050 52878.100.09992211.37871-0.578710-1.487310 46204.530.09332111.258350.2416480.621043 34986.810.0822258.032704-0.132700-0.341050 43700.910.0908447.856176-0.056176-0.144370 下面给出绘制图形的R代码。 > plot(Std_Residuals ~ Leverage, xlab="Leverage", + ylab="Standardized residuals", + xlim=c(0,0.21), ylim=c(-2,2), main = "Residuals vs Leverage") > abline(v = 0.0, h = 0.0, lty=3, col = "gray60") > par(new=TRUE) > lines(lowess(Std_Residuals~Leverage ), col = 'red') 执行上述代码,结果如图55所示,易见与R自动生成的效果一致。 图55标准化残差对杠杆值的散点图 5.3.4总体方差的无偏估计 前面谈到回归模型的基本假定中有这样一条: 随机扰动(噪声)符合一个均值为0,方差为σ2的正态分布,即ui~N0,σ2来表述。随机扰动ui的方差σ2又称为总体方差。由于总体方差σ2未知,而且随机扰动项ui也不可度量,所以只能从ui的估计量——残差ei出发,对总体方差σ2进行估计。可以证明总体方差σ2的无偏估计量为 σ^2=∑e2in-2=∑yi-y^i2n-2 证明: 因为 y-=1n∑yi 即y-是有限个yi的线性组合,所以当yi=w0+w1xi+ui,同样有 y-=w0+w1x-+u- 所以可得 y′i=yi-y-=w0+w1xi+ui-w0+w1x-+u- =w1xi-x-+ui-u-=w1xi+ui-u- 又因为 ei=yi-y^i=yi-w^0-w^1xi=y′i+y--w^0-w^1x′i+x- w^0=y--w^1x-ei=y′i-w^1x′i 所以有 ei=w1x′i+ui-u--w^1x′i=ui-u--β^1-β1x′i 进而有 ∑e2i=∑ui-u--w^1-w1x′i2 =w^1-w12∑x′2i+∑ui-u-2-2w^1-w1∑x′iui-u- 对上式两边同时取期望,则有 E∑e2i=Ew^1-w12∑x′2i+E∑ui-u-2 -2Ew^1-w1∑x′iui-u- 然后对上式右端各项分别进行整理,可得 E∑ui-u-2=E∑u2i-2uiu-+u-2=Enu-2+∑u2i-2u-∑ui =E∑u2i-1n∑ui2=∑Eu2i-1nE∑ui2 =∑Eu2i-1n∑u2i+2∑i≠juiuj =nσ2-1nnσ2-0=n-1σ2 其中用到了ui互不相关以及ui~N0,σ2这两条性质。 一个变量与其均值的离差之总和恒为零,该结论可以简证如下 x-=1n∑xinx-=∑xi∑x-=∑xi∑xi-x-=0 又因为y-是一个常数,所以有 ∑x′iy′i=∑x′iyi-y-=∑x′iyi-y-∑x′i =∑x′iyi-y-∑xi-x-=∑x′iyi 进而得到 w^1=∑x′iyi∑x′2i=∑x′iyi∑x′2i=∑kiyi 其中 ki=x′i∑x′2i 这其实说明w^1是y的一个线性函数; 它是yi的一个加权平均,以ki为权数,从而它是一个线性估计量。同理,w^0也是一个线性估计了。易证ki满足下列性质 ∑ki=∑x′i∑x′2i=1∑x′2i∑x′i=0 ∑k2i=∑x′i∑x′2i2=∑x′2i∑x′2i2=1∑x′2i ∑kix′i=∑kixi=1 于是有 w^1=∑kiyi=∑kiw0+w1xi+ui =w0∑ki+w1∑kixi+∑kiui=w1+∑kiui 即 w^1-w1=∑kiui 以此为基础可以继续前面的整理过程,其中再次用到了ui的互不相关性 Ew^1-w1∑x′iui-u-=E∑kiui∑x′iui-u- =E∑kiui∑x′iui-x′iu- =E∑kiui∑x′iui-u-∑kiui∑x′i =E∑kiui∑x′iui=E∑kix′iu2i=σ2 此外还有 Ew^1-w12∑x′2i=E∑kiui2∑x′2i =E∑x′i∑x′2iui2∑x′2i=E∑x′iui2/∑x′2i=σ2 综上可得 E∑e2i=n-1σ2+σ2-2σ2=n-2σ2 原结论得证,可知σ^2是σ2的无偏估计量。 5.3.5估计参数的概率分布 中央极限定理表明,对于独立同分布的随机变量,随着变量个数的无限增加,其和的分布近似服从正态分布。随机项ui代表了在回归模型中没有单列出来的其他所有影响因素。在众多的影响因素中,每种因素对yi的影响可能都很微弱,如果用ui来表示所有这些随机影响因素之和,则根据中央极限定理,就可以假定随机误差项服从正态分布,即ui~N(0,σ2)。 因为w^0和w^1是yi的线性函数,所以w^0和w^1的分布取决于yi。而yi与随机干扰项ui具有相同类型的分布,所以为了讨论w^0和w^1的概率分布,就必须对ui的分布做出假定。这个假定十分重要,如果没有这一假定,w^0和w^1的概率分布就无法求出,再讨论两者的显著性检验也就无的放矢了。 根据随机项ui的正态分布假定可知,yi服从正态分布,根据正态分布变量的性质,即正态变量的线性函数仍服从正态分布,其概率密度函数由其均值和方差唯一决定。于是可得 w^0~Nw0,σ2∑x2in∑x′2i w^1~Nw1,σ2∑x′2i 并且w^0和w^1的标准差分布为 图56估计量的分布及其偏移 sew^0=σ2∑x2in∑x′2i sew^1=σ2∑x′2i 以w^1的分布为例,如图56所示,w^1是w1的无偏估计量,w^1的分布中心是w1。易见,标准差可以用来衡量估计值接近于其真实值的程度,进而判定估计量的可靠性。 此前,已经证明σ^2是σ2的无偏估计量,那么由此可知w^0和w^1的方差及标准差的估计量分别为 varw^0=σ^2∑x2in∑x′2i,sew^0=σ^∑x2in∑x′2i varw^1=σ^2∑x′2i,sew^1=σ^∑x′2i 例如,在车重与油耗的例子中,一元线性回归的分析结果如下。其中,截距的估计值β^0的标准差为0.506422,斜率的估计值β^1的标准差为0.000387,这两个值已经用方框标出。 > summary(cars.lm) Call: lm(formula = lp100km ~ mass.kg, data = cars) Residuals: Min1QMedian3QMax -0.71186 -0.24574 -0.029380.241930.69276 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.8177680.506422-1.6150.123 mass.kg 0.0080240.00038720.733 1.65e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3891 on 19 degrees of freedom Multiple R-squared: 0.9577, Adjusted R-squared: 0.9554 F-statistic: 429.9 on 1 and 19 DF, p-value: 1.653e-14 标准差可以被用来计算参数的置信区间。例如在本题中,w0的95%的置信区间为 -0.8178±c0.975(t19)×0.5064 =-0.8178±2: 093×0.5064 =(-1.878,0.242) 同理可以计算w1的95%的置信区间为 0.008024±c0.975(t19)×0.000387 =0.008024±2.093×0.000387 =(0.0072,0.0088) 其中,因为残差的自由度为21-2=19,所以数值2.093是自由度为19的t分布值。当然在R中可以通过如下代码来完成上述计算过程。 > confint(cars.lm) 2.5 %97.5 % (Intercept) -1.877722151 0.24218677 mass.kg0.007213806 0.00883382 5.4正态条件下的模型检验 以样本观察值为基础,用最小二乘法求得样本回归直线,从而对总体回归直线进行拟合。但是拟合的程度怎样,必须要进行一系列的统计检验,从而对模型的优劣做出合理的评估,本节就介绍与模型评估检验有关的内容。 5.4.1拟合优度的检验 由样本观察值xi,yi得出的样本回归直线为y^i=w^0+w^1xi,y的第i个观察值yi与样本平均值y-的离差称为yi的总离差,记为y′i=yi-y-,不难看出总离差可以分成两部分,即 y′i=yi-y^i+y^i-y- 其中一部分y^′i=y^i-y-是通过样本回归直线计算的拟合值与观察值的平均值之差,它是由回归直线(即解释变量)所解释的部分。另一部分ei=yi-y^i是观察值与回归值之差,即残差。残差是回归直线所不能解释的部分,它是由随机因素、被忽略掉的因素、观察误差等综合影响而产生的。各变量之间的关系如图57所示。 图57总离差分解 由回归直线所解释的部分y^′i=y^i-y-的绝对值越大,则残差的绝对值就越小,回归直线与样本点xi,yi的拟合就越好。 因为 yi-y-=yi-y^i+y^i-y- 如果用加总y的全部离差来表示显然是不行的,因为 ∑yi-y-=∑yi-∑y-=ny--ny-=0 所以考虑利用加总全部离差的平方和来反映总离差,即 ∑yi-y-2=∑yi-y^i+y^i-y-2 =∑yi-yi2+∑y^i-y-2+2∑yi-y^iy^i-y- 其中 ∑yi-y^iy^i-y-=0 这是因为 ∑yi-y^iy^i-y-=∑eiw^0+w^1xi-y- =w^0-y-∑ei+w^1∑xiei =w^0-y-∑ei+w^1∑xiyi-w^0-w^1xi 注意最小二乘法对于ei有零均值假定,所以对其求和结果仍为零。而上述式子中最后一项为零,则是由最小二乘法推导过程中极值存在条件(令偏导数等于零)所保证的。 于是可得 ∑yi-y-2=∑yi-y^i2+∑y^i-y-2 或者写成 ∑y′2i=∑e2i+∑y^′2i 其中, ∑y′2i=∑yi-y-2 称为总离差平方和(Total Sum of Squares),用SStotal表示 ∑e2i=∑yi-y^i2 称为残差平方和(Residual Sum of Squares),用SSresidual表示 ∑y^′2i=∑y^i-y-2 称为回归平方和(Regression Sum of Squares),用SSregression表示。 总离差平方和可以分解成残差平方和与回归平方和两部分。总离差分解公式还可以写成 SStotal=SSresidual+SSregression 这一公式也是方差分析ANOVA的原理基础,这一点在后续的章节中我们还会详细介绍。 在总离差平方和中,如果回归平方和比例越大,残差平方和所占比例就越小,表示回归直线与样本点拟合得越好; 反之,就表示拟合得不好。把回归平方和与总离差平方和之比定义为样本判定系数,记为 R2=SSregression/SStotal 判断系数R2是一个回归直线与样本观察值拟合优度的数量指标,R2越大则拟合优度就越好; 相反,R2越小,则拟合优度就越差。 注意R中指示判定系数的标签是“Multiple Rsquared”,例如,在前面给出的树高与树龄的例子中,R2=0.9186=91.86%,这表明模型的拟合程度较好。此外,R的输出中还给出了所谓的调整判定系数,调整判定系数是对R2的修正,指示标签为“Adjusted Rsquared”。例如,在树高与树龄的例子中调整判定系数大小为0.9149。 在具体解释调整判定系数的意义之前,还需先考查一下进行线性回归分析时,R中输出的另外一个值——残差标准误差(Residual Standard Error)。在树高与树龄的例子中,R给出的结果数值是0.3368。所谓残差的标准误差其实就是残差的标准差(Residual Standard Deviation)。前面已经证明过,在一元线性回归中,总体方差σ2的无偏估计量为 σ^2=∑e2in-2=∑yi-y^i2n-2 所以残差的标准差为 s=σ^=∑yi-y^i2n-2 如果将这一结论加以推广(即不仅限于一元线性回归),则有 s=σ^=SSresidualn-被估计之参数的数量 因为在一元线性回归中,被估计的参数只有β0和β1两个,所以此时被估计之参数的数量就是2。而在树高与树龄的例子中,研究单元的数量n=24,因此在R中的输出结果上有一句“on 22 degrees of freedom”。 调整判定系数的定义为 1-R2adj=s2/s2y 根据前面给出的公式可知 s2=SSresidualn-p 其中,p是模型中参数的数量。以及 s2y=SStotaln-1 一般认为调整判定系数会比判定系数更好地反映回归直线与样本点的拟合优度。那么其理据何在呢?注意残差ei是扰动项ui的估计值,因为ui的标准差σ无法计算,所以借助ei对其进行估计,而且也可以证明其无偏估计的表达式需要借助自由度来进行修正。另一方面,本书前面也曾经证明过当用样本来估计总体时,方差的无偏估计需要通过除以n-1来进行修正。所以采用上述公式来计算会得到更加准确的结果。 经过简单的代数变换,可得出R2adj的另外一种算式 R2adj=R2-p-1n-p1-R2 对于树高与树龄的例子有 R2adj=0.9186-2-124-21-0.9186≈0.9149 这与R中输出的结果相同。通常情况下,R2adj的值都会比R2的值略小,且两者的差异一般都不大。 5.4.2整体性假定检验 如果随机变量X服从均值为μ、方差为σ2的正态分布,即X~N(μ,σ2),则随机变量Z=(X-μ)/σ是标准正态分布,即Z~N(0,1)。统计理论表明,标准正态变量的平方服从自由度为1的χ2分布,用符号表示为 Z2~χ21 其中,χ2的下标表示自由度为1。与均值、方差是正态分布的参数一样,自由度是χ2分布的参数。在统计学中自由度有各种不同的含义,此处定义的自由度是平方和中独立观察值的个数。 总离差平方和SStotal的自由为n-1,因变量共有n个观察值,由于这n个观察值受∑y′i=∑yi-y-=0的约束,当n-1个观察值确定以后,最后一个观察值就不能自由取值了,因此SStotal的自由为n-1。 回归平方和SSregression的自由度是由自变量对因变量的影响决定的,因此它的自由度取决于解释变量的个数。在一元线性回归模型中,只有一个解释变量,所以SSregression的自由度为1。在多元回归模型中,如果解释变量的个数为k个,则其中SSregression的自由度为k。因为SSregression的自由度与SSresidual的自由度之和等于SStotal的自由度,所以SSresidual的自由度为n-2。 平方和除以相应的自由度称为均方差。因此SSregression的均方差为 ∑y^′2i1=∑yi-y-2=∑w^0+w^1xi-y-2 =∑w^0+w^1x-+x′i-y-2=∑y--w^1+w^1x-+x′i-y-2 =∑w^1x′i2=w^21∑x′2i 而且还有SSresidual的均方差为∑e2i/(n-2)。可以证明,在多元线性回归的条件下(即回归方程中有k个解释变量xi,i=1,2,…,k),有 ∑y^′2i~χ2k ∑e2i~χ2(n-k-1) 根据基本的统计学知识可知,如果Z1和Z2分别是自由度为k1和k2的分布变量,则其均方差之比服从自由度为k1和k2的F分布,即 F=Z1/k1Z2/k2~Fk1,k2 那么 F=∑y^′2i/k∑e2i/(n-k-1)~Fk,n-k-1 下面就利用F统计量对总体线性的显著性进行检验。首先,提出关于k个总体参数的假设 H0: w1=w2=…=wk=0 H0: wi不全为0,i=1,2,…,k 进而根据样本观察值计算并列出方差分析数据如表55所示。 表55方差分析表 方 差 来 源平方和自由度均方差 SSresidual∑y^′2ik∑y^′2i/k SSregression∑e2in-k-1∑e2i/(n-k-1) SStotal∑y′2i 然后在H0成立的前提下计算F统计量 F=∑y′2i/k∑e2i/(n-k-1) 对于给定的显著水平α,查询F分布表得到临界值Fα1,n-k-1,如果F>Fα1,n-k-1,则拒绝原假设,说明犯第一类错误的概率非常之小。也可以通过与这个F统计量对应的P值来判断,说明如果原假设成立,得到此F统计量的概率很小即为P值。这个结果说明我们的回归模型中的解释变量对因变量是有影响的,即回归总体是显著线性的。相反,若F<Fα1,n-k-1,则接受原假设,即回归总体不存在线性关系,或者说解释变量对因变量没有显著的影响关系。 例如,对于树龄与树高的例子,给定α=0.05,可以查表或者在R中输入下列语句得到F0.051,22的值。 > qf(0.05, 1, 22, lower.tail = FALSE) [1] 4.30095 其中,参数lower.tail是一个逻辑值,模型情况下它的值为FALSE,此时给定服从某分布的随机变量X,求得的概率是P[X≤x],如果要求P[X>x],要么用1-P[X≤x],要么就令lower.tail的值为TRUE。 经过简单计算易知∑y′2i=28.1626663,∑e2i=2.496047632。由此便可算得F=248.2238923。当然,R中给出的线性回归分析结果也包含了这个结果。因为F>F0.051,22,所以有理由拒绝原假设H0,即证明回归总体是显著线性的。也可以通过与这个F统计量对应的P值来判断,此时可以在R中使用下面的代码得到相应的P值。 > pf(248.2238923, 1, 22, lower.tail = FALSE) [1] 1.821097e-13 可见,P值远远小于0.05,因此有足够的把握拒绝原假设。 本小节所介绍的其实就是方差分析(ANOVA)的基本步骤。在本书的后续章节中,还将对方差分析做专门介绍。一元线性回归模型中对模型进行整体性检验只用后面介绍的t检验即可。但在多元线性回归模型中,F检验是检验统计假设的非常有用和有效的方法。 5.4.3单个参数的检验 前面介绍了利用R2来估计回归直线的拟合优度,但是R2却不能告诉我们估计的回归系数在统计上是否显著,即是否显著地不为零。实际上确实有些回归系数是显著的,而有些又是不显著的。下面就来介绍具体的判断方法。 本章前面曾经给出了w^0和w^1的概率分布,即 w^0~Nw0,σ2∑x2in∑x′2i w^1~Nw1,σ2∑x′2i 但在实际分析时,由于σ2未知,只能用无偏估计量σ^2来代替。此时,一元线性回归的最小二乘估计量w^0和w^1的标准正态变量服从自由度为n-2的t分布,即 t=w^0-w0sew^0~tn-2 t=w^1-w1sew^1~tn-2 下面以w1为例,演示利用t统计量对单个参数进行检验的具体步骤。首先对回归结果提出如下假设 H0: w1=0 H1: w1≠0 即在原假设条件下,解释变量对因变量没有影响。在备择假设条件下,解释变量对因变量有(正的或者负的)影响,因此备择假设是双边假设。 以原假设H0构造t统计量并由样本观察值计算其结果,则 t=w1sew^1 其中 sew^1=σ^∑x′2i=∑e2in-2∑x′2i 可以通过给定的显著性水平α,检验自由度为n-2的t分布表,得临界值tα2n-2。如果t>tα2n-2,则拒绝H0,此时接受备择假设犯错的概率很小,即说明w1所对应的变量x对y有影响。 相反,若t≤tα2n-2,则无法拒绝H0,即w1与0的差异不显著,说明w1所对应的变量x对y没有影响,变量之间的线性关系不显著。对参数的显著性检验,还可以通过P值来判断,如果相应的P值很小,则可以拒绝原假设,即参数显著不为零。 例如,在树龄与树高的例子中,很容易算得 ∑x′2i=70 于是可得到sew^1=0.3368/70=0.04026,进而有t=0.63429/0.04026=15.75484。相应的P值可以在R中用下列代码算得。 > 2*(1-pt(15.75484,22)) [1] 1.820766e-13 经过计算所得之t值为15.75484,其P值几乎为0。P值越低,拒绝原假设的理由就越充分。现在来看,我们已经有足够的把握拒绝原假设,可见变量之间具有显著的线性关系。 5.5一元线性回归模型预测 预测是回归分析的一个重要应用。这种所谓的预测通常包含两个方面,对于给定的点,一方面要估计它的取值,另一方面还应对可能取值的波动范围进行预测。 5.5.1点预测 对于给定的x=x0,利用样本回归方程可以求出相应的样本拟合值y^0,以此作为因变量个别值y0或其均值Ey0的估计值,这就是所谓的点预测。比如树龄与树高的例子,如果购买了一棵树苗,并且想知道该树的树龄达到4年时,其树高预计为多少。此时你希望求得的值,其实是树龄为4的该种树木的平均树高或者是期望树高。 已知含随机扰动项的总体回归方程为 yi=Eyi+ui=w0+w1xi+ui 当x=x0时,y的个别值为 y0=w0+w1x0+u0 其总体均值为 Ey0=w0+w1x0 样本回归方程在x=x0时的拟合值为 y^0=w^0+w^1x0 对上式两边取期望,得 Ey^0=Ew^0+w^1x0=w0+w1x0=Ey0 这表示在x=x0时,由样本回归方程计算的y^0是个别值y0和总体均值Ey0的无偏估计,所以y^0可以作为y0和Ey0的预测值。 5.5.2区间预测 对于任一给定样本,估计值y^0只能作为y0和Ey0的无偏估计量,不一定能够恰好等于y0和Ey0。也就是说,两者之间存在误差,这个误差就是预测误差。由这个误差开始,期望得到y0和Ey0的可能取值范围,这就是区间预测。 定义误差δ0=y^0-Ey,由于y^0服从正态分布,所以δ0是服从正态分布的随机变量。而且可以得到δ0的数学期望与方差如下 Eδ0=Ey^0-Ey=0 varδ0=Ey^0-Ey2=Ew^0+w^1x0-w0+w1x02 =Ew^0-w02+2w^0-w0w^1-w1+w^1-w12x20 =varw^0+2x0covw^0,w^1+varw^1x20 其中,w^0和w^1的协方差为 covw^0,w^1=Ew^0-w0w^1-w1 =Ey--w^1x--w0w^1-w1 =Ew0+w1x-+u--w^1x--w0w^1-w1 =E-w^1-w1x-+u-w^1-w1 =x-Ew^1-w12+Eu-w^1 因为 Ew^1-w12=varw^1=σ2∑x′2i Eu-w^1=1nE∑ui∑x′i∑x′2iyi =1n∑i=jx′i∑x′2iEuiyi+1n∑i≠jx′i∑x′2iEuiyi =σ2∑x′i∑x′2iEuiyi=0 所以 covw^0,w^1=-x-σ2∑x′2i 于是可得 varδ0=σ2∑x′2in∑x′2i-2σ2x0x-∑x′2i+σ2x20∑x′2i =σ2∑x′2i∑x′2i-nx-n+x-2-2x0x-+x20 =σ2∑x′2i∑x′2in+x0-x-2=σ21n+x0-x-2∑x′2i 由δ0的数学期望与方差可知 δ0~N0,σ21n+x0-x-2∑x′2i 将δ0标准化,则有 δ0σ1n+x0-x-2∑x′2i~N0,1 由于σ未知,所以用σ^来代替,根据抽样分布理论及误差δ0的定义,有 y^0-Ey0σ^1n+x0-x-2∑x′2i~tn-2 那么Ey0的预测区间为 y^0-tα2·σ^1n+x0-x-2∑x′2i≤Ey0≤y^0+tα2·σ^1n+x0-x-2∑x′2i 其中α为显著水平。 在R中可以使用下面的代码来获得总体均值Ey0的预测区间。 > predict(plants.lm, + newdata = data.frame(age = 4), + interval = "confidence") fitlwrupr 1 6.59119 6.442614 6.739767 在此基础上,还可以对总体个别值y0的可能区间进行预测。设误差e0=y0-y^0,由于y^0服从正态分布,所以e0也服从正态分布。而且可以得到e0的数学期望与方差如下。 Ee0=Ey0-y^0=0 vare0=vary0-y^0 由于y^0与y0相互独立,并且 vary0=varw0+w1x0+u0=varu0 vary^0=Ey^0-Ey02=varδ0 所以 vare0=vary0+vary^0=varu0+varδ0 =σ2+σ21n+x0-x-2∑x′2i=σ21+1n+x0-x-2∑x′2i 由e0的数学期望与方差可知 e0~N0,σ21+1n+x0-x-2∑x′2i 将e0标准化,则有 e0σ1+1n+x0-x-2∑x′2i~N0,1 由于σ未知,所以用σ^来代替,根据抽样分布理论及误差e0的定义,有 y0-y^0σ^1+1n+x0-x-2∑x′2i~tn-2 那么y0的预测区间为 y^0-tα2·σ^1+1n+x0-x-2∑x′2i≤y0≤y^0+tα2·σ^1+1n+x0-x-2∑x′2i 在R中可以使用下面的代码来获得总体个别值y0的预测区间。 > predict(plants.lm, newdata = data.frame(age = 4), interval = "prediction") fitlwrupr 1 6.59119 5.877015 7.305366 可见在执行predict函数时,通过选择参数“confidence”或“prediction”即可实现对y0或者y0期望及其置信区间(或称置信带)的估计。而且y0期望的置信区间要比y0的置信区间更窄。