3.2(上) 线性回归模型
这是一篇有关《统计学习基础》,原书名The Elements of Statistical Learning的学习笔记,该书学习难度较高,有很棒的学者将其翻译成中文并放在自己的个人网站上,翻译质量非常高,本博客中有关翻译的内容都是出自该学者的网页,个人解读部分才是自己经过查阅资料和其他学者的学习笔记,结合个人理解总结成的原创内容。
有关ESL更多的学习笔记的markdown文件,可在作者GitHub上查看下载。
原文 | The Elements of Statistical Learning |
---|---|
翻译 | szcf-weiya |
时间 | 2018-08-21 |
解读 | Hytn Chen |
更新 | 2020-02-12 |
翻译原文
!!! note “更新笔记”
@2017-12-27 最近在看 power analysis 的有关内容,R in Action 一书中举了例子来讨论对于不同模型如何进行 power analysis,但其中关于线性模型的部分没有想清楚.后来在这一节的式 (3.13) 找到答案,为丰富本节内容,遂更新与式 3.13 有关的 power analysis.
正如第二章介绍的那样,我们有输入向量 X T = ( X 1 , X 2 , … , X p ) X^T=(X_1,X_2,\ldots,X_p) XT=(X1,X2,…,Xp),而且想要预测实值输出变量 Y Y Y.线性回归模型有如下形式
f ( X ) = β 0 + ∑ j = 1 p X j β j (3.1) f(X)=\beta_0+\sum\limits_{j=1}^pX_j\beta_j \tag{3.1} f(X)=β0+j=1∑pXjβj(3.1)
线性模型要么假设回归函数 E ( Y ∣ X ) \rm{E}(Y\mid X) E(Y∣X) 是线性的,要么假设它本身是一个合理的近似.这里 β j \beta_j βj 是未知的参数或系数,变量 X j X_j Xj 可以有下列不同的来源:
- 定量的输入变量
- 定量输入变量的变换,比如对数变换,平方根变换或者平方变换
- 基函数展开,比如 X 2 = X 1 2 , X 3 = X 1 3 X_2=X_1^2,X_3=X_1^3 X2=X12,X3=X13,产生一个多项式的表示
- 定性输入变量水平 (level) 的数值或“虚拟”编码.举个例子,如果 G G G 是有 5 个水平的输入变量,我们可能构造 X j , j = 1 , … , 5 X_j,j=1,\ldots,5 Xj,j=1,…,5 使得 X j = I ( G = j ) X_j=I(G=j) Xj=I(G=j).通过一系列独立于水平的常数,整个 X j X_j Xj 可以用来表示 G G G 的效应,因为在 ∑ j = 1 5 X j β j \sum_{j=1}^5X_j\beta_j ∑j=15Xjβj 中,其中一个 X j X_j Xj 的系数为 1,其它的都是 0.
- 变量的交叉影响,举个例子, X 3 = X 1 ⋅ X 2 X_3=X_1\cdot X_2 X3=X1⋅X2.
无论 X j X_j Xj 是哪个来源,模型关于参数都是线性的.
一般地,我们有一个用来估计参数
β
\beta
β 的训练数据集合
(
x
1
,
y
1
)
,
…
,
(
x
N
,
y
N
)
(x_1,y_1),\ldots,(x_N,y_N)
(x1,y1),…,(xN,yN).每个
x
i
=
(
x
i
1
,
x
i
2
,
…
,
x
i
p
)
T
x_i=(x_{i1},x_{i2},\ldots,x_{ip})^T
xi=(xi1,xi2,…,xip)T 是第
i
i
i 种情形的特征测量值的向量.最受欢迎的估计方法是 最小二乘 (least squares),我们取参数
β
=
(
β
0
,
β
1
,
…
,
β
p
)
T
\beta=(\beta_0,\beta_1,\ldots,\beta_p)^T
β=(β0,β1,…,βp)T 使得下式的残差平方和最小
R
S
S
(
β
)
=
∑
i
=
1
N
(
y
i
−
f
(
x
i
)
)
2
=
∑
i
=
1
N
(
y
i
−
β
0
−
∑
j
=
1
p
x
i
j
β
j
)
2
(3.2)
\begin{aligned} \rm{RSS}(\beta)&=\sum\limits_{i=1}^N(y_i-f(x_i))^2\\ &=\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2\tag{3.2} \end{aligned}
RSS(β)=i=1∑N(yi−f(xi))2=i=1∑N(yi−β0−j=1∑pxijβj)2(3.2)
从统计学的观点来看,如果训练观测值 ( x i , y i ) (x_i,y_i) (xi,yi) 是从总体独立随机抽取的,则该标准是很合理的.即使 x i x_i xi's 不是随机选取的,如果在给定输入 x i x_i xi 的条件下 y i y_i yi 条件独立,这个准则也是有效的.图 3.1 展示了在充满数据对 ( X , Y ) (X,Y) (X,Y) 的 I R p + 1 \rm{IR}^{p+1} IRp+1 维空间的最小二乘拟合的几何意义.注意到 ( 3.2 ) (3.2) (3.2) 对模型 (3.1) 的有效性没有作假设,而是简单地找到对数据最好的线性拟合.无论数据是怎样产生的,最小二乘拟合直观上看是满意的,这个准则衡量了平均的拟合误差.
图 3.1 X ∈ R 2 X\in R^2 X∈R2 中的最小二乘拟合.我们寻找关于 X X X 并且使 Y Y Y 的残差最小的线性函数.
那我们应该怎样最小化
(
3.2
)
(3.2)
(3.2)?记
X
\mathbf{X}
X 为
N
×
(
p
+
1
)
N\times (p+1)
N×(p+1) 的矩阵,矩阵每一行为一个输入向量(在第一个位置为 1),类似地令
y
\mathbf{y}
y 为训练集里的
N
N
N 维输出向量.然后我们可以将残差平方和写成如下形式
R
S
S
(
β
)
=
(
y
−
X
β
)
T
(
y
−
X
β
)
(3.3)
\rm{RSS}(\beta)=(\mathbf{y}-\mathbf{X}\beta)^T(\mathbf{y}-\mathbf{X}\beta)\tag{3.3}
RSS(β)=(y−Xβ)T(y−Xβ)(3.3)
这是含 p + 1 p+1 p+1 个参数的二次函数.关于 β \beta β 求导有
∂ R S S ∂ β = − 2 X T ( y − X β ) ∂ 2 R S S ∂ β ∂ β T = 2 X T X (3.4) \begin{array}{ll} \dfrac{\partial \rm{RSS}}{\partial \beta}&=-2\mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta)\\ \dfrac{\partial^2 \rm{RSS}}{\partial \beta\partial \beta^T}&=2\mathbf{X}^T\mathbf{X} \end{array} \tag{3.4} ∂β∂RSS∂β∂βT∂2RSS=−2XT(y−Xβ)=2XTX(3.4)
假设
X
\mathbf{X}
X 列满秩,因此
X
T
X
\mathbf{X}^T\mathbf{X}
XTX 是正定的,我们令一阶微分为 0,即
X
T
(
y
−
X
β
)
=
0
(3.5)
\mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta)=0\tag{3.5}
XT(y−Xβ)=0(3.5)
得到唯一解
β
^
=
(
X
T
X
)
−
1
X
T
y
(3.6)
\hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.6}
β^=(XTX)−1XTy(3.6)
在输入向量
x
0
x_0
x0 下的预测值由
f
^
(
x
0
)
=
(
1
:
x
0
)
T
β
^
\hat{f}(x_0)=(1:x_0)^T\hat{\beta}
f^(x0)=(1:x0)Tβ^ 给出;在训练输入值下的拟合值为
y
^
=
X
β
^
=
X
(
X
T
X
)
−
1
X
T
y
(3.7)
\hat{y}=\mathbf{X}\hat{\beta}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.7}
y^=Xβ^=X(XTX)−1XTy(3.7)
其中, y ^ i = f ^ ( x i ) \hat y_i=\hat f(x_i) y^i=f^(xi).在式 ( 3.7 ) (3.7) (3.7) 中出现的矩阵 H = X ( X T X ) − 1 X T \mathbf{H=X(X^TX)^{-1}X^T} H=X(XTX)−1XT 有时称之为“帽子”矩阵,因为其为 y \mathbf y y 戴了一顶帽子.
【Hytn注:这里的 β ^ \hat \beta β^其实就是由多个训练数据所拟合得出的参数了,当给定多个训练点且满足正定条件的时候, β ^ \hat \beta β^可以唯一确定,所以就是根据拟合得到的参数 β ^ \hat \beta β^,结合训练数据的输入特征 X \rm{X} X,得出其对应的估计输出向量 y ^ \hat y y^。】
图 3.2 含两个预测的最小二乘回归的 N N N 维几何图形.输出向量 y \mathbf{y} y 正交投影在输入向量 x 1 \mathbf{x}_1 x1 和 x 2 \mathbf{x}_2 x2 张成的超平面中.投影向量 y ^ \hat{\mathbf{y}} y^ 表示最小二乘法的预测值.
图 3.2 展示了在 I R N \rm{IR}^N IRN 中最小二乘估计的另一种几何表示.我们记 X \mathbf{X} X 的列向量为 x 0 , x 1 , … , x p \rm{x}_0,\rm{x}_1,\ldots,\rm{x}_p x0,x1,…,xp,其中 x 0 ≡ 1 \rm{x}_0 \equiv \rm{1} x0≡1.下文中第一列都是这样假设.这些向量张成了 I R N \rm{IR}^N IRN 的子空间,也被称作 X \mathbf{X} X 的列空间.我们通过选择 β ^ \hat{\beta} β^ 来最小化 R S S ( β ) = ∥ ( y − X β ) ∥ 2 \rm{RSS}(\beta)=\lVert(\mathbf{y}-\mathbf{X}\beta)\rVert^2 RSS(β)=∥(y−Xβ)∥2,则残差向量 y − y ^ \mathbf{y-\hat{y}} y−y^ 正交于子空间. ( 3.5 ) (3.5) (3.5)式描述了这种正交,然后估计值 y ^ \hat{\rm{y}} y^ 因此可以看成是 y \mathbf{y} y 在子空间中的正交投影.帽子矩阵 H \mathbf{H} H 计算了正交投影,因此也被称作投影矩阵.
可能会出现 X \mathbf{X} X 的列不是线性独立的,则 X \mathbf{X} X 不是满秩的.举个例子,如果两个输入是完全相关的,(比如, x 2 = 3 x 1 x_2=3x_1 x2=3x1).则矩阵 X T X \mathbf{X^TX} XTX 是奇异的,并且最小二乘的系数 β ^ \hat{\beta} β^ 不是唯一的.然而,拟合值 y ^ = X β ^ \mathbf{\hat{y}}=\mathbf{X}\hat{\beta} y^=Xβ^ 仍然是 y \mathbf{y} y 在列空间 X \mathbf{X} X 的投影;用 X \mathbf{X} X 的列向量来表示这种投射的方式不止一种.当一个或多个定性输入用一种冗余的方式编码时经常出现非满秩的情形.通过重编码或去除 X \mathbf{X} X 中冗余的列等方式可以解决非唯一表示的问题.大多数回归软件包会监测这些冗余并且自动采用一些策略去除冗余项.信号和图像分析中经常发生秩缺失,其中输入的个数 p p p 可以大于训练的情形个数 N N N.在这种情形下,特征经常通过滤波来降维或者对拟合进行正则化.(5.2.3节和第 18 章)
截至目前我们已经对数据的真实分布做的假设非常少.为了约束 β ^ \hat{\beta} β^ 的取样特点,我们现在假设观测值 y i y_i yi 不相关,且有固定的方差 σ 2 \sigma^2 σ2,并且 x i x_i xi 是固定的(非随机).最小二乘法的参数估计的方差-协方差阵可以很容易从式 ( 3.6 ) (3.6) (3.6) 导出:
V a r ( β ^ ) = ( X T X ) − 1 σ 2 (3.8) \rm{Var}(\hat{\beta})=(\mathbf{X}^T\mathbf{X})^{-1}\sigma^2\tag{3.8} Var(β^)=(XTX)−1σ2(3.8)
一般通过下式来估计方差
σ
2
\sigma^2
σ2
σ
^
2
=
1
N
−
p
−
1
∑
i
=
1
N
(
y
i
−
y
^
i
)
2
\hat{\sigma}^2=\frac{1}{N-p-1}\sum\limits_{i=1}^N(y_i-\hat{y}_i)^2
σ^2=N−p−11i=1∑N(yi−y^i)2
之所以分母是 N − p − 1 N-p-1 N−p−1 而不是 N N N,是因为 N − p − 1 N-p-1 N−p−1 使得 σ ^ 2 \hat{\sigma}^2 σ^2 为无偏估计: E ( σ ^ 2 ) = σ 2 E(\hat{\sigma}^2)=\sigma^2 E(σ^2)=σ2
为了对参数和模型进行推断,需要一些额外的假设.我们现在假设式
(
3.1
)
(3.1)
(3.1) 是关于均值的正确模型;则
Y
Y
Y 的条件期望关于
X
1
,
X
2
,
…
,
X
p
X_1,X_2,\ldots,X_p
X1,X2,…,Xp 是线性的.我们也假设
Y
Y
Y 与其期望的偏差是可加的且是正态的.因此
Y
=
E
(
Y
∣
X
1
,
…
,
X
p
)
+
ε
=
β
0
+
∑
j
=
1
p
X
j
β
j
+
ε
(3.9)
\begin{aligned} Y&=\rm{E}(Y\mid X_1,\ldots,X_p)+\varepsilon\\ &=\beta_0+\sum\limits_{j=1}^pX_j\beta_j+\varepsilon\tag{3.9} \end{aligned}
Y=E(Y∣X1,…,Xp)+ε=β0+j=1∑pXjβj+ε(3.9)
其中误差 ϵ \epsilon ϵ 是期望值为 0 方差为 σ 2 \sigma^2 σ2 的高斯随机变量,记作 ε ∼ N ( 0 , σ 2 ) \varepsilon\sim N(0,\sigma^2) ε∼N(0,σ2)
由式
(
3.9
)
(3.9)
(3.9),可以很简单地证明
β
^
∼
N
(
β
,
(
X
T
X
)
−
1
σ
2
)
(3.10)
\hat{\beta}\sim N(\beta,(\mathbf{X}^T\mathbf{X})^{-1}\sigma^2) \tag{3.10}
β^∼N(β,(XTX)−1σ2)(3.10)
个人解读
对于参数 β ^ \hat \beta β^的最终证明究竟如何一步步得来?
首先式(3.8)证明过程如下:
β
^
=
(
X
T
X
)
−
1
X
T
(
X
β
+
ϵ
)
=
(
X
T
X
)
−
1
X
T
X
β
+
(
X
T
X
)
−
1
X
T
ϵ
=
β
+
(
X
T
X
)
−
1
X
T
ϵ
\begin{aligned}\hat{\beta} &=\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T}(X \beta+\epsilon) \\&=\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \mathbf{X} \beta+\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \epsilon \\&=\beta+\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \epsilon\end{aligned}
β^=(XTX)−1XT(Xβ+ϵ)=(XTX)−1XTXβ+(XTX)−1XTϵ=β+(XTX)−1XTϵ
由于原文中声明误差误差
ϵ
\epsilon
ϵ 是期望值为 0 方差为
σ
2
\sigma^2
σ2 的高斯随机变量,记作
ε
∼
N
(
0
,
σ
2
)
\varepsilon\sim N(0,\sigma^2)
ε∼N(0,σ2)所以
E
[
β
^
]
=
β
+
(
X
T
X
)
−
1
X
T
0
=
β
\begin{aligned} E[\hat{\beta}] &=\beta+\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \mathbf{0} \\ &=\beta \end{aligned}
E[β^]=β+(XTX)−1XT0=β
β ^ ∼ N ( β , ( X T X ) − 1 X T σ 2 ) \hat{\beta} \sim N\left(\beta,\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \sigma^{2}\right) β^∼N(β,(XTX)−1XTσ2)
而对于方差的证明有两条路,第一条如下:
Var
(
β
^
)
=
(
X
T
X
)
−
1
X
T
X
(
X
T
X
)
−
1
σ
2
=
(
X
T
X
)
−
1
σ
2
⇒
β
^
∼
N
(
β
,
(
X
T
X
)
−
1
σ
2
)
\begin{aligned} \operatorname{Var}(\hat{\beta}) &=\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \mathbf{X}^{T} \mathbf{X}\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \sigma^{2} \\ &=\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \sigma^{2} \\ \Rightarrow & \hat{\beta} \sim N\left(\beta,\left(\mathbf{X}^{T} \mathbf{X}\right)^{-1} \sigma^{2}\right) \end{aligned}
Var(β^)⇒=(XTX)−1XTX(XTX)−1σ2=(XTX)−1σ2β^∼N(β,(XTX)−1σ2)
第二条根据公式:
V a r ( β ^ ) = E ( β ^ 2 ) − [ E ( β ^ ) ] 2 \rm{Var}(\hat \beta) = \rm{E}({\hat \beta}^2) - [\rm{E}(\hat \beta)]^2 Var(β^)=E(β^2)−[E(β^)]2
其中 β ^ 2 \hat \beta^2 β^2展开式如下:
β ^ 2 = β 2 + 2 ( X T X ) − 1 X T ϵ β + ( X T X ) − 1 ϵ 2 \hat \beta^2 = \beta^2 + 2(X^TX)^{-1}X^T\epsilon\beta + (X^TX)^{-1}\epsilon^2 β^2=β2+2(XTX)−1XTϵβ+(XTX)−1ϵ2
其中第二项的期望为0,第三项的期望是 ( X T X ) − 1 σ 2 (X^TX)^{-1}\sigma^2 (XTX)−1σ2,将结果整合进行消项之后可得最终结果,与原文结论一致。
卡方分布,n个服从标准正态分布的随机变量的平方和构成新的随机变量,这个随机变量的分布规律称为卡方分布。下文就使用卡方分布来检验统计量,比较期望结果和实际结果之间的差别,学者Wang CS在其博客中对卡方分布进行了具体描述,表明了这个检验统计量提供了一种期望值与观察值之间差异的度量方法,最后反映在 χ 2 \chi^2 χ2的数值大小上。
那么对于这个数值,究竟大到什么程度差异才算显著?这里摘录博客中的描述:
需要根据自由度,设定的显著性水平查找卡方分布表判定,就是根据三要素:一个公式,一张分布表,一张概率密度图。
具体的卡方分布计算实例在博客中也有例举,更加详细的实例和讲解在这篇名为t分布, 卡方x分布,F分布的博客中也有例举,在此不进行展开。
继续原文
这是一个有上述均值向量和方差-协方差矩阵的多元正态分布.同时有
(
N
−
p
−
1
)
σ
^
2
∼
σ
2
χ
N
−
p
−
1
2
(3.11)
(N-p-1)\hat \sigma^2\sim \sigma^2\chi^2_{N-p-1}\tag{3.11}
(N−p−1)σ^2∼σ2χN−p−12(3.11)
是一个自由度为 N − p − 1 N-p-1 N−p−1 的卡方分布.另外, β ^ \hat{\beta} β^ 和 σ ^ 2 \hat\sigma^2 σ^2 是统计独立的.我们利用这些分布性质得到假设检验以及对于参数 β j \beta_j βj 的置信区间
图 3.3 t 30 t_{30} t30, t 100 t_{100} t100以及标准正态三种分布的尾概率 Pr ( ∣ Z ∣ > z ) \Pr(|Z|>z) Pr(∣Z∣>z).图中显示了在 p = 0.05 p=0.05 p=0.05 和 0.01 0.01 0.01 水平下检验显著性的合适分位数.当 N N N 大于 100 时 t t t 分布与标准正态分布之间的差异很小.
为了检验系数
β
j
\beta_j
βj 的这一假设,我们构造标准化因数或者 Z-分数.
z
j
=
β
^
j
σ
^
v
j
(3.12)
z_j=\dfrac{\hat{\beta}_j}{\hat{\sigma}\sqrt{v_j}}\tag{3.12}
zj=σ^vjβ^j(3.12)
其中 v j v_j vj 是 ( X T X ) − 1 (\mathbf{X}^T\mathbf{X})^{-1} (XTX)−1 的第 j j j 个对角元.零假设为 β j = 0 \beta_j=0 βj=0, z j z_j zj 分布为 t N − p − 1 t_{N-p-1} tN−p−1(自由度为 N − p − 1 N-p-1 N−p−1 的 t t t 分布),因此当 z j z_j zj 的绝对值较大时会拒绝零假设.如果用已知的 σ \sigma σ 值替换 σ ^ \hat{\sigma} σ^,则 z j z_j zj 服从标准正态分布. t t t 分布和标准正态分布在尾概率之间的差异随着样本规模增大可以忽略,因此我们一般使用正态的分位数(图 3.3).
(Hytn注:卡方分布自由度的计算分为两种情况,第一种:拟合优度检验。就是如果式子中包含n个独立的随机变量,和由它们所构成的k个样本统计量,则该表达式的自由度为n-k。比如统计量是n个随机变量的平均数,统计量只有一个,因此自由度为n-1。第二种:两个变量的独立性检验。若列联表为h行k列,则自由度=(h-1) x (k-1))
我们经常需要同时检验多个系数 (groups of coefficients) 的显著性.举个例子,检验有 k k k 个水平的分类变量是否要从模型中排除,我们需要检验用来表示水平的虚拟变量的系数是否可以全部设为 0.这里我们利用 F F F 统计量
F = ( R S S 0 − R S S 1 ) / ( p 1 − p 0 ) R S S 1 / ( N − p 1 − 1 ) (3.13) F=\dfrac{(\rm{RSS}_0-\rm{RSS}_1)/(p_1-p_0)}{\rm{RSS}_1/(N-p_1-1)}\tag{3.13} F=RSS1/(N−p1−1)(RSS0−RSS1)/(p1−p0)(3.13)
其中 R S S 1 \rm{RSS}_1 RSS1 是有 p 1 + 1 p_1+1 p1+1 个参数的大模型的最小二乘法拟合的残差平方和, R S S 0 \rm{RSS}_0 RSS0 是有 p 0 + 1 p_0+1 p0+1 参数的小模型的最小二乘法拟合的残差平方和,其有 p 1 − p 0 p_1-p_0 p1−p0 个参数约束为 0. F F F 统计量衡量了在大模型中每个增加的系数对残差平方和的改变,而且用 σ 2 \sigma^2 σ2 的估计值进行标准化.在高斯分布以及小模型的零假设为正确的情况下, F F F 统计量服从 F p 1 − p 0 , N − p 1 − 1 F_{p_1-p_0,N-p_1-1} Fp1−p0,N−p1−1 分布.可以证明 ( 3.12 ) (3.12) (3.12) 式中的 z j z_j zj 等于从模型中去除单系数 β j \beta_j βj 的 F F F 统计量(习题 3.1),当 N N N 足够大时, F p 1 − p 0 , N − p 1 − 1 F_{p_1-p_0,N-p_1-1} Fp1−p0,N−p1−1 近似 χ p 1 − p 0 2 \chi^2_{p_1-p_0} χp1−p02.
!!! info “weiya 注:Ex. 3.1”
已解决,详见 Issue 69: Ex. 3.1.
!!! note “weiya 注”
在 R 语言中,可以采用 pwr
包来进行 power analysis.对于线性模型(比如多重回归),pwr.f2.test(u=, v=, f2=, sig.level=, power=)
可以用来进行 power analysis,其中 u
和 v
分别是分子和分母的自由度,f2
是 有效大小 (effect size).
当需要衡量一个变量集对输出变量的影响,f2
计算公式为,
f
2
=
R
2
1
−
R
2
,
f^2 = \frac{R^2}{1-R^2}\,,
f2=1−R2R2,
当需要衡量一个变量集优于另一个变量集的影响时,
f
2
=
R
A
B
2
−
R
A
2
1
−
R
A
B
2
,
f^2 = \frac{R^2_{AB}-R^2_A}{1-R_{AB}^2}\,,
f2=1−RAB2RAB2−RA2,
其中
R
A
2
R_A^2
RA2 是变量集
A
A
A 解释的方差,
R
A
B
2
R_{AB}^2
RAB2 是变量集
A
A
A 和
B
B
B 共同解释的方差.
这里举一个实际例子来说明如何在线性回归模型中进行 power analysis.
假设我们要研究老板的领导风格对公司员工满意度的影响,同时员工满意度还与薪酬有关,其中薪酬用 3 个变量来衡量,而老板领导风格用 4 个变量来衡量.根据以往经验,薪酬能够解释员工满意度方差的 30%,现在问题是确定多少的调查个体,使得显著度为 .05,而且 power 为 .90.求解代码为`pwr.f2.test(u=3, f2=(.35-.30)/(1-.35), sig.level=0.05, power=0.90)`,输出`v=184.2426`.
其中`u`和`v`分别被称为分子、分母自由度.结合 \eqref{3.13} 式,我们可以看出 $p_1-p_0=u, N-p_1-1=v$,并且这里 $p_0=4, p_1=7$.换个角度看,其实式 \eqref{3.13} 就是用于 power analysis 的检验统计量,原假设为可以将领导力水平的四个变量系数设为 0,若 $F$ 足够大,则会拒绝原假设.
类似地,我们可以分离式 ( 3.10 ) (3.10) (3.10) 中的 β j \beta_j βj 得到 β j \beta_j βj 的置信水平为 1 − 2 α 1-2\alpha 1−2α 的置信区间
( β ^ j − z 1 − α v j 1 / 2 σ ^ , β ^ j + z 1 − α v j 1 / 2 σ ^ ) (3.14) (\hat{\beta}_j-z^{1-\alpha}v_j^{1/2}\hat{\sigma},\hat{\beta}_j+z^{1-\alpha}v_j^{1/2}\hat{\sigma})\tag{3.14} (β^j−z1−αvj1/2σ^,β^j+z1−αvj1/2σ^)(3.14)
这里,
z
1
−
α
z^{1-\alpha}
z1−α 是正态分布的
1
−
α
1-\alpha
1−α 分位数.
z
1
−
0.025
=
1.96
,
z
1
−
0.05
=
1.645
,
e
t
c
.
\begin{aligned} z^{1-0.025}&=1.96,\\ z^{1-0.05}&=1.645,etc. \end{aligned}
z1−0.025z1−0.05=1.96,=1.645,etc.
因此报告 β ^ ± 2 ⋅ s e ( β ^ ) \hat{\beta}\pm 2\cdot se(\hat{\beta}) β^±2⋅se(β^) 的标准做法是约 95 % 95\% 95% 的置信区间.即便没有高斯误差的假设,区间也近似正确,因为当样本规模 N ∼ ∞ N\sim \infty N∼∞ 时,覆盖范围近似为 1 − 2 α 1-2\alpha 1−2α.
运用类似的方法我们可以得到全体系数向量 β \beta β 的近似置信集,记作
C β = { β ∣ ( ( β ^ − β ) T X T X ( β ^ − β ) ≤ σ ^ 2 χ p + 1 2 ( 1 − α ) } (3.15) C_{\beta} = \{\beta \mid ((\hat{\beta}-\beta)^T\mathbf{X^TX}(\hat{\beta}-\beta)\le \hat{\sigma}^2\chi^{2^{(1-\alpha)}}_{p+1} \}\tag{3.15} Cβ={β∣((β^−β)TXTX(β^−β)≤σ^2χp+12(1−α)}(3.15)
其中, χ ℓ 2 ( 1 − α ) \chi_\ell^{2^{(1-\alpha)}} χℓ2(1−α) 是 ℓ \ell ℓ 个自由度的卡方分布的 1 − α 1-\alpha 1−α 的分位数:举个例子, χ 2 ( 1 − 0.05 ) _ 5 = 11.1 , χ 2 ( 1 − 0.1 ) _ 5 = 9.2 \chi^{2^{(1-0.05)}}\_5=11.1,\chi^{2^{(1-0.1)}}\_5=9.2 χ2(1−0.05)_5=11.1,χ2(1−0.1)_5=9.2.这个关于 β \beta β 的置信集导出关于真函数 f ( x ) = x T β f(x)=x^T\beta f(x)=xTβ 对应的置信集,记作 x T β ∣ β ∈ C β \\{x^T\beta\mid \beta \in C_\beta\\} xTβ∣β∈Cβ(练习 3.2;并见 5.2.2 节中函数置信带的例子图 5.4)
!!! info “weiya 注:Ex. 3.2”
已解决,详见 Issue 120: Ex. 3.2.
例子:前列腺癌
这个例子的数据来自 Stamey et al. (1989)[^1].他们检验即将接受彻底前列腺切除术的男性的前列腺特定抗原水平和一系列临床措施之间的相关性.变量有癌体积的对数值 (lcavol)、前列腺重量的对数值 (lweight)、良性前列腺增生数量 (lbph)、精囊浸润 (svi)、包膜浸透的对数值 (lcp)、Gleason 得分 (gleason)、Gleason 得分为 4 或 5 的比例 (pgg45).表 3.1 中给出的预测变量的相关性矩阵显示了很多强相关性.第 1 章的图 1.1 是散点图矩阵,显示了变量两两间的图象.我们看到 svi 是二进制变量,gleason 是有序分类变量.举个例子,我们看到 lcavol 和 lcp 都与响应变量 lpsa 有强的相关性.我们需要联合拟合这些影响来解开预测变量和响应变量之间的关系.
表 3.1 前列腺癌数据中预测变量的相关系数;表 3.2 前列腺癌数据的线性模型拟合. Z Z Z 分数为系数除以标准误差 \eqref{3.12}.大体上, Z Z Z 分数绝对值大于 2 表示在 p = 0.05 p=0.05 p=0.05 的水平下显著不为零.
我们对前列腺特定抗原的对数值 lpsa 和经过一次标准化有单位方差的预测变量进行线性模型拟合.我们随机将数据集分离得到规模为 67 的训练集以及规模为 30 的测试集.我们对训练集应用最小二乘估计得到表 3.2 中的估计值、标准误差以及
Z
Z
Z 分数.
Z
Z
Z 分数定义式为 (3.12),并且衡量了从模型中去掉该变量的影响.在
5
%
5\%
5% 的水平下,
Z
Z
Z 分数的绝对值大于 2 近似显著.(在我们的例子中,我们有 9 个参数,
t
67
−
9
t_{67-9}
t67−9 的分布的 0.025 尾分位数是
±
2.002
\pm 2.002
±2.002!)预测变量 lcavol
显示了最强的影响,lweight 和 svi 的影响同样显著.注意到一旦 lcavol
在模型中,lcp
便不显著(当使用一个没有 lcavol
的模型,lcp
是强显著的).我们也可以用
F
F
F 统计量 \eqref{3.13} 一次除去一些项.举个例子,我们考虑除掉表 3.2 中所有不显著的项,age
,lcp
,gleason
以及 pgg45
.我们得到
F = ( 32.81 − 29.43 ) / ( 9 − 5 ) 29.43 / ( 67 − 9 ) = 1.67 (3.16) F=\dfrac{(32.81-29.43)/(9-5)}{29.43/(67-9)}=1.67\tag{3.16} F=29.43/(67−9)(32.81−29.43)/(9−5)=1.67(3.16)
其中, p p p 值为 0.17( Pr ( F 4 , 58 > 1.67 ) = 0.17 \Pr(F_{4,58}>1.67)=0.17 Pr(F4,58>1.67)=0.17),因此不是显著的.
测试数据的平均预测误差为 0.521.相反地,使用 lpsa
训练数据的平均值预测的测试误差为 1.057,称作“基本误差阶”.因此线性模型将基本误差阶降低了大概
50
%
50\%
50%.下文中我们将要回到这个例子来比较不同的选择和收缩方法.
Guass-Markov 定理
统计学中一个很有名的结论称参数 β \beta β 的最小二乘估计在所有的线性无偏估计中有最小的方差.我们将要严谨地说明这一点,并且说明约束为无偏估计不是明智的选择.这个结论导致我们考虑本章中像岭回归的有偏估计.我们考虑任意参数 θ = a T β \theta = a^T\beta θ=aTβ 的线性组合,举个例子,预测值 f ( x 0 ) = x 0 T β f(x_0)=x_0^T\beta f(x0)=x0Tβ 是这种形式. a T β a^T\beta aTβ 的最小二乘估计为
θ ^ = a T β ^ = a T ( X T X ) − 1 X T y (3.17) \hat{\theta}=a^T\hat{\beta}=a^T\mathbf{(X^TX)^{-1}X^Ty}\tag{3.17} θ^=aTβ^=aT(XTX)−1XTy(3.17)
考虑固定
X
\mathbf{X}
X,则上式是关于响应变量
y
\mathbf{y}
y 的线性函数
c
0
T
y
\mathbf{c}_0^T\mathbf{y}
c0Ty.如果我们假设线性模型是正确的,则
a
T
β
^
a^T\hat{\beta}
aTβ^ 是无偏的,因为
E
(
a
T
β
^
)
=
E
(
a
T
(
X
T
X
)
−
1
X
T
y
)
=
a
T
(
X
T
X
)
−
1
X
T
X
β
=
a
T
β
(3.18)
\begin{aligned} \rm{E}(a^T\hat{\beta})&=E(a^T\mathbf{(X^TX)^{-1}X^Ty})\\ &=a^T\mathbf{(X^TX)^{-1}X^TX}\beta\\ &=a^T\beta\tag{3.18} \end{aligned}
E(aTβ^)=E(aT(XTX)−1XTy)=aT(XTX)−1XTXβ=aTβ(3.18)
Gauss-Markov 定理说如果
a
T
β
a^T\beta
aTβ 存在其它无偏估计
θ
~
=
c
T
y
\tilde{\theta}=\mathbf{c}^T\mathbf{y}
θ~=cTy,即
E
(
c
T
y
)
=
a
T
β
\rm{E}(\mathbf{c}^T\mathbf{y})=a^T\beta
E(cTy)=aTβ,则
V
a
r
(
a
T
β
^
)
≤
V
a
r
(
c
T
y
)
(3.19)
Var(a^T\hat{\beta})\le Var(\mathbf{c}^T\mathbf{y})\tag{3.19}
Var(aTβ^)≤Var(cTy)(3.19)
该证明(习题3.3a)运用三角不等式.为了简化我们已经用单个参数 a T β a^T\beta aTβ 估计来叙述结果,但若额外增加习题3.3b中的定义,则可以用整个参数向量 β \beta β 来说明Gauss-Markov 定理.
!!! note “weiya注:Ex. 3.3”
已完成,详见 Issue 70: Ex. 3.3
考虑 θ \theta θ 的估计值 θ ~ \tilde{\theta} θ~ 的均方误差
M S E ( θ ~ ) = E ( θ ~ − θ ) 2 = V a r ( θ ~ ) + [ E ( θ ~ ) − θ ] 2 (3.20) \begin{aligned} \rm{MSE}(\tilde{\theta})&=\rm{E}(\tilde{\theta}-\theta)^2\\ &=\rm{Var}(\tilde{\theta})+[\rm{E}(\tilde{\theta})-\theta]^2\tag{3.20} \end{aligned} MSE(θ~)=E(θ~−θ)2=Var(θ~)+[E(θ~)−θ]2(3.20)
第一项为方差,第二项为平方偏差.Gauss-Markov 定理表明最小二乘估计在所有无偏线性估计中有最小的均方误差.然而,或许存在有较小均方误差的有偏估计.这样的估计用小的偏差来换取方差大幅度的降低.实际中也会经常使用有偏估计.任何收缩或者将最小二乘的一些参数设为 0 的方法都可能导致有偏估计.我们将在这章的后半部分讨论许多例子,包括 变量子集选择 和 岭回归.从一个更加实际的观点来看,许多模型是对事实的曲解,因此是有偏的;挑选一个合适的模型意味着要在偏差和方差之间创造某种良好的平衡.我们将在第 7 章中详细讨论这些问题.
正如在第 2 章中讨论的那样,均方误差与预测的正确性紧密相关.考虑在输入 x 0 x_0 x0 处的新的响应变量
Y 0 = f ( x 0 ) + ε 0 (3.21) Y_0=f(x_0)+\varepsilon_0 \tag{3.21} Y0=f(x0)+ε0(3.21)
其估计量 f ~ ( x 0 ) = x 0 T β ~ \tilde{f}(x_0)=x_0^T\tilde{\beta} f~(x0)=x0Tβ~ 的期望预测误差为
E ( Y 0 − f ~ ( x 0 ) ) 2 = σ 2 + E ( x 0 T β ~ − f ( x 0 ) ) 2 = σ 2 + M S E ( f ~ ( x 0 ) ) (3.22) \begin{aligned} \rm{E}(Y_0-\tilde{f}(x_0))^2 &=\sigma^2+E(x_0^T\tilde{\beta}-f(x_0))^2\\ &= \sigma^2+MSE(\tilde{f}(x_0))\tag{3.22} \end{aligned} E(Y0−f~(x0))2=σ2+E(x0Tβ~−f(x0))2=σ2+MSE(f~(x0))(3.22)
因此,预测误差的估计值和均方误差只有常数值 σ 2 \sigma^2 σ2 的差别, σ 2 \sigma^2 σ2 表示了新观测 y 0 y_0 y0 的方差.