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