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

泊松回归

定义

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

• e的概念
  
  $lim_{ntoinfty}(1-frac{1}{n})^n=e \ lim_{ntoinfty}(1-frac{lambda}{n})^n=e^lambda \$

  n→∞lim​(1−n1​)n=en→∞lim​(1−nλ​)n=eλ

• 二项式分布
  
  $P(X=k)=binom{n}{k}p^k(1-p)^{n-k}$

  P(X=k)=(kn​)pk(1−p)n−k
  如果我们令
  $p=lambda/n$

  p=λ/n,当
  $ntoinfty$

  n→∞时
  $P$

  P的极限是:
  
  $begin{aligned} lim_{ntoinfty}P(X=k)&=lim_{ntoinfty}binom{n}{k}p^k(1-p)^{n-k} \ &=lim_{ntoinfty}frac{n!}{(n-k)!k!} left( frac{lambda}{n} right)^kleft( 1-frac{lambda}{n} right)^{n-k} \ &=lim_{ntoinfty}left[frac{n!}{n^k(n-k)!}right] left( frac{lambda^k}{k!}right) underbrace{ left( 1-frac{lambda}{n} right) ^n }_{toexp(-lambda)}underbrace{left( 1-frac{lambda}{n} right) ^{-k} }_{to1}\ &=lim_{ntoinfty} 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}$

  n→∞lim​P(X=k)​=n→∞lim​(kn​)pk(1−p)n−k=n→∞lim​(n−k)!k!n!​(nλ​)k(1−nλ​)n−k=n→∞lim​[nk(n−k)!n!​](k!λk​)→exp(−λ)(1−nλ​)n​​→1(1−nλ​)−k​​=n→∞lim​[(1−n1​)(1−n2​)⋯(1−nk−1​)](k!λk​)(1−n1​)n(1−nλ​)−k=(k!λk​)exp{−λ}​
  这也说明当
  $ntoinfty$

  n→∞时,二项式分布可以近似至泊松分布

• 泊松分布与泊松回归
  
  $P(X=k)=frac{lambda^k}{k!}e^{-lambda},k=0,1,cdots$

  P(X=k)=k!λk​e−λ,k=0,1,⋯
  其中有
  $E(X)=lambda,Var(X)=lambda$

  E(X)=λ,Var(X)=λ,为导出泊松回归的形式,我们记上式为
  $Y$

  Y,则有:
  
  $Y=frac{lambda^k}{k!}e^{-lambda},k=0,1,cdots$

  Y=k!λk​e−λ,k=0,1,⋯
  做对数变换:
  
  $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^{log(Y;lambda)}=f(Y;lambda)=exp{klog(lambda)-lambda-log(k!) }$

  elog(Y;λ)=f(Y;λ)=exp{klog(λ)−λ−log(k!)}
  其中
  $log(Y)$

  log(Y)称为链接函数,由此得到了泊松回归模型的形式.那么为什么说泊松回归模型是一种广义线性模型(GLM)呢？首先考虑线性回归模型:
  
  $y=Xbeta+e \ eta=Xbeta\ mu=E(y)$

  y=Xβ+eη=Xβμ=E(y)
  在广义线性模型中,我们不再要求
  $y$

  y服从
  $N(mu,sigma^2)$

  N(μ,σ2)而是服从于指数分布族,例如Normal;Poisson;Gamma;Binomial.现证明泊松回归属于指数分布族,首先考虑指数分布族的定义：
  
  $f(y;theta,phi)=expleft{ frac{(ytheta-b(theta))}{a(phi)} +c(y,phi)right}$

  f(y;θ,ϕ)=exp{a(ϕ)(yθ−b(θ))​+c(y,ϕ)}
  现在我们回顾一下泊松回归模型的形式:
  
  $f(Y;lambda)=exp{klog(lambda)-lambda-log(k!) }$

  f(Y;λ)=exp{klog(λ)−λ−log(k!)}
  记
  $theta=log(lambda),lambda=exp{theta}$

  θ=log(λ),λ=exp{θ},则有：
  
  $begin{aligned} f(Y=k;theta,phi)&=exp{ytheta-exp{theta}-log(k!) } \ &=expleft{ frac{ytheta-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)中抽取容量为
  $n$

  n的随机样本,此时有似然函数:
  
  $L(beta)=prod_{i=1}^nf(beta;y_i,x_i)=prod_{i=1}^nexp{ y_iX_i^{rm{T}}beta_i-exp(X_i^{rm{T}}beta_i)-log(y_i!)}$

  L(β)=i=1∏n​f(β;yi​,xi​)=i=1∏n​exp{yi​XiT​βi​−exp(XiT​βi​)−log(yi​!)}
  对数似然:
  
  $l(beta)=sum_{i=1}^n[y_ix_ibeta_i-exp(x_ibeta_i)-log(y_i!)]$

  l(β)=i=1∑n​[yi​xi​βi​−exp(xi​βi​)−log(yi​!)]
  对
  $beta_i$

  βi​求偏导:
  
  $frac{partial l(beta)}{partialbeta_i}=sum_{i=1}^n[y_ix_i-x_iexp(x_ibeta_i)]=0$

  ∂βi​∂l(β)​=i=1∑n​[yi​xi​−xi​exp(xi​βi​)]=0
  上式方程组一般情况下并没有解析解,但我们可以通过牛顿拉夫逊迭代法求解:
  
  $F(beta)= begin{bmatrix} sum_{i=1}^n [y_i-exp(x_ibeta_i)] \ sum_{i=1}^n [y_ix_{i1}-x_{i1}exp(x_ibeta_i)] \ vdots\ sum_{i=1}^n [y_ix_{iq}-x_{iq}exp(x_ibeta_i)] end{bmatrix}$

  F(β)=⎣⎢⎢⎢⎡​∑i=1n​[yi​−exp(xi​βi​)]∑i=1n​[yi​xi1​−xi1​exp(xi​βi​)]⋮∑i=1n​[yi​xiq​−xiq​exp(xi​βi​)]​⎦⎥⎥⎥⎤​
  其中
  $beta_0=1,x_0=1$

  β0​=1,x0​=1.则
  $F(beta_i)$

  F(βi​)关于
  $beta$

  β的雅克比矩阵:
  
  $J(beta)=frac{partial^2l(beta)}{partialbeta_kpartialbeta_j}=-sum_{i=1}^nx_ix_{ik}exp(X^Tbeta),k=0,1,cdots,q;j=0,1,cdots,q.$

  J(β)=∂βk​∂βj​∂2l(β)​=−i=1∑n​xi​xik​exp(XTβ),k=0,1,⋯,q;j=0,1,⋯,q.
  此时有Newton-Raphson算法：
  
  $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$

  ε,不断计算上述过程直至
  $|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);