9  回归

9.1 引言

工程和科学中的许多问题都涉及到:确定一组变量之间的关系。例如,在化学反应中,我们可能对化学反应的产物、反应时的温度以及所用催化剂的用量之间的关系感兴趣。了解这种不同变量之间的关系可以使我们预测不同温度和催化剂用量下的反应产物。

在许多情况下,存在一个单一的 响应变量response variable)——也称之为 因变量dependent variable)—— Y因变量 Y 依赖于一组 输入变量input variables)——也称之为 自变量independent variables)——x1,,xr 的值。最简单的 因变量 Y自变量 x1,,xr 之间的关系类型是线性关系。即,对于一组常数 β0,β1,,βr,先行关系的等式为:

(9.1)Y=β0+β1x1++βrxr

如果 Yxi 之间的关系如 所示,那么一旦学习到 βi,我们就可以准确预测任意一组输入值的 Y。然而,实际上,几乎无法实现 100% 的精确预测,我们最多可以期望 在随机误差下(subject to random error)是有效的。此时, 将表示为:

(9.2)Y=β0+β1x1++βrxr+e

其中 e 表示随机误差,并假设随机误差是均值为 0 的随机变量。实际上, 的另一种表达为:

(9.3)E[Y|x]=β0+β1x1++βrxr

其中 x=(x1,,xr)自变量E[Y|x] 是在给定输入 x 下的结果的期望。

我们称 线性回归方程linear regression equation)。 描述了 因变量 Y 对一组 自变量 x1,,xr 的回归。我们称 β0,β1,,βr回归系数regression coefficients),并且通常需要通过估计得到 回归系数。如果回归方程中只包含一个 自变量——即 r=1,我们称之为 简单回归方程simple regression equation);相比之下,如果回归方程包含多个 自变量,我们称之为 多元回归方程multiple regression equation)。

因此,简单线性回归模型假设 Y 的均值与单个 自变量 的值之间存在线性关系,即:

(9.4)Y=α+βx+e

其中,x 是自变量的值,也称为 输入水平input level),Y 是对应的输出,e 代表均值为 0 的随机误差。

术语说明

在回归分析中,而对于因变量和自变量的取值会有多种描述,例如对于:

  • 因变量的取值:输出,结果,响应值(response values
  • 自变量的取值:输入,输入水平(input level

我们在后续的章节中,会根据具体的场景来穿插使用这些术语。

例 9.1 考虑以下 10 对数据 (xi,yi),i=1,,10,其中 y 是实验的百分比产量,x 是实验温度。

ixiyiixiyi1100456150682110527160753120548170764130639180925140621019088

如上数据的散点图(scatter diagram)如 所示。

代码
library(ggplot2)
x <- seq(100, 190, 10)
y <- c(45, 52, 54, 63, 62, 68, 75, 76, 92, 88)
df <- data.frame(x = x, y = y)

ggplot(df, aes(x = x, y = y)) + 
  geom_point() + 
  scale_x_continuous(breaks = seq(100, 200, 20), limits = c(100, 200)) +
  scale_y_continuous(breaks = seq(40, 100, 10), limits = c(40, 100))
图 9.1: 散点图

正如 的散点图所显示的那样,受随机误差的影响,yx 之间似乎存在线性关系。因此,看起来,对于本例的数据而言,简单线性回归模型似乎是合适的。

9.2 回归参数的最小二乘估计

假设对于输入值 xii=1,,n)的输出为 Yi,我们将这些数据作为观察值并用于估计简单线性回归模型()中的 αβ

为了确定 αβ 的估计,我们作如下推理:如果 Aα 的估计,Bβ 的估计,那么输入变量 xi 的输出估计将是 A+Bxi。由于实际输出是 Yi ,因此实际输出和输出估计之间的平方差为 (YiABxi)2,因此如果 ABαβ 的估计,那么输出估计值和实际输出值之间的平方差的和(SS: the sum of the squared)为:

(9.5)SS=i=1n(YiABxi)2

最小二乘法(the least squares)选择让 SS 达到最小值的 AB 以作为 αβ 的估计。为了确定 AB,我们首先对 SS 进行微分运算:

SSA=2i=1n(YiABxi)

SSB=2i=1nxi(YiABxi)

令偏导数为零,得到令 AB 最小化的方程:

(9.6)i=1nYi=nA+Bi=1nxii=1nxiYi=Ai=1nxi+Bi=1nxi2

我们称 正规方程the normal equations)。如果令:

Y=iYi/n,x=ixi/n

那么我们可以将第一个 正规方程 写为:

(9.7)A=YBx

A 的代入到第二个 正规方程 得到:

ixiYi=(YBx)nx+Bixi2

即,B(ixi2nx2)=ixiYinxY

即,B=ixiYinxYixi2nx2

因此,根据 nY=i=1nYi,我们可以证明

命题 9.1 在数据集 (xi,yi),i=1,,n 上的 βα 的最小二乘法估计量分别为:

B=i=1nxiYixi=1nYii=1nxi2nx2

A=YBx

我们称直线 A+Bx 为估计线性方程(the estimated regression line)。

我们可以使用 R 来获得数据对 (x1,y1),...,(xn,yn) 的估计线性方程。

x <- c(x_1, ..., x_n)
y <- c(y_1, ..., y_n)
df <- data.frame(x = x, y = y)
1fit <- lm(y ~ x, df)
2fit
1
定义模型为 lm(y~x),该模型为线性模型并且假定 yx 的线性函数加上一个随机误差,并且我们将该模型命名为 fit
2
输出模型 fit 的估计线性方程。

例如:

x <- c(1, 2, 3, 4, 5, 6, 7)
y <- c(3, 2, 5, 6, 4, 8, 9)
df <- data.frame(x = x, y = y)
fit <- lm(y ~ x, df)
fit

Call:
lm(formula = y ~ x, data = df)

Coefficients:
(Intercept)            x  
      1.143        1.036  

所以,根据如上的结果,估计线性方程为:y=1.143+1.036x

要获得估计线性方程图,首先输入 plot(x, y),这将生成散点图。然后输入 abline(fit),这将在散点图上添加估计线性方程的直线。

x <- c(1, 2, 3, 4, 5, 6, 7)
y <- c(3, 2, 5, 6, 4, 8, 9)
df <- data.frame(x = x, y = y)
fit <- lm(y ~ x, df)
plot(x, y)
abline(fit)

例 9.2 生产某种合成纤维所用的原材料储存在一个没有湿度控制的地方。该储存地点 15 天内的相对湿度和原材料样品的水分含量数据如下所示:

相对湿度(%) 46 53 29 61 36 39 47 49 52 38 55 32 57 54 44
原材料含水量(%) 12 15 7 17 10 11 11 12 14 9 16 8 18 14 12

将该线性回归模型命名为 moisture,并且利用 R 得到:

x <- c(46, 53, 29, 61, 36, 39, 47, 49, 52, 38, 55, 32, 57, 54, 44)
y <- c(12, 15,  7, 17, 10, 11, 11, 12, 14,  9, 16,  8, 18, 14, 12)
df <- data.frame(x = x, y = y)
moisture <- lm(y ~ x, df)
moisture

Call:
lm(formula = y ~ x, data = df)

Coefficients:
(Intercept)            x  
    -2.5105       0.3232  

数据的散点图和估计线性方程图为:

x <- c(46, 53, 29, 61, 36, 39, 47, 49, 52, 38, 55, 32, 57, 54, 44)
y <- c(12, 15,  7, 17, 10, 11, 11, 12, 14,  9, 16,  8, 18, 14, 12)
df <- data.frame(x = x, y = y)
moisture <- lm(y ~ x, df)

plot(x, y)
abline(moisture)

9.3 估计量的分布

为了确定估计量 AB 的分布,不能仅仅假设随机误差的均值为 0,我们有必要对随机误差做出额外的假设。通常,我们假设随机误差是均值为 0,方差为 σ2 的独立的、正态分布随机变量。也就是说,如果输入值 xi 对应的输出为 Yi,那么 Y1,,Yn 相互独立,并且:

YiN(α+βxi,σ2)

注意,上述假设认为随机误差的方差不依赖于输入值,而是一个常数。并且,我们假设 σ2 的值是未知的,我们必须通过数据对 σ2 进行估计。

根据 β 的最小二乘估计量 B 可以表示为:

(9.8)B=i(xix)Yiixi2nx2

因此,B 是独立、正态分布随机变量 Yi,i=1,,n 的线性组合,所以 B 本身也服从正态分布。根据 B 的均值和方差为:

E[B]=i(xix)E[Yi]ixi2nx2=i(xix)(α+βxi)ixi2nx2=αi(xix)+βixi(xix)ixi2nx2=βixi2xixiixi2nx2i(xix)=0=β

因此,E[B]=β,所以 Bβ 的无偏估计量。现在,我们计算 B 的方差。

(9.9)Var(B)=Var(i=1n(xix)Yi)(i=1nxi2nx2)2=i=1n(xix)2Var(Yi)(i=1nxi2nx2)2 by independence=σ2i=1n(xix)2(i=1nxi2nx2)2=σ2i=1nxi2nx2i=1n(xix)2=i=1nxi2nx2

根据 以及 A=i=1nYinBx 表明 A 也可以表示为独立正态分布随机变量 Yi,i=1,,n 的线性组合,因此 A 也服从正态分布。A 的均值为:

E[A]=i=1nE[Yi]nxE[B]=i=1nα+βxinxβ=α+βxxβ=α

因此,A 也是 α 的无偏估计量。首先将 A 表示为 Yi 的线性组合来计算 A 的方差,其结果为(详细的推导步骤作为练习留给读者):

(9.10)Var(A)=σ2i=1nxi2n(i=1nxi2nx2)

YiABxi,i=1,,n,表示实际值(即 Yi)与其最小二乘估计量(即 A+Bxi)之间的差异,我们称之为 残差residuals)。可以用残差的平方和:

SSR=i=1n(YiABxi)2

来估计随机误差的方差 σ2。实际上,可以证明

SSRσ2χn22

也就是说,SSRσ2 服从自由度为 n2 的卡方分布,这意味着

E[SSRσ2]=n2

或者

E[SSRn2]=σ2

因此,SSRn2σ2 的无偏估计量。此外,还可以证明 SSRAB 是独立的。

备注

关于为什么 SSRσ2 服从自由度为 n2 的卡方分布并且与 AB 独立的合理性论证如下。

因为 Yi 是独立正态随机变量,因此 (YiE[Yi])Var(Yi)i=1,,n)服从标准正态分布,所以:

i=1n(YiE[Yi])2Var(Yi)=i=1n(Yiαβxi)2σ2χn2

现在,如果我们用估计量 AB 替代 αβ,那么将失去 2 个自由度,因此 SSRσ2 服从自由度为 n2 的卡方分布并不令人惊讶。

事实上,SSRAB 的独立性与一个基本结果非常相似,即在正态分布抽样中,XS2 是独立的()。实际上,XS2 的独立性表明,如果 Y1,,Yn 是均值为 μ、方差为 σ2 的正态分布样本,那么平方和 i=1n(Yiμ)2/σ2 服从自由度为 n 的卡方分布。如果用均值的估计量 Y 替代 μ 得到新的平方和 i(YiY)2/σ2,那么这个量(等于 (n1)S2/σ2)将与 Y 独立,并且将服从自由度为 n1 的卡方分布。由于 SSR/σ2 是通过在平方和 i=1n(Yiαβxi)2/σ2 中用估计量 AB 来分别替代 αβ,因此这个量可能与 AB 是独立的。

Yi 是正态分布随机变量时,最小二乘估计量也是最大似然估计量。为了验证这一点,注意到 Y1,,Yn 的联合密度由下式给出:

fY1,,Yn(y1,,yn)=i=1nfYi(yi)=i=1n12πσe(yiαβxi)2/2σ2=1(2π)n/2σnei=1n(yiαβxi)2/2σ2

因此,αβ 的最大似然估计量恰好是使 i=1n(yiαβxi)2 最小的值,即最小二乘估计量。

符号说明

如果我们令:

SxY=i=1n(xix)(YiY)=i=1nxiYinxY

Sxx=i=1n(xix)2=i=1nxi2nx2

SYY=i=1n(YiY)2=i=1nYi2nY2

则最小二乘估计量可以表示为:

B=SxYSxx

A=YBx

因此,我们得到如下的残差的平方和的计算等式:

(9.11)SSR=SxxSYYSxY2Sxx

对本节的内容进行了总结。

命题 9.2 假设因变量 Yi,i=1,,n 是均值为 α+βxi、方差为 σ2 的独立的、正态分布随机变量。βα 的最小二乘估计量为:

B=SxYSxx,A=YBx

其分布为:

AN(α,σ2ixi2nSxx)

BN(β,σ2Sxx)

此外,如果我们令:

SSR=i(YiABxi)2

表示残差平方和,那么:

SSRσ2χn22

并且 SSR 与最小二乘估计量 AB 是独立的。同时,SSR 可以通过如下的方式来计算:

SSR=SxxSYY(SxY)2Sxx

可以使用 R 计算 SxYSxxSYYSSR

x <- c(x_1, ..., x_n)
y <- c(y_1, ..., y_n)
Sxy <- sum(x * y) - n * mean(x) * mean(y)
Sxx <- sum(x * x) - n * mean(x)^2
Syy <- sum(y * y) - n * mean(y)^2
SSR <- (Sxx * Syy - Sxy^2) / Sxx

习题 9.1 下面的数据中,x 为某种产品的湿混合物的水分,Y 为最终产品的密度。

xi yi
5 7.4
6 9.3
7 10.6
10 15.4
12 18.1
15 22.2
18 24.1
20 24.8

为这些数据拟合一条线性曲线,并确定回归方程的残差的平方和 SSR

解 9.1. 使用 R 计算 SSR

x <- c(5,6,7,10,12,15,18,20)
y <- c(7.4, 9.3, 10.6, 15.4, 18.1, 22.2, 24.1, 24.8)
Sxy = sum(x * y) - 8 * mean(x) * mean(y)
Sxx = sum(x * x) - 8 * mean(x)^2
Syy = sum(y * y) - 8 * mean(y)^2
SSR = (Sxx * Syy - Sxy^2) / Sxx
SSR
[1] 9.469758

使用 R 求线性回归方程:

x <- c(5,6,7,10,12,15,18,20)
y <- c(7.4, 9.3, 10.6, 15.4, 18.1, 22.2, 24.1, 24.8)
df <- data.frame(x = x, y = y)
model <- lm(y ~ x, df)
model

Call:
lm(formula = y ~ x, data = df)

Coefficients:
(Intercept)            x  
      2.463        1.206  
plot(x, y)
abline(model)
图 9.2:

所示,估计的线性回归方程为:

y=2.463+1.206x

9.4 关于回归参数的统计推断

使用 ,设计回归参数的假设检验和置信区间是一件非常简单的事情。

9.4.1 关于 β 的统计推断

β=0 是简单线性回归模型 Y=α+βx+e 的一个重要的假设。该假设的重要性源自这样的事实,即 β=0 相当于声明模型的输出均值不依赖于模型的输入,或者说没有对输入变量进行任何的回归。

为了检验:

H0:β=0vs.H1:β0

我们注意到,根据 有:

(9.12)Bβσ2/Sxx=SxxBβσN(0,1)

同时,残差平方和 SSRB 之间相互独立,且:

SSRσ2χn22

因此,根据 t 分布随机变量的定义,

(9.13)Sxx(Bβ)/σSSRσ2(n2)=(n2)SxxSSR(Bβ)tn2

也就是说,(n2)SxxSSR(Bβ) 服从自由度为 n2t 分布。因此,如果 H0 为真(即 β=0),则:

(n2)SxxSSRBtn2

H0:β=0 的假设检验 H0 在显著性水平 γ 下的检验为:

拒绝: H0(n2)SxxSSR|B|>tγ/2,n2

接受: H0

为了计算该检验,首先计算检验统计量 (n2)SXX/SSR|B| 的值(v),然后如果如下所示的 p-value 小于给定的显著性水平则拒绝 H0

p-value=P{|Tn2|>v}=2P{Tn2>v}

其中 Tn2 服从自由度为 n2 的 t 分布随机变量。

可以使用 R 来计算如上的假设检验。例如,如果模型的名称为 fit,那么 summary(fit) 将提供所需的 p-value

习题 9.2 一个人声称他的汽车的燃油消耗量不依赖于汽车的行驶速度。为了检验这个假设的合理性,在不同的速度(45~70 英里/小时)测试了汽车。不同速度下的每加仑英里的测试数据如下:

速度 每加仑英里数
45 24.2
50 25.0
55 23.3
60 22.0
65 21.5
70 20.6
75 19.8

这些数据是否拒绝了每加仑汽油的英里数不受汽车行驶速度影响的说法?

解 9.2. 假设 Y 为汽车每加仑汽油的英里数,x 为汽车的行驶速度,并且 Yx 符合简单线性回归模型:

Y=α+βx+e

因此,题目中的假设意味着回归系数 β=0

使用 R,我们有:

x <- c(45, 50, 55, 60, 65, 70, 75)
y <- c(24.2, 25, 23.3, 22, 21.5, 20.6, 19.8)
miles <- lm(y ~ x)
miles

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
      32.54        -0.17  
summary(miles)

Call:
lm(formula = y ~ x)

Residuals:
        1         2         3         4         5         6         7 
-0.692857  0.957143  0.107143 -0.342857  0.007143 -0.042857  0.007143 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 32.54286    1.27059  25.612 1.69e-06 ***
x           -0.17000    0.02089  -8.138 0.000455 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5527 on 5 degrees of freedom
Multiple R-squared:  0.9298,    Adjusted R-squared:  0.9158 
F-statistic: 66.23 on 1 and 5 DF,  p-value: 0.0004548

summary(miles) 首先给出七个数据的残差 yi(A+Bxi),其中 i=1,,7。接下来的几行给出了截距 α 和斜率 β 的估计值分别为 32.54286 和 -0.17。然后,给出了这些估计值的标准误差。接着给出了用于检验回归参数为零的 t-统计量的值。例如,检验 β=0t-值为 -8.138。随后给出了回归参数为零的检验的 p-value,例如:α=0p-value1.69×106β=0p-value 为 0.000455,这表明在实践中,无论显著性水平为何值,两个假设都将被拒绝。残差标准误差的值为 SSRn2。由于 SSR/(n2)σ2 的估计量,因此残差标准误差是 σ 的估计量。

根据 可以很容易地得到 β 的置信区间估计。实际上,从 可以推导出,对于任意 a,其中 0<a<1,有

P{ta/2,n2<(n2)SxxSSR(Bβ)<ta/2,n2}=1a

即:

P{BSSR(n2)Sxxta/2,n2<β<B+SSR(n2)Sxxta/2,n2}=1a

于是,我们得到了以下的结果。

β 的置信区间

β100(1a)% 的置信区间估计是:

(BSSR(n2)Sxxta/2,n2,B+SSR(n2)Sxxta/2,n2)

备注

因为 σ2 未知,所以不能直接使用 Bβσ2/SxxN(0,1) 来对 β 进行统计推断。相反,我们采用了 中提到的 σ2 的估计量 SSR/(n2),这导致统计量的分布从标准正态分布变为自由度为 n2t 分布。

例 9.3 我们还可以使用 R 来获得 αβ 的置信区间。如果模型的名称是 name,那么 R 命令 confint(name, level = m) 会返回模型 name 中的 αβ100m% 置信区间。例如,假设我们想要计算模型 miles 的参数的 95% 置信区间。如果我们输入 confint(miles, level = 0.95),那么 R 会返回 αβ 的 95% 置信区间。

x <- c(45, 50, 55, 60, 65, 70, 75)
y <- c(24.2, 25, 23.3, 22, 21.5, 20.6, 19.8)
miles <- lm(y ~ x)
confint(miles)
                 2.5 %     97.5 %
(Intercept) 29.2766922 35.8090220
x           -0.2236954 -0.1163046

所以 αβ 的 95% 置信区间分别为:(29.2767,35.8090)(0.2237,0.1163)

9.4.2 均值回归

术语 回归regression)最初由 Francis Galton 在描述遗传定律时使用。Galton 认为这些定律导致了种群中的极端值“向均值回归”(regress toward the mean)。在“向均值回归”中,Galton 的意思是:在某个特征上具有极端值的人的孩子会倾向于比他们的父母具有更不极端的值。

如果我们假设在某特征上,后代 (Y) 和父代 (x) 之间存在线性回归关系,那么当回归参数 β 介于 0 和 1 之间时,将发生均值回归。也就是说,如果

E[Y]=α+βx

并且 0<β<1,那么当 x 很大时,E[Y] 将小于 x;当 x 很小时,E[Y] 将大于 x。可以通过代数方法或者绘制两条直线 y=α+βxy=x 来轻松验证如上的陈述。

在直线图中我们可以发现,当 0<β<1 时,对于较小的 x 值,y=α+βxy=x 的上方,而对于较大的 x 值,则位于其下方。

例 9.4 为了演示 Galton 的均值回归理论,英国统计学家 Karl Pearson 绘制了 10 名随机选择的男孩与其父亲的身高的对比图。具体数据(以英寸为单位)如下:

父亲的身高:60, 62, 64, 65, 66, 67, 68, 70, 72, 74

儿子的身高:63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70

如上数据的散点图如 所示。

代码
father_height <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
son_height <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
plot(father_height, son_height)
图 9.3: 父子之间的身高对比图

注意,虽然数据表明身高较高的父亲其儿子的身高也往往较高,但也表明较高或者较低身高的父亲的儿子的身高往往比他们更“平均”——也就是说,孩子的身高有“向均值回归”的现象。

我们将通过将 均值回归 作为 备则假设 来确定上述数据是否足以证明 均值回归 的存在。也就是说,我们将使用上述数据来检验:

H0:β1vsH1:β<1

即:

H0:β=1vsH1:β<1

根据 ,当 β=1 时,检验统计量:

TS=8SxxSSR(B1)t8

如果 TS=v,则原假设 H0:β1p-value 为:

p-value=P{T8v}

其中,T8 为自由度为 8 的 t-分布。使用 R 来计算如上的值:

x <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
y <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
height = lm(y ~ x)

summary(height)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6738 -0.2374 -0.0852  0.2835  0.6742 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 35.97681    2.20760   16.30 2.02e-07 ***
x            0.46457    0.03298   14.08 6.27e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4321 on 8 degrees of freedom
Multiple R-squared:  0.9612,    Adjusted R-squared:  0.9564 
F-statistic: 198.4 on 1 and 8 DF,  p-value: 6.273e-07

由此可知,β 的估计值 B=0.46457,并且 β 的估计的标准误差为 0.03298,β=0 时的统计检验 8SxxSSRB 为 14.08。因此,对于 H0:β=1 时的统计检验值 TS=8SxxSSR(B1)=14.0814.080.46457=16.2276

pt(-16.2276, 8)
[1] 1.045567e-07

因此,统计检验量的值是 -16.2276,相应的 p-value1.045569×107。因此,在任何显著性水平上,原假设 $$ 将被拒绝,从而建立了均值回归。(参见 ,在 的散点图上增加了直线 y=x)。

代码
father_height <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
son_height <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
plot(father_height, son_height, xlim=c(60, 74), ylim=c(60, 74))
abline(a=0, b=1, col="red")
图 9.4: 对于比较小的 x,y > x;对于比较大的 x,y < x

一种现代生物学解释认为,均值回归现象大致可以解释为:后代会随机获取父母的一半的基因,因此,假设某人的身高非常高,则其后代很可能由于偶然原因拥有了更少的“高个子”基因。

尽管均值回归现象最重要的应用是关于后代与其父母的生物学特征之间的关系,但这一现象在我们拥有指向相同变量的两组数据场景时也会出现。

例 9.5  

表 9.1: 美国西北部 12 个县在 1988 年和 1989 年的机动车死亡人数
县市 1988 年死亡人数 1989 死亡人数
1 121 104
2 96 91
3 85 101
4 113 110
5 102 117
6 118 108
7 90 96
8 84 102
9 107 114
10 112 96
11 95 488
12 101 106

的数据涉及美国西北部 12 个县在 1988 年和 1989 年发生的机动车死亡人数。

代码
death_1988 <- c(121, 96, 85, 113, 102, 118, 90, 84, 107, 112, 95, 101)
death_1989 <- c(104, 91, 101, 110, 117, 108, 96, 102, 114, 96, 488, 106)
plot(death_1988, death_1989, xlim = c(80, 130), ylim = c(80, 130))
图 9.5: 1988 年和 1989 年的死亡人数对比图

可以看出,1989 年大多数县的死亡人数减少了,这些县在 1988 年有大量的机动车死亡人数。同样,那些在 1988 年死亡人数较少的县,1989 年的死亡人数似乎有所增加。因此,我们可以认为均值回归现象确实存在。实际上,可以用 R 计算 的估计回归方程:

y=74.589+0.276x

这也表明,β 的估计值确实小于 1。

在考虑上述数据中的均值回归现象的原因时,我们必须要小心。例如,我们会很自然的假设,在 1988 年机动车死亡人数较多的县可能通过改善道路安全或让人们更多地意识到不安全驾驶的潜在危险来做出大力改进,以减少机动车死亡人数。同时,我们可能假设在 1988 年死亡人数最少的那些县可能“沾沾自喜”,并且没有做出更多的努力来进一步改善机动车的死亡人数,因此导致第二年死亡人数的增加。

尽管这种假设可能是正确的,但重要的是要认识到,即使没有任何一个县做出了任何非同寻常的举措,也可能会发生均值回归现象。事实上,很有可能在 1988 年那些死亡人数较多的县在那一年非常不幸,因此下一年的减少只是其更正常结果的回归。例如,如果抛 10 次硬币,如果出现了 9 次正面朝上,那么再抛这枚硬币 10 次,则正面向上的次数很可能会少于 9 次。同样,那些在 1988 年死亡人数较少的县在那一年可能非常“幸运”,而 1989 年的更正常结果将导致死亡人数的增加。

错误地认为均值回归是由于某种外部影响导致,但实际上均值回归只是由于统计中的随机性(chance occurs frequently)而导致的,这通常被称为 回归谬误the regression fallacy)。

9.4.3 关于 α 的统计推断

确定 α 的置信区间和假设检验的方法与 β 的方法完全相同。具体而言,可利用 证明:

(9.14)n(n2)SxxSSRixi2(Aα)tn2

α 的置信区间估计

α100(1α)% 的置信区间为:

A±tα/2,n2SSRixi2n(n2)Sxx

关于 α 的假设检验可以很容易地由 得到,其推导过程将留作习题。

9.4.4 关于 α+βx0 的统计推断

通常我们对利用数据对 (xi,Yi),i=1,,n 来估计 α+βx0 感兴趣(在给定输入 x0 下估计平均输出)。

对于点估计量而言,那么我们自然会想到 A+Bx0,因为

E[A+Bx0]=E[A]+x0E[B]=α+βx0

所以 A+Bx0α+βx0 的无偏估计量。

然而,如果我们需要一个置信区间,或对 α+βx0 进行某种假设检验,那么首先必须确定估计量 A+Bx0 的概率分布。

使用 中的 B 的表达式,可以得到:

B=c(i=1n(xix)Yi)

其中,c=1i=1nxi2nx2=1Sxx

由于 A=YBx,所以

A+Bx0=i=1nYinB(xx0)=i=1nYi[1nc(xix)(xx0)]

由于 Yi 是独立的正态分布随机变量,上述等式表明 A+Bx0 可以表示为独立正态分布随机变量的线性组合,因此上述等式本身也服从正态分布。我们已经知道了 A+Bx0 的均值,现在只需要计算它的方差:

Var(A+Bx0)=i=1n[1nc(xix)(xx0)]2Var(Yi)=σ2i=1n[1n2+c2(xix0)2(xix)22c(xix)(xx0)n]=σ2[1n+c2(x0x)2i=1n(xix)22c(xx0)i=1n(xix)n]i=1n(xix)2=1c=Sxx,i=1n(xix)=0=σ2[1n+(xx0)2Sxx]

于是,我们得到:

(9.15)A+Bx0N(α+βx0,σ2[1n+(xx0)2Sxx])

根据 SSRσ2χn22 且其与 A+Bx0 独立,有:

(9.16)A+Bx0(α+βx0)1n+(x0x)2SxxSSRn2tn2

现在,可以利用 得到 α+βx0 的置信区间估计。

α+βx0 的置信区间估计

100(1α)% 的置信度下,α+βx0 将位于:

A+Bx0±(1n+(x0x)2SxxSSRn2)tα/2,n2

习题 9.3 使用 的数据,确定在 95% 的置信度下,父亲身高为 68 英寸下的所有男性的平均身高的置信区间。

解 9.3. 我们使用 R 来计算该置信区间:

x <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
y <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
Sxy <- sum(x * y) - 10 * mean(x) * mean(y)
Sxx <- sum(x * x) - 10 * mean(x)^2
Syy <- sum(y * y) - 10 * mean(y)^2
Syy
[1] 38.529
SSR <- (Sxx * Syy - Sxy^2) / Sxx 
SSR
[1] 1.493578
B <- Sxy / Sxx
A <- mean(y) - B * mean(x)
a <- sqrt((1 / 10 + (68 - mean(x))^2 / Sxx) * SSR / 8)
l <- (A + B * 68) - a * qt (0.975, 8)
u <- (A + B * 68) + a * qt (0.975, 8)
c(l,u)
[1] 67.23944 67.89552

所以,95% 置信区间为:α+βx0(67.23944,67.89552)

9.4.5 输出的预测区间

通常情况下,估计输出的实际值比估计输出的均值更重要。例如,如果在温度为 x0 时进行实验,那么我们可能对预测 Y(x0) 更感兴趣(即实验的实际输出),而不是估计输出的期望 E[Y(x0)]=α+βx0。另一方面,如果在温度为 x0 时进行一系列实验,那么我们可能会更想估计输出的期望值,即 α+βx0

首先,假设我们对单个值(而不是区间)感兴趣,并以此作为 Y(x0) 的预测值,即在输入为 x0 下的输出值。显然,Y(x0) 的最佳预测值是其均值 α+βx0。实际上,这并不是那么显而易见,因为人们可能认为随机变量的最佳预测值是:

  1. 均值——这最小化了预测值与实际值之间的平方差;
  2. 中位数——这最小化了预测值与实际值之间的绝对差;
  3. 众数——最可能出现的值。

然而,由于正态分布随机变量的均值、中位数和众数相等,并且输出服从正态分布,因此此时如上的 3 种估计量并没有区别。由于 αβ 是未知的,使用其估计量 AB 似乎是合理的,因此可以使用 A+Bx0 作为新的输入 x0 的新输出预测值。

假设我们现在不关心确定某个值来预测输出,而是希望找到一个预测区间,以使得在给定的置信度下,该预测区间将包含输出值。为了获得这样的区间,让 Y 表示输入为 x0 时的输出预测值,并考虑在输出中减去其预测值的概率分布,即 YABx0。因为,

YN(α+βx0,σ2)

并且,如 所示,

A+Bx0N(α+βx0,σ2[1n+(x0x)2Sxx])

因此,由于 Y 与用于确定 AB 的数据值 Y1,Y2,,Yn 独立,所以 YA+Bx0 独立,因此

YABx0N(0,σ2[1+1n+(x0x)2Sxx])

即,

(9.17)YABx0σn+1n+(x0x)2SxxN(0,1)

因为 SSRAB 独立(当然,也与 Y 独立),并且

SSRσ2χn22

于是得到,

YABx0n+1n+(x0x)2SxxSSRn2tn2

因此,对于任意的 a0<a<1

P{ta/2,n2<YABx0n+1n+(x0x)2SxxSSRn2<ta/2,n2}=1a

即我们构建了以下的结论。

在输入为 x0 处时输出的预测区间

基于输入值 xi 对应的输出值 Yii=1,2,,n:在 100(1a)% 的置信度下,x0 处的输出值 Y 将包含在如下的区间内

A+Bx0±ta/2,n2n+1n+(x0x)2SxxSSRn2

例 9.6 中,假设我们希望得到一个区间,并且我们“有 95% 的信心”认为该区间将包含父亲的身高为 68 英寸的男性的身高。简单计算后,现在得到预测区间

Y(68)67.568±1.050

即,在 95% 的置信度下,父亲身高为 68 英寸的男性的身高将在 66.518 和 68.618 之间。

备注
  • 通常,置信区间预测区间 之间的区别会存在一些混淆。置信区间是在给定置信度下包含固定参数的区间。预测区间是在给定置信度下包含随机变量的区间。
  • 如果待预测的输入数据与获得回归方程的数据的差异较大,那么我们不应该使用该回归方程来预测该输入数据的输出。例如,不应使用 的数据获得的回归方程来预测父亲身高为 42 英寸的男性的身高。
表 9.2: 计算不同推断对象时所用的分布
推断对象 使用的分布
β (n2)SxxSSR(Bβ)tn2
α n(n2)Sxxxi2SSR(Aα)tn2
α+βx0 A+Bx0αβx01n+(x0x)2SxxSSRn2tn2
Y(x0) Y(x0)ABx01+1n+(x0x)2SxxSSRn2tn2
置信区间与预测区间的不同

置信区间对应的是回归方程输出的平均值即 α+βx0)的置信度区间,而预测区间是回归方程的输出值 α+βx0+e)(而不是平均值)的置信区间。

根据 可知,Y(x0) 的区间比 α+βx0 的区间要宽。更直观的,我们可以使用 R 计算 的回归方程的置信区间和预测区间。

x <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
y <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
df <- data.frame(x = x, y = y)
fit <- lm(y ~ x, df)

new_x <- seq(min(x), max(x), length.out = 100)
conf_interval <- predict(fit, newdata = data.frame(x = new_x), interval = "confidence")
pred_interval <- predict(fit, newdata = data.frame(x = new_x), interval = "prediction")

plot(x, y, main = "The Confidence Interval & Predict Interval of the Regression")
abline(fit, col = "blue", lwd = 2)
lines(new_x, conf_interval[, "lwr"], col = "red", lwd = 1, lty = 2)
lines(new_x, conf_interval[, "upr"], col = "red", lwd = 1, lty = 2)
lines(new_x, pred_interval[, "lwr"], col = "green", lwd = 1, lty = 3)
lines(new_x, pred_interval[, "upr"], col = "green", lwd = 1, lty = 3)
legend("topleft", legend = c("Regression Line", "Confidence Interval", "Prediction Interval"),
  col = c("blue", "red", "green"), lty = c(1, 2, 3), lwd = 2)

9.4.6 总结

  • 模型Y=α+βx+e,其中 eN(0,σ2)
  • 数据(xi,Yi)i=1,2,,n

9.5 决定系数和样本相关系数

变异量

变异量the amount of variation)是指一组数据中的变异性,即数据的波动或分散程度。这种波动可以是由多种因素引起的,包括个体差异、测量误差、随机误差等。在后续的翻译章节中,我们会根据具体的场景交叉使用 变异量数据波动,如无特殊说明,这些指的都是一个概念。

同时,我们会避免使用 数据差异difference) 来指代 变异量 的概念,因为在数学中,数据差异 更多指代两个数据之间的差值,因此和我们所讨论的数据的波动和分布程度并不是一个概念。

针对一组输入值 x1,,xn,假设我们想要测量这组输入值对应的响应值(response valueY1,,Yn 之间的数据波动(或 变异量)。在统计学中,测量 Y1,,Yn 之间的数据波动(变异量)的标准方法为:

SYY=i=1n(YiY)2

例如,如果所有的 Yi 相等,那么 Yi 都等于 Y,所以 SYY 将等于 0。

Yi 之间的波动来源于两个因素。

  1. 首先,由于输入值 xi 是不同的,因此,响应变量 Yi 具有不同的响应均值,这将导致 Yi 的值存在一定的波动。
  2. 其次,波动还来自这样一个事实:除了输入值 xi 的差异外,每个响应变量 Yi 都具有 σ2 的方差。因此,Yi 不会和xi 处的预测值完全相等。

让我们考虑这样一个问题:响应值的波动在多大程度上是由不同输入值引起的,又在多大程度上是由响应值的固有方差引起的(已经考虑了输入值的差异)。为了解答这个问题,我们关注到 SSR 的定义:

SSR=i=1n(YiABxi)2

因此,SSR 度量了响应变量中,除了不同输入值带来的变异量之外的、剩余的变异量。

因此,SYYSSR 表示:响应变量中,由不同输入值能够解释的变异量。因此,我们令 R2 表示为:响应变量中,由不同输入值能够解释的变异量比例。我们称 R2决定系数coefficient of determination),其定义为:

R2=SYYSSRSYY=1SSRSYY

解释变异量

解释变异量Explained Variation)是回归分析中的一个重要概念,用于衡量回归模型中自变量能够解释的因变量的变异量。换句话说,解释变异量是由回归模型的预测结果解释的因变量的总变异的一部分。

决定系数 R2 的取值范围在 0 和 1 之间。接近 1 的 R2 表示响应数据的大部分波动性可以由不同的输入值来解释,而接近 0 的 R2 则表示响应数据的波动几乎不能够由不同的输入值来解释。

例 9.7 中的数据揭示了父子身高之间的关系,我们可以使用 R 计算得到:

x <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
y <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
Sxy <- sum(x * y) - 10 * mean(x) * mean(y)
Sxx <- sum(x * x) - 10 * mean(x)^2
Syy <- sum(y * y) - 10 * mean(y)^2
SSr <- (Sxx * Syy - Sxy^2) / Sxx
Syy
[1] 38.529
SSr
[1] 1.493578
R2 = 1 - SSr / Syy
R2
[1] 0.961235

我们可以使用 R 中的 summary() 来获取 R2 的计算结果。

x <- c(60, 62, 64, 65, 66, 67, 68, 70, 72, 74)
y <- c(63.6, 65.2, 66, 65.5, 66.9, 67.1, 67.4, 68.3, 70.1, 70)
df <- data.frame(x = x, y = y)
fit <- lm(y ~ x, df)

summary(fit)

Call:
lm(formula = y ~ x, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6738 -0.2374 -0.0852  0.2835  0.6742 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 35.97681    2.20760   16.30 2.02e-07 ***
x            0.46457    0.03298   14.08 6.27e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4321 on 8 degrees of freedom
Multiple R-squared:  0.9612,    Adjusted R-squared:  0.9564 
F-statistic: 198.4 on 1 and 8 DF,  p-value: 6.273e-07

根据如上结果,我们可以看到:Multiple R-squared: 0.9612Multiple R-squaredR2。换句话说,在 10 个人的身高波动中,有 96% 的身高波动是由其父亲的身高带来的(解释的),剩下的(未解释的)4% 的数据波动是由于儿子的身高波动带来的(即随机误差的方差 σ2)。

通用用 决定系数 R2 的值作为回归模型对数据拟合好坏的指标,R2 值越接近 1 则表示拟合效果越好,而 R2 值越接近 0 则表示拟合效果越差。换句话说,如果回归模型能够解释响应数据中的绝大部分波动(variation in the response data),那么就认为该模型可以较好的拟合数据。

回想一下我们在 中定义的数据对 (xi,Yi),i=1,,n样本相关系数 r

r=i=1n(xix)(YiY)i=1n(xix)2i=1n(YiY)2

样本相关系数 r 提供了一种衡量 xY 之间相关性强度的方式。当 r 接近 +1 时,表示较大的 x 对应着较大的 Y,并且较小的 x 对应着较小的 Y;当 r 接近 -1 时,表示较大的 x 对应着较小的 Y,并且较小的 x 对应着较大的 Y

在本章中,r=SxYSxxSYY,根据 有:SSR=SxxSYYSxY2Sxx,于是有:

r2=SxY2SxxSYY=SxxSYYSSRSxxSxxSYY=1SSRSYY=R2

即:

|r|=R2

因此,除了 r 的符号表示它是正相关还是负相关外,样本相关系数决定系数 的平方根是一致的。r 的符号与参数 B 的符号相同。

如上的内容为 样本相关系数 提供了额外的意义。例如,如果一个数据集的 样本相关系数 r 等于 0.9,那么这意味着对于这些数据的 简单线性回归模型 可以解释 响应变量 81% 的数据波动(variation)(因为 R2=0.92=0.81)。换句话说,对于响应变量的数据波动而言,81% 可以由不同的输入值来解释。

9.6 残差分析:评估模型

可以通过分析数据的散点图来初步确定简单线性回归模型

Y=α+βx+e,e(0,σ2)

是否适用于特定场景。事实上,数据散点图的分析通常足以让我们确认回归模型是否正确。当散点图本身不能确定简单线性模型是否合适时,就需要计算模型的 最小二乘估 计量 AB,并分析残差 Yi(A+Bxi), i=1,,n。残差分析(analysis of residuals)需要首先对残差 正态化normalizing)或者 标准化standardizing),残差除以 SSR/(n2)(即 Yi 的标准差的估计值,)得到的指标量

Yi(A+Bxi)SSR/(n2),i=1,,n

称为 标准化残差Standardized Residuals)。

当简单回归模型正确时,标准化残差 应近似为相互独立的、标准正态分布随机变量,因此 标准化残差 应围绕 0 随机分布,且有大约 95% 的 标准化残差 应位于 (2,2) 区间内(因为 P{1.96<Z<1.96}=0.95)。此外,标准化残差图不应显示出任何明显的模式(pattern)。实际上,任何明显的模式都暗示着我们需要怀疑假定的简单线性回归模型的有效性。

展示了三个不同的散点图及其相关的标准化残差图。 中的散点图及其标准化残差的随机性表明线性模型的拟合比较好。 中的残差图显示出一个明显的模式(pattern),即残差随着输入变量值的增加时呈现出先减小后增加的趋势。这通常意味着需要更高阶(不仅仅是线性)的项来描述输入和响应之间的关系。在散点图中,也体现了如上的结论。 中的标准化残差图也显示出一种模式,即残差的值及其平方随着输入变量值的增加而增加。这通常意味着 Yi 的方差并不是常数,而是随着输入变量值的增加而增加。

(a) A
(b) B
(c) C
图 9.6

9.7 转换为线性回归

在许多情况下,响应均值(输出值的均值)明显不是输入水平(输入值)的线性函数。在这种情况下,如果可以确定输入水平和其响应均值之间的关系形式,那么有时可以通过变量变换将其转换为线性形式。例如,在某些应用中,已知信号在 t 时刻处的振幅 W(t) 近似满足:

W(t)cedt

对上式去对数运算,则可以表示为:

logW(t)logcdt

如果令 Y=logW(t)α=logcβ=d,则模型 W(t)cedt 可以转换为:

Y=α+βt+e

然后,我们可以利用最小二乘法来估计简单线性回归模型的回归参数 αβ,于是 W(t)cedt 可以用如下方式来预测:

W(t)eA+Bt

习题 9.4 下表给出了在不同温度(摄氏度)下进行实验时所使用的化学物质的百分比。使用这些数据估计:如果实验在 350 度下进行,那么需要使用的化学物质的百分比是多少?

温度 百分比
0.061
10° 0.113
20° 0.192
30° 0.259
40° 0.339
50° 0.401
60° 0.461
80° 0.551

解 9.4.

代码
x <- c(5, 10, 20, 30, 40, 50, 60, 80)
y <- c(0.061, 0.113, 0.192, 0.259, 0.339, 0.401, 0.461, 0.551)
plot(x, y)
图 9.7:

P(x) 为在温度 10x 时所使用的化学物质的百分比。尽管如 所示,P(x) 看起来大致呈线性关系,但我们可以通过考虑 xP(x) 之间的非线性关系来改进模型的拟合效果。具体而言,我们考虑以下形式的关系:

1P(x)c(1d)x

也就是说,假设在温度为 x 的实验中,所用化学物质的百分比大约以指数速率减少。对上式进行对数运算后,可以写成:

log(1P(x))log(c)+xlog(1d)

令:

Y=log(1P)α=logcβ=log(1d)

我们得到了一般的回归方程:

Y=α+βx+e

为了看出数据是否支持这个模型,我们需要画出 xlog(1P) 的散点图。 给出了转换后的数据, 给出了对应的散点图。

表 9.3
温度 log(1P)
0.063
10° 0.120
20° 0.213
30° 0.300
40° 0.414
50° 0.512
60° 0.618
80° 0.801
代码
x <- c(5, 10, 20, 30, 40, 50, 60, 80)
y <- c(0.061, 0.113, 0.192, 0.259, 0.339, 0.401, 0.461, 0.551)
y1 <- -log(1 - y)
plot(x, y1)
图 9.8

使用 R 得到 αβ 的最小二乘估计值():

x <- c(5, 10, 20, 30, 40, 50, 60, 80)
y <- c(0.061, 0.113, 0.192, 0.259, 0.339, 0.401, 0.461, 0.551)
y1 <- -log(1 - y)

B <- (sum(x * y1) - mean(x) * sum(y1)) / (sum(x^2) - length(x) * mean(x)^2)
A <- mean(y1) - B * mean(x)

B
[1] 0.009889817
A
[1] 0.01544615

于是得到 cd 的估计值:

c^=eA=0.98471d^=eB=0.9901

最后得到估计的函数关系为:

P^=10.9847(0.9901)x

如上模型的残差 PP^ 所示。

表 9.4
x P P^ PP^
5 0.061 0.063 -0.002
10 0.113 0.109 0.040
20 0.192 0.193 -0.001
30 0.259 0.269 -0.010
40 0.339 0.339 0.000
50 0.401 0.401 0.000
60 0.461 0.458 0.003
80 0.551 0.556 -0.005

9.8 加权最小二乘法

在回归模型

Y=α+βx+e

中,我们经常发现 响应变量 的方差不是常数,而是依赖于其 输入水平input level)。如果已经知道这些方差——至少应该知道这些方差的系数——那么可以通过最小化加权平方和来估计回归参数 αβ。具体来说,如果

Var(Yi)=σ2wi

那么估计量 AB 应该为下式的最小化:

i[Yi(A+Bxi)]2Var(Yi)=1σ2iwi(YiABxi)2

分别对 AB 进行偏微分计算,然后令偏微分为 0,得到令 AB 最小化的方程:

(9.18)iwiYi=Aiwi+BiwixiiwixiYi=Aiwixi+Biwixi2

很容易就可以求解这些方程以得到最小二乘估计量。

例 9.8 为了理解为什么应该通过最小化加权平方和而不是普通平方和来获得估计量,需要考虑以下情况。假设 X1,,Xn 是具有均值 μ 和方差 σ2 的独立正态分布随机变量。假设 Xi 不可直接观测,而只能观测到

Y1=X1++Xk,Y2=Xk+1++Xn,k<n

我们如何基于 Y1Y2 估计 μ

虽然 μ 的最佳估计量显然是 X=i=1nXi/n=(Y1+Y2)/n,但让我们看看普通最小二乘估计量会是什么。因为

E[Y1]=kμ,E[Y2]=(nk)μ

μ 的最小二乘估计量是使得以下表达式最小化的 μ 的值:

(Y1kμ)2+(Y2[nk]μ)2

对其求导并令其等于零,我们得到最小二乘估计量 μ(记作 μ^)应满足:

2k(Y1kμ^)2(nk)[Y2(nk)μ^]=0

即:

[k2+(nk)2]μ^=kY1+(nk)Y2

所以:

μ^=kY1+(nk)Y2k2+(nk)2

因为

E[μ^]=kE[Y1]+(nk)E[Y2]k2+(nk)2=k2μ+(nk)2μk2+(nk)2=μ

因此,我们看到虽然普通最小二乘估计量是 μ 的无偏估计量,但它却不是 X 的最佳估计量。

现在,我们通过最小化加权平方和来确定 μ 的估计量。也就是说,让我们确定 μ 的值(记作 μw)以使得以下表达式最小化:

(Y1kμ)2Var(Y1)+[Y2(nk)μ]2Var(Y2)

因为

Var(Y1)=kσ2,Var(Y2)=(nk)σ2

这相当于需要选择令下式最小化的 μ

(Y1kμ)2k+[Y2(nk)μ]2nk

对其求导并令其等于零,我们得到 μw 应满足:

2k(Y1kμw)k2(nk)[Y2(nk)μw]nk=0

即:

Y1+Y2=nμw

所以:

μw=Y1+Y2n

即加权最小二乘估计量实际上就是首选的估计量 (Y1+Y2)/n=X

备注
  1. 假设数据是正态分布的,加权最小二乘估计量就是最大似然估计量。这是因为数据 Y1,,Yn 的联合概率密度为:

    fY1,,Yn(y1,,yn)=i=1n12π(σ/wi)e(yiαβxi)2/(2σ2/wi)=w1wn(2π)n/2σnei=1nwi(yiαβxi)2/2σ2

    因此,αβ 的最大似然估计量正是使的加权平方和 i=1nwi(yiαβxi)2 最小的 αβ 的值。

  2. 加权平方和也可以看作是回归方程 Y=α+βx+e 乘以 w,即:

    Yw=αw+βxw+ew

    在这个方程中,误差项 ew 的均值为 0,且方差是常数。因此,αβ 的最小二乘估计量应是使以下表达式最小的 AB 的值:

    i(YiwiAwiBxiwi)2=iwi(YiABxi)2

  3. 加权最小二乘法主要关注那些具有最大权重的数据对(会让权重最大的数据项的误差的方差最小)。

因为加权最小二乘法需要知道在任意输入水平处的响应变量的方差(至少需要知道其系数),因此加权最小二乘法可能看起来不太有用。但是,通常可以通过分析模型生成的数据来确定这些值。以下两个例子将对此进行说明。

习题 9.5 下表表示某城市市区内的通勤时间。自变量或输入变量是通勤的距离。

通勤距离(英里) 0.5 1 1.5 2 3 4 5 6 8 10
通勤时间(分钟) 15.0 15.1 16.5 19.9 27.7 29.7 26.7 35.9 42 49.4

假设通勤时间 Y 和通勤距离 x 之间的关系为线性关系:

Y=α+βx+e

我们应如何估计 αβ

为了使用加权最小二乘法,我们需要知道 Y 的方差。现在,我们提出一个论点,即 Var(Y) 应与 x 成正比。

解 9.5. d 表示一个街区的长度,于是一个通勤距离为 x 的行程将由 x/d 个街区组成。如果我们让 Yi 表示第 i 个街区所需的时间,那么总通勤时间可以表示为:

Y=Y1+Y2++Yx/d

在许多应用中,合理的假设是 Yi 是具有相同方差的独立随机变量,因此

Var(Y)=Var(Y1)++Var(Yx/d)=(x/d)Var(Y1)=xσ2

其中 σ2=Var(Y1)/d

因此,估计量 AB 应该令以下表达式最小化:

i(YiABxi)2xi

使用上面的数据,且权重为 wi=1/xi 的最小二乘方程为:

104.22=5.34A+10B277.9=10A+41B

求解得:

A=12.561,B=3.714

展示了数据点和其估计回归线 12.561+3.714x 的图形。对我们的解的定性分析可以知:在输入水平较小时,回归的效果最好,这主要是因为权重与输入水平成反比。

代码
miles <- c(0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10)
time <- c(15.0, 15.1, 16.5, 19.9, 27.7, 29.7, 26.7, 35.9, 42, 49.4)
plot(miles, time)
abline(a = 12.561, b = 3.714)
图 9.9

例 9.9 考虑在交通繁忙的高速公路上发生的事故数量 Y 与在高速公路上行驶的车辆数量 x 之间的关系。稍微思考后,大多数人可能会认为其比较合适的关系应该为线性模型:

Y=α+βx+e

然而,由于没有任何先验(priori)理由认为 Var(Y) 和输入水平 x 之间无关,因此我们不清楚是否可以使用普通最小二乘法来估计 αβ。实际上,我们将论证:应该采用加权最小二乘法且权重为 1/x,即我们应该选择 AB 以最小化

i(YiABxi)2xi

背后的理由是,假设 Y 近似服从泊松分布是合理的。因为我们可以想象每辆车 x 都存在一个小概率会导致事故,因此对于较大的 x,事故的数量应该近似为泊松随机变量。由于泊松随机变量的方差等于其均值,我们可以看到:

Var(Y)E[Y]=α+βxβx

备注
  1. 在响应变量的方差依赖于输入水平时,另一种常用的技术是通过适当的变换来稳定方差。例如,如果 Y 是一个均值为 λ 的泊松随机变量,那么可以证明(参见第 2 点),无论 λ 的值是多少 Y 的近似方差都是 0.25。基于这一事实,可以尝试将 E[Y] 建模为输入水平的线性函数,即:

    Y=α+βx+e

  2. Y 是具有均值 λ 的泊松分布时,证明 Var(Y)0.25。考虑 g(y)=yλ 处的泰勒级数(Taylor series)展开。忽略二阶导数项以后的所有项,我们得到:

    (9.19)g(y)g(λ)+g(λ)(yλ)+g(λ)(yλ)22

    由于

    g(λ)=12λ1/2,g(λ)=14λ3/2

    于是,在 y=Y 处,我们得到

    Yλ+12λ1/2(Yλ)18λ3/2(Yλ)2

    因为 E[Yλ]=0E[(Yλ)2]=Var(Y)=λ,所以:

    E[Y]λ18λ

    因此:

    (E[Y])2λ+164λ14λ14

    因此:

    Var(Y)=E[Y](E[Y])2λ(λ14)=14

9.9 多项式回归

当无法通过线性关系充分拟合响应变量 Y 和自变量 x 之间的函数关系的情况下,有时可以通过考虑多项式关系来获得合理的拟合。也就是说,可以尝试将数据集拟合到以下形式的函数关系:

Y=β0+β1x+β2x2++βrxr+e

其中 β0,β1,,βr 是需要估计的回归系数。如果数据集由 n(xi,Yi) 组成,其中 i=1,,n,那么 β0,,βr 的最小二乘估计量——我们称之为 B0,,Br——就是使下式最小化的值:

i=1n(YiB0B1xiB2xi2Brxir)2

为了确定这些估计量,我们分别求解 B0Br 的偏导数,并令这些偏导数等于 0。在执行这些操作,并重新排列所得的方程之后,我们得到的最小二乘估计 B0,B1,,Br 满足以下由 r+1 个线性方程组(称为 正规方程normal equations))成的方程组:

i=1nYi=B0n+B1i=1nxi+B2i=1nxi2++Bri=1nxiri=1nxiYi=B0i=1nxi+B1i=1nxi2+B2i=1nxi3++Bri=1nxir+1i=1nxi2Yi=B0i=1nxi2+B1i=1nxi3+B2i=1nxi4++Bri=1nxir+2i=1nxirYi=B0i=1nxir+B1i=1nxir+1++Bri=1nxi2r

在对一组数据进行多项式拟合时,通常可以通过研究散点图来确定多项式的必要阶数。我们强调,应始终使用最小可能阶数来充分描述数据的关系。例如,尽管通常可以找到一个阶数为 n 的多项式以通过所有 n(xi,Yi),但我们很难对这种拟合寄予太大的信心。

与线性回归相比,使用多项式拟合来预测那些远离拟合数据(i=1,,n)的 x0 的响应值是一件更为危险的事情。多项式拟合可能只在 xii=1,,n)周围的某一区域有效,并且这个区域并不包括 x0

习题 9.6 为如下的数据进行多项式拟合。

x Y
1 20.6
2 30.8
3 55.0
4 71.4
5 97.3
6 131.8
7 156.3
8 197.3
9 238.7
10 291.7

解 9.6.

代码
x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(20.6, 30.8, 55.0, 71.4, 97.3, 131.8, 156.3, 197.3, 238.7, 291.7)
plot(x, y)
图 9.10

数据的散点图()暗示着 Yx 之间存在二次关系:

Y=β0+β1x+β2x2+e

ixi=55,ixi2=385,ixi3=3025,ixi4=25,333iYi=1291.1,ixiYi=9549.3,ixi2Yi=77,758.9

所以,最小二乘估计就是如下的方程组的解:

(9.20)1291.1=10B0+55B1+385B29549.3=55B0+385B1+3025B277758.9=385B0+3025B1+25333B2

解如上的方程组可以得到最小二乘估计为:

B0=12.59326,B1=6.326172,B2=2.122818

因此,估计的二次回归方程为:

Y=12.59+6.33x+2.12x2

数据的散点图和拟合的二次回归方程图如 所示。

代码
x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(20.6, 30.8, 55.0, 71.4, 97.3, 131.8, 156.3, 197.3, 238.7, 291.7)
plot(x, y)

beta0 <- 12.59
beta1 <- 6.33
beta2 <- 2.12
curve(beta0 + beta1 * x + beta2 * x^2, from=0, to=10, col="red", lwd=2, add=TRUE)
图 9.11

注记

可以使用 R 的 lm() 实现 的多项式拟合。

x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(20.6, 30.8, 55.0, 71.4, 97.3, 131.8, 156.3, 197.3, 238.7, 291.7)
model <- lm(y ~ x + I(x^2))

summary(model)

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5482 -2.4869 -0.4486  2.7006  4.8739 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.6433     4.3477   2.908   0.0227 *  
x             6.2971     1.8158   3.468   0.0104 *  
I(x^2)        2.1250     0.1609  13.209 3.33e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.697 on 7 degrees of freedom
Multiple R-squared:  0.9987,    Adjusted R-squared:  0.9984 
F-statistic:  2745 on 2 and 7 DF,  p-value: 7.368e-11
plot(x, y)
curve(12.6433 + 6.2971 * x + 2.1250 * x^2, from=0, to=10, col="blue", lwd=2, add=TRUE)

备注

在矩阵表示法中, 可以表示为:

[1291.19549.377758.9]=[1055385553853025385302525,333][B0B1B2]

其解为:

[B0B1B2]=[1055385553853025385302525,333]1[1291.19549.377,758.9]

9.10 Multiple linear regression

9.11 用于二元输出数据的逻辑回归模型

在本节中,我们考虑实验结果要么是成功要么是失败的实验。我们假设可以在不同的水平(x)执行这些实验,并且在水平 x 处实验成功的概率为 p(x),其中 <x<。如果 p(x) 具有以下形式:

p(x)=ea+bx1+ea+bx

那么,我们称这些实验来自于一个 逻辑回归模型logistic regression model),而称 p(x)逻辑回归函数logistics regression function)。

  • 如果 b>0,那么 p(x)=1[e(a+bx)+1] 是一个递增函数,随着 x 而收敛到 1。
  • 如果 b<0,那么 p(x) 是一个递减函数,随着 x 而收敛到 0。
  • b=0 时,p(x) 是一个常数。

给出了 逻辑回归函数 的图形,从图中可以看出,这些曲线呈现 S 形状。

图 9.12: 逻辑回归函数

p(x) 写为 p(x)=1[11+ea+bx] 并对 x 求导,得到:

xp(x)=bea+bx(1+ea+bx)2=bp(x)[1p(x)]

因此,p(x) 的变化率取决于 x,并且在 p(x) 接近 0.5 时的那些 x 处,p(x) 的变化率最大。例如,在 p(x)=0.5 时,xp(x)=0.25b,而在 p(x)=0.8 时,x 处,xp(x)=0.16b

如果我们将 o(x) 定义为实验在 x 处成功时的 赔率),那么:

o(x)=p(x)1p(x)=ea+bx

因此,当 b>0 时,赔率 随着 x 的增加而指数级增加;当 b<0 时,赔率 随着 x 的增加而指数级减少。取 赔率 取对数运算,我们有:

log(o(x))=a+bx

这表明,对数赔率(log odds)——即 logit 函数——是一个线性函数:

log(o(x))=a+bx

逻辑回归函数的参数 ab 是未知的,需要对其进行估计。可以通过最大似然法来估计逻辑回归函数的参数。也就是说,假设在 x1,,xk 进行实验,令 Yi(成功时为 1,失败时为 0) 表示在 xi 处进行实验时的实验结果。然后,使用伯努利概率密度函数(即单次试验的概率密度函数为二项概率密度函数):

P(Yi=yi)=[p(xi)]yi[1p(xi)]1yi=(ea+bxi1+ea+bxi)yi(11+ea+bxi)1yi,yi=0,1

因此,在 xi 处实验结果为 yi 的概率(对于所有 i=1,,k) 为:

P(Yi=yi,i=1,,k)=i(ea+bxi1+ea+bxi)yi(11+ea+bxi)1yi=i(ea+bxi)yi1+ea+bxi

等式两边取对数得到:

log(P{Yi=yi,i=1,,k})=i=1kyi(a+bxi)i=1klog(1+ea+bxi)

现在可以通过数值方法求解使得上述似然函数取得最大值的 ab 的值,从而获得最大似然估计。然而,由于似然函数是非线性的,所以通常需要迭代的方法来计算最大似然估计;因此,通常需要使用专门的软件来获得逻辑回归函数的参数估计。

glm()

glm() 函数用于拟合 广义线性模型Generalized Linear Models, GLM)。glm() 是一个非常强大的函数,可以处理各种类型的回归模型,包括线性回归、逻辑回归、泊松回归等。

可以通过 family 参数来控制 glm() 的回归类型,例如:

  • family = gaussian(默认)用于线性回归
  • family = binomial 用于逻辑回归
  • family = poisson 用于泊松回归

可以使用 R 确定逻辑回归模型中 ab 的估计。可以使用 R 中的 glm(y ~ x, family = "binomial") 确定逻辑回归模型。例如,假设某个实验进行了 7 次,在 x=(5,9,13,22,22,24,30) 时,其对应的结果为 (0,1,0,1,0,1,1),其中 1 表示该实验成功,0 表示实验失败。则估计的逻辑回归方程如下:

x <- c(5, 9, 3, 12, 22, 24, 30)
y <- c(0, 1, 0, 1, 0, 1, 1)
model <- glm(y ~ x, family = "binomial")
model

Call:  glm(formula = y ~ x, family = "binomial")

Coefficients:
(Intercept)            x  
    -1.2678       0.1103  

Degrees of Freedom: 6 Total (i.e. Null);  5 Residual
Null Deviance:      9.561 
Residual Deviance: 8.025    AIC: 12.03
注意

在原书中,此处使用的是 glm(y ~ x),如果没有指定 family 参数,则 glm() 实际执行的是线性回归,而不是逻辑回归。这里要特别注意。

也就是说,在逻辑回归模型 p(x)=ea+bx1+ea+bx 中,ab 的估计值分别是:-1.2678 和 0.1103。

当响应数据为二元数据时,虽然逻辑回归模型是最常用的模型,但通常也会使用其他模型。例如,在输入水平为 x 时,假设 p(x)(实验成功的概率)是一个随 x 的增加而增加的函数是合理的。在这种情况下,常常假设 p(x) 具有特定的概率分布函数。实际上,当 b>0 时,逻辑回归模型就满足这种形式。此时,因为 p(x) 就是参数为 μ=a/bν=1/b 的逻辑分布随机变量的分布函数(?sec-5_9)。

可用于二元响应数据的另一种模型是 probit model,在这种模型中,假设对于某些常数 α,β>0

p(x)=Φ(α+βx)=12πα+βxey2/2dy

换句话说,p(x) 等于标准正态分布随机变量小于 α+βx 的概率。

例 9.10 常见的假设是:当一个动物暴漏在剂量水平为 x 下的某化学物质时,动物是否会生病。该假设使用了阈值模型(threshold model),阈值模型认为每个动物都有一个随机阈值,当剂量水平超过该阈值时,动物就会生病。有时,可以用指数分布来作为阈值分布。例如,Freedman 和 Zeisel 在论文(From Mouse to Man: The Quantitative Assessment of Cancer Risks,” Statistical Science, 1988, 3, 1, 3–56)中的模型假设:一只暴露于 x 剂量(以 ppm 计量)的 DDT(二氯二苯基三氯乙烷,一种合成的有机氯杀虫剂,最早于 1940 年代被引入,用于控制蚊子、苍蝇等传播疾病的昆虫,并在农业中广泛用于保护农作物免受害虫侵害) 下的老鼠将会以如下的概率患上肝癌:

p(x)=1eαx,x>0

由于指数分布的“无记忆性”,这相当于假设:如果老鼠在接收到部分剂量 x 后仍然健康,那么剂量水平 x 带来的风险与之前没有接收到任何剂量时一样。

根据 Freedman 和 Zeisel 的研究,暴露于 250 ppm DDT 下的 111 只老鼠中,有 84 只发展成癌症。因此,可以通过以下公式估计 α

1e250α^=84111

即:

α^=log(27/111)250=0.005655

Problems