一般我们了解的线性模型是针对连续性变量,并且服从正态分布的,但是在实际应用上显得非常的局限。因为我们我看到的数据很多都是离散的,而且不是服从正态分布的。针对这种情况,对传统线性模型进行推广,行成了现在的广义线性模型。广义线性模型使得变量从正态分布拓展到指数分布族,从连续型变量拓展到离散型变量,这就使得在现实中有着很好的运用,下面开始介绍广义线性模型。
广义线性模型(GLM)定义
由以下三部分组成:
1 随机部分
随机样本
Y
1
,
Y
2
,
.
.
.
,
Y
n
Y_{1},Y_{2},...,Y_{n}
Y1,Y2,...,Yn服从的分布来自指数分布族,即
Y
i
Y_{i}
Yi的分布形式为
f
(
y
i
;
θ
i
,
ϕ
)
=
e
x
p
{
y
i
θ
i
−
b
(
θ
i
)
a
(
ϕ
)
+
c
(
y
i
,
ϕ
)
}
f\left ( y_{i};\theta _{i},\phi \right )=exp\left \{ \frac{y_{i}\theta _{i}-b\left ( \theta _{i} \right )}{a\left ( \phi \right )}+c\left ( y_{i},\phi \right ) \right \}
f(yi;θi,ϕ)=exp{a(ϕ)yiθi−b(θi)+c(yi,ϕ)}
其中参数
θ
i
\theta_{i}
θi叫做正则参数,并且随着指数
i
(
i
=
1
,
2
,
.
.
.
,
n
)
i(i = 1,2,...,n)
i(i=1,2,...,n)而变化,但是扰乱因子
ϕ
\phi
ϕ 是个常数。
2 系统部分
对于第
i
i
i个观测
Y
i
Y_{i}
Yi,我们有一个称为系统部分的线性预测值,即所研究变量的线性组合,即
η
i
=
x
i
T
β
=
∑
j
=
1
p
x
i
j
β
j
,
i
=
1
,
2
,
.
.
.
n
\eta _{i}=x_{i}^{T}\beta =\sum_{j=1}^{p}x_{ij}\beta _{j},i=1,2,...n
ηi=xiTβ=j=1∑pxijβj,i=1,2,...n
3 连接函数
有一个单调可微函数
g
(
)
g\left ( \right )
g()称为连接函数,它将随机部分的期望和系统部分连接起来,通过下面的等式
g
(
μ
i
)
=
η
i
=
x
i
T
β
,
i
=
1
,
2
,
.
.
.
n
,
g\left ( \mu_{i} \right )=\eta _{i}=x_{i}^{T}\beta ,i=1,2,...n,
g(μi)=ηi=xiTβ,i=1,2,...n,其中
μ
i
=
E
(
Y
i
)
\mu_{i}=E\left ( Y_{i} \right )
μi=E(Yi)是
Y
i
Y_{i}
Yi的期望。
矩阵表示:
η
=
[
η
1
η
2
⋮
η
n
]
n
×
1
,
μ
=
[
μ
1
μ
2
⋮
μ
n
]
n
×
1
,
X
=
[
x
1
′
x
2
′
⋮
x
n
′
]
n
×
p
\eta =\begin{bmatrix}\eta _{1}\\ \eta _{2}\\ \vdots \\ \eta _{n}\end{bmatrix}_{n\times 1},\mu=\begin{bmatrix}\mu _{1}\\ \mu _{2}\\ \vdots \\ \mu _{n}\end{bmatrix}_{n\times 1},X=\begin{bmatrix}x_{1}^{'}\\ x_{2}^{'}\\ \vdots \\ x_{n}^{'}\end{bmatrix}_{n\times p}
η=⎣⎢⎢⎢⎡η1η2⋮ηn⎦⎥⎥⎥⎤n×1,μ=⎣⎢⎢⎢⎡μ1μ2⋮μn⎦⎥⎥⎥⎤n×1,X=⎣⎢⎢⎢⎡x1′x2′⋮xn′⎦⎥⎥⎥⎤n×p
那么连接函数可以用矩阵形式表示
g
(
μ
)
=
η
=
X
β
g\left ( \mu \right )=\eta =X\beta
g(μ)=η=Xβ
连接函数介绍
1、正如
G
L
M
s
GLMs
GLMs的定义所指出的那样,连接函数是单调可微的,用于连接随机部分的期望和系统部分的线性预测值
2、选择与分布相关的自然参数作为连接函数,在这种情况下,它被称为点则连接。具体而言,如果连接函数
g
(
)
g()
g()采用与自然参数相同的函数形式,那么它被称为点则连接函数。
3、点则连接的优点是它可以带来非常好的统计特性,并且使用起来很方便。例如,对于最常用的分布,我们有以下点则连接函数。
Normal | μ = η \mu =\eta μ=η | (identity-link) |
---|---|---|
Poisson | l o g μ = η log\mu =\eta logμ=η | (log-link) |
Bernoulli | l o g π 1 − π = η log\frac{\pi}{1-\pi}=\eta log1−ππ=η | (logistic-link) |
Binomial | l o g π 1 − π = η log\frac{\pi}{1-\pi}=\eta log1−ππ=η | (logistic-link) |
4、然而,点则连接并不是连接函数的唯一选择。其他可能的
G
L
M
s
GLMs
GLMs连接函数包括
(a)二项分布的Probit连接:
η
=
Φ
−
1
(
π
)
\eta =\Phi ^{-1}\left ( \pi \right )
η=Φ−1(π);
0
<
π
<
1
0<\pi<1
0<π<1,其中
Φ
(
)
\Phi()
Φ()叫做累计分布函数(不是点则连接呦)
(b)补充的二项分布的log-log连接
η
=
l
o
g
{
−
l
o
g
(
1
−
π
)
}
,
0
<
π
<
1
\eta =log\left \{ -log\left ( 1-\pi \right ) \right \},0<\pi<1
η=log{−log(1−π)},0<π<1
(c)属于任何幂族分布的连接
η
=
{
μ
λ
,
i
f
λ
≠
0
l
o
g
μ
,
i
f
λ
=
0
\eta =\left\{\begin{matrix}\mu ^{\lambda }, if \lambda \neq 0 & \\ log\mu ,if \lambda =0&\end{matrix}\right.
η={μλ,ifλ=0logμ,ifλ=0
最大似然估计(MLE)的一般原则
假设我们对未知参数
θ
\theta
θ的对数似然函数,比如说
l
(
θ
)
l\left ( \theta \right )
l(θ)我们想找出
θ
\theta
θ的最大似然估计(MLE)
θ
^
\hat{\theta }
θ^,即
θ
^
≡
a
r
g
m
a
x
θ
⊂
Ω
{
l
(
θ
)
}
\hat{\theta }\equiv arg \underset{\theta \subset \Omega }{max}\left \{ l\left ( \theta \right ) \right \}
θ^≡argθ⊂Ωmax{l(θ)}
这是估计方程的解
∂
l
(
θ
)
∂
θ
=
0
\frac{\partial l\left ( \theta \right )}{\partial \theta }=0
∂θ∂l(θ)=0
1、在这种情况下,例如,对于正态分布参数
θ
\theta
θ的最大似然估计
θ
^
\hat{\theta }
θ^可以有一个明确的数学表达式
(
μ
=
1
n
∑
i
=
1
n
l
n
x
i
)
(\mu =\frac{1}{n}\sum_{i=1}^{n}lnx_{i})
(μ=n1∑i=1nlnxi)
2、一般来说,最大似然估计
θ
\theta
θ没有没有明确的数学解。相反,需要某些数值优化方法。
3、统计学中最常用的两种优化方法是Newton-Raphson算法和Fisher得分算法,他们都涉及计算
l
(
θ
)
l\left ( \theta \right )
l(θ)对
θ
\theta
θ的2次偏导数。
Newton-Raphson算法
该算法计算最大似然估计
θ
^
\hat{\theta }
θ^,通过下面的迭代:
θ
m
=
θ
m
−
1
+
[
−
l
′
′
(
θ
(
m
−
1
)
)
]
−
1
[
l
′
(
θ
(
m
−
1
)
)
]
(
1
)
\theta ^{m}=\theta ^{m-1}+\left [ -l^{''} \left ( \theta ^{\left ( m-1 \right )} \right )\right ]^{-1}\left [ l^{'} \left ( \theta ^{\left ( m-1 \right )} \right )\right ](1)
θm=θm−1+[−l′′(θ(m−1))]−1[l′(θ(m−1))](1)
其中
m
=
1
,
2
,
.
.
.
m=1,2,...
m=1,2,...这里,
l
′
(
θ
(
m
−
1
)
)
=
∂
l
(
θ
)
∂
θ
∣
θ
=
θ
(
m
−
1
)
l^{'}\left ( \theta ^{\left ( m-1 \right )} \right )=\frac{\partial l\left ( \theta \right )}{\partial \theta }|_{\theta =\theta ^{\left ( m-1 \right )}}
l′(θ(m−1))=∂θ∂l(θ)∣θ=θ(m−1)
l
′
′
(
θ
(
m
−
1
)
)
=
∂
2
l
(
θ
)
∂
θ
∂
θ
T
∣
θ
=
θ
(
m
−
1
)
l^{''}\left ( \theta ^{\left ( m-1 \right )} \right )=\frac{\partial ^{2}l\left ( \theta \right )}{\partial \theta \partial \theta ^{T}}|_{\theta =\theta ^{\left ( m-1 \right )}}
l′′(θ(m−1))=∂θ∂θT∂2l(θ)∣θ=θ(m−1)是
p
×
1
p\times 1
p×1和
p
×
p
p\times p
p×p的向量和矩阵
(
p
是
θ
的
维
数
)
(p是\theta的维数)
(p是θ的维数)
注1
l
′
(
θ
)
被
称
为
θ
的
得
分
函
数
。
−
l
′
′
(
θ
)
被
称
为
θ
的
观
测
信
息
矩
阵
l^{'}\left ( \theta \right )被称为\theta的得分函数。-l^{''}\left ( \theta \right )被称为\theta的观测信息矩阵
l′(θ)被称为θ的得分函数。−l′′(θ)被称为θ的观测信息矩阵
注2 算法(1)需要初始值,例如
θ
0
\theta ^{0}
θ0,以开始迭代过程。初始值的选择通常需要经验。
注3 算法(1)迭代直到收敛。例如,当迭代结果满足
∥
θ
(
m
)
−
θ
(
m
−
1
)
∥
∥
θ
(
m
−
1
)
∥
≤
1
0
−
5
\frac{\left \| \theta ^{\left ( m \right )}-\theta ^{\left ( m-1 \right )} \right \|}{\left \| \theta ^{\left ( m-1 \right )}\right \|}\leq 10^{-5}
∥∥θ(m−1)∥∥∥∥θ(m)−θ(m−1)∥∥≤10−5则迭代停止。
θ
(
m
)
\theta ^{\left ( m \right )}
θ(m)可以认为是最大似然估计
θ
^
\hat{\theta }
θ^。
Fisher得分算法
Fisher得分算法与Newton-Raphson算法相同,只是(1)式中的观测矩阵
−
l
′
′
(
θ
)
-l^{''}\left ( \theta \right )
−l′′(θ)被Fisher信息矩阵所代替
I
(
θ
)
=
E
[
−
l
′
′
(
θ
)
]
=
−
∫
l
′
′
(
θ
∣
Y
)
f
Y
(
Y
∣
θ
)
d
Y
I\left ( \theta \right )=E\left [ -l^{''}\left ( \theta \right ) \right ]=-\int l^{''}\left ( \theta |Y \right )f_{Y}\left ( Y|\theta \right )dY
I(θ)=E[−l′′(θ)]=−∫l′′(θ∣Y)fY(Y∣θ)dY
注释 Fisher得分算法和Newton-Raphson算法一般收敛于同一解。对于前者,在某些情况下,在信息矩阵的解析式上可能比后者更简单。例如Fisher信息矩阵可能是对角阵或者块对角阵,二观测信息矩阵可能不是。其次作为副产物,这两种算法都提供了极大似然估计的协方差矩阵。
广义线性模型(GLMs)中的最大似然估计(MLE)
首先,GLMs中的对数似然函数具有这样的形式
l
=
∑
i
=
1
n
l
o
g
f
(
y
i
;
θ
i
,
ϕ
)
=
∑
i
=
1
n
(
y
i
θ
i
−
b
(
θ
i
)
)
a
i
(
ϕ
)
+
∑
i
=
1
n
c
(
y
i
,
ϕ
)
l=\sum_{i=1}^{n}logf\left ( y_{i};\theta _{i} ,\phi \right )=\sum_{i=1}^{n}\frac{\left ( y_{i}\theta _{i}-b\left ( \theta _{i} \right ) \right )}{a_{i}\left ( \phi \right )}+\sum_{i=1}^{n}c\left ( y_{i} ,\phi \right )
l=i=1∑nlogf(yi;θi,ϕ)=i=1∑nai(ϕ)(yiθi−b(θi))+i=1∑nc(yi,ϕ)其中,
β
=
(
β
1
,
β
2
,
.
.
.
,
β
p
)
T
\beta =\left ( \beta _{1} ,\beta _{2},..., \beta _{p}\right )^{T}
β=(β1,β2,...,βp)T,
β
j
\beta_{j}
βj的得分函数为
U
j
=
∂
l
∂
β
j
=
∑
i
=
1
n
(
y
i
−
b
′
(
θ
i
)
)
a
i
(
ϕ
)
∂
θ
i
∂
β
j
=
∑
i
=
1
n
(
y
i
−
μ
)
a
i
(
ϕ
)
∂
θ
i
∂
β
j
(
2
)
U_{j}=\frac{\partial l}{\partial \beta _{j}}=\sum_{i=1}^{n}\frac{\left ( y_{i}-b^{'}\left ( \theta _{i} \right ) \right )}{a_{i}\left ( \phi \right )}\frac{\partial \theta _{i}}{\partial \beta _{j}}=\sum_{i=1}^{n}\frac{\left ( y_{i}-\mu \right )}{a_{i}\left ( \phi \right )}\frac{\partial \theta _{i}}{\partial \beta _{j}}(2)
Uj=∂βj∂l=i=1∑nai(ϕ)(yi−b′(θi))∂βj∂θi=i=1∑nai(ϕ)(yi−μ)∂βj∂θi(2)
其中,
μ
i
=
E
(
Y
i
)
=
b
′
(
θ
i
)
,
V
a
r
(
Y
i
)
=
b
′
′
(
θ
i
)
a
(
ϕ
)
\mu _{i}=E\left ( Y_{i} \right )=b^{'}\left ( \theta _{i} \right ),Var\left ( Y_{i} \right )=b^{''}\left ( \theta _{i} \right )a\left ( \phi \right )
μi=E(Yi)=b′(θi),Var(Yi)=b′′(θi)a(ϕ),我们使用链式法则进行差异化
∂
θ
i
∂
β
j
=
∂
θ
i
∂
μ
i
∂
μ
i
∂
β
j
\frac{\partial \theta _{i}}{\partial \beta _{j}}=\frac{\partial \theta _{i}}{\partial \mu _{i}}\frac{\partial \mu _{i}}{\partial \beta _{j}}
∂βj∂θi=∂μi∂θi∂βj∂μi因为
∂
θ
i
∂
μ
i
=
1
∂
μ
i
∂
θ
i
=
1
b
′
′
(
θ
i
)
=
a
i
(
ϕ
)
b
′
′
(
θ
i
)
a
i
(
ϕ
)
=
a
i
(
ϕ
)
V
a
r
(
Y
i
)
\frac{\partial \theta _{i}}{\partial \mu _{i}}=\frac{1}{\frac{\partial \mu _{i}}{\partial \theta _{i}}}=\frac{1}{b^{''}\left ( \theta _{i} \right )}=\frac{a_{i}\left ( \phi \right )}{b^{''}\left ( \theta _{i} \right )a_{i}\left ( \phi \right )}=\frac{a_{i}\left ( \phi \right )}{Var\left ( Y_{i} \right )}
∂μi∂θi=∂θi∂μi1=b′′(θi)1=b′′(θi)ai(ϕ)ai(ϕ)=Var(Yi)ai(ϕ)
并且
∂
μ
i
∂
β
j
=
∂
μ
i
∂
η
i
∂
η
i
∂
β
j
=
∂
μ
i
∂
η
i
x
i
j
\frac{\partial \mu _{i}}{\partial \beta _{j}}=\frac{\partial \mu _{i}}{\partial \eta _{i}}\frac{\partial \eta _{i}}{\partial \beta _{j}}=\frac{\partial \mu _{i}}{\partial \eta _{i}}x_{ij}
∂βj∂μi=∂ηi∂μi∂βj∂ηi=∂ηi∂μixij其中
x
i
j
x_{ij}
xij是
x
i
的
第
j
个
分
量
x_{i}的第j个分量
xi的第j个分量,我们知道
∂
θ
i
∂
β
j
=
a
i
(
ϕ
)
V
a
r
(
Y
i
)
∂
μ
i
∂
η
j
x
i
j
\frac{\partial \theta _{i}}{\partial \beta _{j}}=\frac{a_{i}\left ( \phi \right )}{Var\left ( Y_{i} \right )}\frac{\partial \mu _{i}}{\partial \eta _{j}}x_{ij}
∂βj∂θi=Var(Yi)ai(ϕ)∂ηj∂μixij因此(2)式就化为了
U
j
=
∑
i
=
1
n
[
(
y
i
−
μ
i
)
V
a
r
(
Y
i
)
x
i
j
(
∂
μ
i
∂
η
i
)
]
=
∑
i
=
1
n
(
y
i
−
μ
i
)
g
′
(
μ
i
)
V
i
x
i
j
(
3
)
U_{j}=\sum_{i=1}^{n}\left [ \frac{\left ( y_{i}-\mu _{i} \right )}{Var\left ( Y_{i} \right )}x_{ij} \left ( \frac{\partial \mu _{i}}{\partial \eta _{i}} \right )\right ]=\sum_{i=1}^{n}\frac{\left ( y_{i} -\mu _{i}\right )}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{ij}(3)
Uj=i=1∑n[Var(Yi)(yi−μi)xij(∂ηi∂μi)]=i=1∑ng′(μi)Vi(yi−μi)xij(3)其中
V
i
=
V
a
r
(
Y
i
)
V_{i}=Var\left ( Y_{i} \right )
Vi=Var(Yi),并且
∂
μ
i
∂
η
i
=
1
∂
η
i
∂
μ
i
=
1
g
′
(
μ
i
)
\frac{\partial \mu _{i}}{\partial \eta _{i}}=\frac{1}{\frac{\partial \eta _{i}}{\partial \mu _{i}}}=\frac{1}{g^{'}\left ( \mu _{i}\right )}
∂ηi∂μi=∂μi∂ηi1=g′(μi)1由于
η
i
=
g
(
μ
i
)
\eta _{i}=g\left ( \mu _{i} \right )
ηi=g(μi),因此
β
\beta
β的得分向量是
U
≡
U
(
β
)
=
∑
i
=
1
n
(
y
i
−
μ
i
)
g
′
(
μ
i
)
V
i
x
i
(
4
)
U\equiv U\left ( \beta \right )=\sum_{i=1}^{n}\frac{\left ( y_{i}-\mu _{i} \right )}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{i}(4)
U≡U(β)=i=1∑ng′(μi)Vi(yi−μi)xi(4)另一方面,(3)式对
β
j
\beta_{j}
βj求偏导得
∂
2
l
∂
β
j
∂
β
k
=
∂
U
j
∂
β
k
=
∑
i
=
1
n
(
−
∂
μ
i
∂
β
k
)
1
g
′
(
μ
i
)
V
i
x
i
j
+
∑
i
=
1
n
(
y
i
−
μ
i
)
∂
[
1
g
′
(
μ
i
)
V
i
]
∂
β
k
x
i
j
(
5
)
\frac{\partial ^{2}l}{\partial \beta _{j}\partial \beta _{k}}=\frac{\partial U_{j}}{\partial \beta _{k}}=\sum_{i=1}^{n}\left ( -\frac{\partial \mu _{i}}{\partial \beta _{k}} \right )\frac{1}{g^{'}\left ( \mu _{i} \right )V_{i}}x_{ij}+\sum_{i=1}^{n}\left ( y_{i} -\mu _{i}\right )\frac{\partial \left [ \frac{1}{g^{'}\left ( \mu _{i} \right )V_{i}} \right ]}{\partial \beta _{k}}x_{ij}(5)
∂βj∂βk∂2l=∂βk∂Uj=i=1∑n(−∂βk∂μi)g′(μi)Vi1xij+i=1∑n(yi−μi)∂βk∂[g′(μi)Vi1]xij(5)由于
E
(
Y
i
−
μ
i
)
=
0
E\left ( Y_{i} -\mu _{i}\right )=0
E(Yi−μi)=0,所以(5)式的第二项在进行期望时就消失了。即Fisher信息阵的矩阵形式就变成了
I
(
β
)
=
E
(
∂
2
l
∂
β
∂
β
T
)
=
∑
i
=
1
n
1
g
′
(
μ
i
)
2
V
i
x
i
j
x
i
k
I\left ( \beta \right )=E\left ( \frac{\partial ^{2}l}{\partial \beta \partial \beta ^{T}} \right )=\sum_{i=1}^{n}\frac{1}{g^{'} \left ( \mu _{i} \right )^{2}V_{i}}x_{ij}x_{ik}
I(β)=E(∂β∂βT∂2l)=i=1∑ng′(μi)2Vi1xijxik因此当我们表示
W
i
=
1
g
′
(
μ
i
)
2
V
i
W_{i}=\frac{1}{g^{'}\left ( \mu _{i} \right )^{2}V_{i}}
Wi=g′(μi)2Vi1,并且
W
=
d
i
a
g
(
W
1
,
W
2
,
.
.
.
,
W
n
)
=
(
W
1
0
⋯
0
0
W
2
⋯
0
⋮
0
⋱
0
0
⋯
0
W
n
)
W=diag\left ( W_{1},W_{2},...,W_{n} \right )=\begin{pmatrix} W_{1} &0 & \cdots & 0\\ 0& W_{2} & \cdots & 0\\ \vdots & 0 & \ddots & 0\\ 0& \cdots &0 & W_{n} \end{pmatrix}
W=diag(W1,W2,...,Wn)=⎝⎜⎜⎜⎛W10⋮00W20⋯⋯⋯⋱0000Wn⎠⎟⎟⎟⎞则Fisher信息阵就可以表示为
I
(
β
)
=
X
T
W
X
I\left ( \beta \right )=X^{T}WX
I(β)=XTWX令
D
=
d
i
a
g
(
g
′
(
μ
1
)
,
g
′
(
μ
2
)
,
.
.
.
,
g
′
(
μ
n
)
)
D=diag\left ( g^{'}\left ( \mu _{1} \right ) ,g^{'}\left ( \mu _{2} \right ),...,g^{'}\left ( \mu _{n} \right )\right )
D=diag(g′(μ1),g′(μ2),...,g′(μn)),这样(4)式就可以写成
U
=
U
(
β
)
=
X
T
W
D
(
y
−
μ
)
U=U\left ( \beta \right )=X^{T}WD\left ( y-\mu \right )
U=U(β)=XTWD(y−μ)
计算最大似然估计(MLE)参数 β \beta β的算法
假设我们有一个估计
β
(
m
−
1
)
\beta ^{\left ( m-1 \right )}
β(m−1),基于这个估计我们计算
μ
(
m
−
1
)
=
μ
(
β
(
m
−
1
)
)
,
W
(
m
−
1
)
=
W
(
β
(
m
−
1
)
)
\mu ^{\left ( m-1 \right )}=\mu \left ( \beta ^{\left ( m-1 \right )} \right ),W^{\left ( m-1\right )}=W\left ( \beta ^{\left ( m-1 \right )} \right )
μ(m−1)=μ(β(m−1)),W(m−1)=W(β(m−1))
并且有
D
(
m
−
1
)
=
D
(
β
(
m
−
1
)
)
D^{\left ( m-1 \right )}=D\left ( \beta ^{\left ( m-1 \right )} \right )
D(m−1)=D(β(m−1)) 那么Fisher得分算法就会显示
β
\beta
β的下一次迭代
β
(
m
)
=
β
(
m
−
1
)
+
[
I
(
β
(
m
−
1
)
)
]
−
1
[
U
(
β
(
m
−
1
)
)
]
=
β
(
m
−
1
)
+
[
X
T
W
(
M
−
1
)
X
]
−
1
[
X
T
W
(
M
−
1
)
D
(
M
−
1
)
(
y
−
μ
(
m
−
1
)
)
]
\beta ^{\left ( m \right )}=\beta ^{\left ( m-1 \right )}+\left [ I\left ( \beta ^{\left ( m-1 \right )} \right ) \right ]^{-1}\left [ U\left ( \beta ^{\left ( m-1 \right )} \right ) \right ]=\beta ^{\left ( m-1 \right )}+\left [ X^{T} W^{\left ( M-1 \right )}X\right ]^{-1}\left [ X^{T} W^{\left ( M-1 \right )}D^{\left ( M-1 \right )}\left ( y-\mu ^{\left ( m-1 \right )} \right )\right ]
β(m)=β(m−1)+[I(β(m−1))]−1[U(β(m−1))]=β(m−1)+[XTW(M−1)X]−1[XTW(M−1)D(M−1)(y−μ(m−1))]可以写成
β
(
m
)
=
[
X
T
W
X
]
−
1
X
T
W
(
m
−
1
)
[
X
β
(
m
−
1
)
+
D
(
m
−
1
)
(
y
−
μ
(
m
−
1
)
)
]
\beta ^{\left ( m \right )}=\left [ X^{T}WX \right ]^{-1}X^{T}W^{\left ( m-1 \right )}\left [ X\beta ^{\left ( m-1 \right )}+D^{\left ( m-1 \right )} \left ( y-\mu ^{\left ( m-1 \right )} \right )\right ]
β(m)=[XTWX]−1XTW(m−1)[Xβ(m−1)+D(m−1)(y−μ(m−1))]令
Z
(
m
−
1
)
=
X
β
(
m
−
1
)
+
D
(
m
−
1
)
(
y
−
μ
(
m
−
1
)
)
Z^{\left ( m-1 \right )}=X\beta ^{\left ( m-1 \right )}+D^{\left ( m-1 \right )}\left ( y-\mu ^{\left ( m-1 \right )} \right )
Z(m−1)=Xβ(m−1)+D(m−1)(y−μ(m−1))然后它又可以被写成
β
(
m
)
=
(
X
(
T
)
W
(
m
−
1
)
X
)
−
1
X
T
W
m
−
1
Z
(
m
−
1
)
(
6
)
\beta ^{\left ( m \right )}=\left ( X^{\left ( T \right )}W^{\left ( m-1 \right )}X \right )^{-1}X^{T}W^{m-1}Z^{\left ( m-1 \right )}(6)
β(m)=(X(T)W(m−1)X)−1XTWm−1Z(m−1)(6)
注释 (6)式意味着,给定参数
β
\beta
β的解,我们需要计算“工作权重矩阵”
W
W
W和“工作响应向量”
Z
Z
Z,然后利用广义加权最小二乘法得到
β
\beta
β的更新解。
广义线性模型实例解析
下表中的ARPI事物数据在协变量X的不同处观察到Y,并且数据是服从Poisson分布的。我们利用GLM来解决这个问题。
Y i Y_{i} Yi | 2 | 3 | 6 | 7 | 8 | 9 | 10 | 12 | 15 |
---|---|---|---|---|---|---|---|---|---|
x i x_{i} xi | -1 | -1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
数据即探索Y和X之间的关系。设
Y
i
Y_{i}
Yi为变量
y
y
y的第
i
i
i个数,表示
E
(
Y
i
)
=
μ
i
E\left ( Y_{i} \right )=\mu _{i}
E(Yi)=μi。我们通过建立关系
g
(
μ
i
)
=
x
i
′
β
g\left ( \mu _{i} \right )=x_{i}^{'}\beta
g(μi)=xi′β 对于这个Poisson数据集,点则连接是对数连接函数。
l
o
g
μ
i
=
β
0
+
β
1
x
i
=
(
1
,
x
i
)
(
β
0
β
1
)
=
x
i
T
β
log\mu _{i}=\beta _{0}+\beta _{1}x_{i}=\left ( 1,x_{i} \right )\begin{pmatrix} \beta _{0}\\ \beta _{1}\end{pmatrix}=x_{i}^{T}\beta
logμi=β0+β1xi=(1,xi)(β0β1)=xiTβ接下来我们要来求
W
和
Z
W和Z
W和Z的表达式。
我们已知的条件有
g
′
(
μ
i
)
=
1
μ
i
g^{'}\left ( \mu _{i} \right )=\frac{1}{\mu _{i}}
g′(μi)=μi1,对于Poisson分布显然有
V
i
=
E
(
Y
i
)
=
μ
i
V_{i}=E\left ( Y_{i} \right )=\mu _{i}
Vi=E(Yi)=μi,所以
W
i
=
[
(
g
′
(
μ
i
)
2
)
V
i
]
−
1
=
e
x
p
{
x
i
T
β
}
W_{i}=\left [ \left ( g^{'}\left ( \mu _{i} \right ) ^{2}\right )V_{i} \right ]^{-1}=exp\left \{ x_{i}^{T} \beta \right \}
Wi=[(g′(μi)2)Vi]−1=exp{xiTβ}并且
Z
i
=
x
i
T
β
+
g
′
(
μ
i
)
(
y
i
−
μ
i
)
=
x
i
T
β
+
(
y
i
−
μ
i
)
μ
i
Z_{i}=x_{i}^{T}\beta +g^{'}\left ( \mu _{i} \right )\left ( y_{i} -\mu _{i}\right )=x_{i}^{T}\beta +\frac{\left ( y_{i}-\mu _{i} \right )}{\mu _{i}}
Zi=xiTβ+g′(μi)(yi−μi)=xiTβ+μi(yi−μi)
我们选择
β
的
初
始
值
β
0
=
2
,
β
1
=
1
\beta的初始值\beta_{0}=2,\beta_{1}=1
β的初始值β0=2,β1=1。结合Fisher迭代算法,代入数据。这个过程一直持续到收敛。结果如下表
m | 0 | 1 | 2 | 3 | 4 |
---|---|---|---|---|---|
β 0 m \beta _{0}^{m} β0m | 2 | 1.9150 | 1.8902 | 1.8892 | 1.8892 |
β 1 m \beta _{1}^{m} β1m | 1 | 0.7235 | 1.8902 | 1.8892 | 1.8892 |
因此 β 的 M L E 是 β 0 = 1.8892 , β 1 = 0.6697 \beta的MLE是\beta _{0}=1.8892,\beta _{1}=0.6697 β的MLE是β0=1.8892,β1=0.6697
R语言代码
y <- c(2,3,6,7,8,9,10,12,15);
x <- c(-1,-1,0,0,0,0,1,1,1)
X <- cbind(rep(1,9),x); beta_0 <- c(2,1)
for (i in 1:100){
beta <- beta_0
eta <- X %*% beta
mu <- exp(eta)
W <- diag(as.vector(mu))
Z <- X %*% beta + ((y-mu)*mu^(-1))
XWX <- t(X) %*% W %*% X
XWZ <- t(X) %*% W %*% Z
Cov <- solve(XWX)
beta_0 <- Cov %*% XWZ}
testdata<-data.frame(y,x)
summary(glm(y~x,family=poisson,data=testdata))