广义线性模型之泊松回归

2021年6月27日 3点热度 0条评论 来源: Infinity343

最近在研究GWPR,参考了很多广义线性模型,特别是泊松回归的相关内容,知识琐碎且繁杂,做个笔记。

泊松回归

定义

泊松回归(Poisson regression)是用来为计数资料和列联表建模的一种回归分析.泊松回归假设反应变量Y是泊松分布,并假设它期望值的对数可被未知参数的线性组合建模.泊松回归模型有时(特别是当用作列联表模型时)又被称作对数-线性模型.需要注意的是,对数线性模型和泊松回归模型并不完全相同,通常对数线性回归的响应变量是连续的,而泊松回归则是离散的.再给出泊松回归模型的形式之前,我们先考虑几个概念:

  • e的概念
    lim ⁡ n → ∞ ( 1 − 1 n ) n = e lim ⁡ n → ∞ ( 1 − λ n ) n = e λ \lim_{n\to\infty}(1-\frac{1}{n})^n=e \\ \lim_{n\to\infty}(1-\frac{\lambda}{n})^n=e^\lambda \\ nlim(1n1)n=enlim(1nλ)n=eλ

  • 二项式分布
    P ( X = k ) = ( n k ) p k ( 1 − p ) n − k P(X=k)=\binom{n}{k}p^k(1-p)^{n-k} P(X=k)=(kn)pk(1p)nk
    如果我们令 p = λ / n p=\lambda/n p=λ/n,当 n → ∞ n\to\infty n P P P的极限是:
    lim ⁡ n → ∞ P ( X = k ) = lim ⁡ n → ∞ ( n k ) p k ( 1 − p ) n − k = lim ⁡ n → ∞ n ! ( n − k ) ! k ! ( λ n ) k ( 1 − λ n ) n − k = lim ⁡ n → ∞ [ n ! n k ( n − k ) ! ] ( λ k k ! ) ( 1 − λ n ) n ⏟ → exp ⁡ ( − λ ) ( 1 − λ n ) − k ⏟ → 1 = lim ⁡ n → ∞ [ ( 1 − 1 n ) ( 1 − 2 n ) ⋯ ( 1 − k − 1 n ) ] ( λ k k ! ) ( 1 − 1 n ) n ( 1 − λ n ) − k = ( λ k k ! ) exp ⁡ { − λ } \begin{aligned} \lim_{n\to\infty}P(X=k)&=\lim_{n\to\infty}\binom{n}{k}p^k(1-p)^{n-k} \\ &=\lim_{n\to\infty}\frac{n!}{(n-k)!k!} \left( \frac{\lambda}{n} \right)^k\left( 1-\frac{\lambda}{n} \right)^{n-k} \\ &=\lim_{n\to\infty}\left[\frac{n!}{n^k(n-k)!}\right] \left( \frac{\lambda^k}{k!}\right) \underbrace{ \left( 1-\frac{\lambda}{n} \right) ^n }_{\to\exp(-\lambda)}\underbrace{\left( 1-\frac{\lambda}{n} \right) ^{-k} }_{\to1}\\ &=\lim_{n\to\infty} \left[ \left( 1-\frac{1}{n} \right) \left( 1-\frac{2}{n} \right)\cdots \left( 1-\frac{k-1}{n} \right)\right]\left( \frac{\lambda^k}{k!} \right) \left( 1-\frac{1}{n} \right)^n \left( 1-\frac{\lambda}{n} \right)^{-k} \\ &=\left( \frac{\lambda^k}{k!} \right)\exp\{-\lambda\} \end{aligned} nlimP(X=k)=nlim(kn)pk(1p)nk=nlim(nk)!k!n!(nλ)k(1nλ)nk=nlim[nk(nk)!n!](k!λk)exp(λ) (1nλ)n1 (1nλ)k=nlim[(1n1)(1n2)(1nk1)](k!λk)(1n1)n(1nλ)k=(k!λk)exp{ λ}
    这也说明当 n → ∞ n\to\infty n时,二项式分布可以近似至泊松分布

  • 泊松分布与泊松回归
    P ( X = k ) = λ k k ! e − λ , k = 0 , 1 , ⋯ P(X=k)=\frac{\lambda^k}{k!}e^{-\lambda},k=0,1,\cdots P(X=k)=k!λkeλ,k=0,1,
    其中有 E ( X ) = λ , V a r ( X ) = λ E(X)=\lambda,Var(X)=\lambda E(X)=λ,Var(X)=λ,为导出泊松回归的形式,我们记上式为 Y Y Y,则有:
    Y = λ k k ! e − λ , k = 0 , 1 , ⋯ Y=\frac{\lambda^k}{k!}e^{-\lambda},k=0,1,\cdots Y=k!λkeλ,k=0,1,
    做对数变换:
    l o g ( Y ; λ ) = l o g ( e − λ λ k ) − l o g ( k ! ) = l o g ( e − λ ) + l o g ( λ k ) − l o g ( k ! ) = − λ + k l o g ( λ ) − l o g ( k ! ) \begin{aligned} log(Y;\lambda)&=log(e^{-\lambda}\lambda^k)-log(k!)\\ &=log(e^{-\lambda})+log(\lambda^k)-log(k!)\\ &=-\lambda+klog(\lambda)-log(k!) \\ \end{aligned} log(Y;λ)=log(eλλk)log(k!)=log(eλ)+log(λk)log(k!)=λ+klog(λ)log(k!)
    再做指数变换:
    e l o g ( Y ; λ ) = f ( Y ; λ ) = exp ⁡ { k l o g ( λ ) − λ − l o g ( k ! ) } e^{log(Y;\lambda)}=f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} elog(Y;λ)=f(Y;λ)=exp{ klog(λ)λlog(k!)}
    其中 l o g ( Y ) log(Y) log(Y)称为链接函数,由此得到了泊松回归模型的形式.那么为什么说泊松回归模型是一种广义线性模型(GLM)呢?首先考虑线性回归模型:
    y = X β + e η = X β μ = E ( y ) y=X\beta+e \\ \eta=X\beta\\ \mu=E(y) y=Xβ+eη=Xβμ=E(y)
    在广义线性模型中,我们不再要求 y y y服从 N ( μ , σ 2 ) N(\mu,\sigma^2) N(μ,σ2)而是服从于指数分布族,例如Normal;Poisson;Gamma;Binomial.现证明泊松回归属于指数分布族,首先考虑指数分布族的定义:
    f ( y ; θ , ϕ ) = exp ⁡ { ( y θ − b ( θ ) ) a ( ϕ ) + c ( y , ϕ ) } f(y;\theta,\phi)=\exp\left\{ \frac{(y\theta-b(\theta))}{a(\phi)} +c(y,\phi)\right\} f(y;θ,ϕ)=exp{ a(ϕ)(yθb(θ))+c(y,ϕ)}
    现在我们回顾一下泊松回归模型的形式:
    f ( Y ; λ ) = exp ⁡ { k l o g ( λ ) − λ − l o g ( k ! ) } f(Y;\lambda)=\exp\{klog(\lambda)-\lambda-log(k!) \} f(Y;λ)=exp{ klog(λ)λlog(k!)}
    θ = l o g ( λ ) , λ = exp ⁡ { θ } \theta=log(\lambda),\lambda=\exp\{\theta\} θ=log(λ),λ=exp{ θ},则有:
    f ( Y = k ; θ , ϕ ) = exp ⁡ { y θ − exp ⁡ { θ } − l o g ( k ! ) } = exp ⁡ { y θ − exp ⁡ { θ } 1 + ( − ) l o g ( k ! ) } \begin{aligned} f(Y=k;\theta,\phi)&=\exp\{y\theta-\exp\{\theta\}-log(k!) \} \\ &=\exp\left\{ \frac{y\theta-\exp\{\theta\}}{1}+(-)log(k!)\right\} \end{aligned} f(Y=k;θ,ϕ)=exp{ yθexp{ θ}log(k!)}=exp{ 1yθexp{ θ}+()log(k!)}
    所以泊松回归模型也属于指数分布族

参数估计

  • 似然函数
    再给出了模型的形式以后,进一步需要对参数进行估计,采用极大似然估计法,从总体 ( Y , X ) (Y,X) (Y,X)中抽取容量为 n n n的随机样本,此时有似然函数:
    L ( β ) = ∏ i = 1 n f ( β ; y i , x i ) = ∏ i = 1 n exp ⁡ { y i X i T β i − exp ⁡ ( X i T β i ) − log ⁡ ( y i ! ) } L(\beta)=\prod_{i=1}^nf(\beta;y_i,x_i)=\prod_{i=1}^n\exp\{ y_iX_i^{\rm{T}}\beta_i-\exp(X_i^{\rm{T}}\beta_i)-\log(y_i!)\} L(β)=i=1nf(β;yi,xi)=i=1nexp{ yiXiTβiexp(XiTβi)log(yi!)}
    对数似然:
    l ( β ) = ∑ i = 1 n [ y i x i β i − exp ⁡ ( x i β i ) − log ⁡ ( y i ! ) ] l(\beta)=\sum_{i=1}^n[y_ix_i\beta_i-\exp(x_i\beta_i)-\log(y_i!)] l(β)=i=1n[yixiβiexp(xiβi)log(yi!)]
    β i \beta_i βi求偏导:
    ∂ l ( β ) ∂ β i = ∑ i = 1 n [ y i x i − x i exp ⁡ ( x i β i ) ] = 0 \frac{\partial l(\beta)}{\partial\beta_i}=\sum_{i=1}^n[y_ix_i-x_i\exp(x_i\beta_i)]=0 βil(β)=i=1n[yixixiexp(xiβi)]=0
    上式方程组一般情况下并没有解析解,但我们可以通过牛顿拉夫逊迭代法求解:
    F ( β ) = [ ∑ i = 1 n [ y i − exp ⁡ ( x i β i ) ] ∑ i = 1 n [ y i x i 1 − x i 1 exp ⁡ ( x i β i ) ] ⋮ ∑ i = 1 n [ y i x i q − x i q exp ⁡ ( x i β i ) ] ] F(\beta)= \begin{bmatrix} \sum_{i=1}^n [y_i-\exp(x_i\beta_i)] \\ \sum_{i=1}^n [y_ix_{i1}-x_{i1}\exp(x_i\beta_i)] \\ \vdots\\ \sum_{i=1}^n [y_ix_{iq}-x_{iq}\exp(x_i\beta_i)] \end{bmatrix} F(β)=i=1n[yiexp(xiβi)]i=1n[yixi1xi1exp(xiβi)]i=1n[yixiqxiqexp(xiβi)]
    其中 β 0 = 1 , x 0 = 1 \beta_0=1,x_0=1 β0=1,x0=1.则 F ( β i ) F(\beta_i) F(βi)关于 β \beta β的雅克比矩阵:
    J ( β ) = ∂ 2 l ( β ) ∂ β k ∂ β j = − ∑ i = 1 n x i x i k exp ⁡ ( X T β ) , k = 0 , 1 , ⋯   , q ; j = 0 , 1 , ⋯   , q . J(\beta)=\frac{\partial^2l(\beta)}{\partial\beta_k\partial\beta_j}=-\sum_{i=1}^nx_ix_{ik}\exp(X^T\beta),k=0,1,\cdots,q;j=0,1,\cdots,q. J(β)=βkβj2l(β)=i=1nxixikexp(XTβ),k=0,1,,q;j=0,1,,q.
    此时有Newton-Raphson算法:
    β ( m + 1 ) = β ( m ) − [ J ( β ( m ) ) ] − 1 F ( β ( m ) ) , m = 0 , 1 , ⋯ \beta^{(m+1)}=\beta^{(m)}-[J(\beta^{(m)})]^{-1}F(\beta^{(m)}),m=0,1,\cdots β(m+1)=β(m)[J(β(m))]1F(β(m)),m=0,1,
    在这个迭代过程中,需要给定 β \beta β的初值和精度 ε \varepsilon ε,不断计算上述过程直至 ∣ β ( m + 1 ) − β ( m ) ∣ < ε |\beta^{(m+1)}-\beta^{(m)}|<\varepsilon β(m+1)β(m)<ε收敛后结束.
    同时附上MATLAB代码
function F = PoissionRegressopt(b,Y,X)
n = length(Y);
F = 0;
for k = 1:n
    F = F + Y(k)*X(k,:)*b - exp(X(k,:)*b);% - factorial(Y(k));
end
F = - F;


function F = PoissionF(b,Y,X)
n = length(Y);
F = zeros(size(b));
for k = 1:n
    F = F + Y(k)*X(k,:)'- exp(X(k,:)*b)*X(k,:)';
end


function JM  = PoissionJM(b,Y,X)
n = length(Y);
JM  = zeros(size(b,1));
for k = 1:n
    JM = JM +  exp(X(k,:)*b)*X(k,:)'*X(k,:);
end

function [ bm fv1,fv2] = PoissionNR(bm0,Y,X)
itermax = 30;         
errstol = 1e-4;        
iters = 0;             
deltabm = ones(size(bm0));  
bm1 = bm0 + deltabm;    
while (iters<itermax)||(max(abs(deltabm))>errstol)
    deltabm  = pinv(PoissionJM(bm0,Y,X))*PoissionF(bm0,Y,X);
    bm1 = bm0 + deltabm;
    bm0  = bm1;   iters = iters +1;
end
bm = bm0;
fv1 = PoissionF(bm,Y,X);
fv2 = PoissionRegressopt(bm,Y,X);

    原文作者:Infinity343
    原文地址: https://blog.csdn.net/qq_44638724/article/details/106183848
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系管理员进行删除。