泊松回归统计模型_泊松回归模型

背景知识

广义线性模型（GLM）是线性模型的扩展，通过联结函数建立响应变量的数学期望值与线性组合的预测变量之间的关系。其特点是不强行改变数据的自然度量，数据可以具有非线性和非恒定方差结构。是线性模型在研究响应值的非正态分布以及非线性模型简洁直接的线性转化时的一种发展。

• 广义线性模型定义

  

• 指数分布族定义

  

• 指数分布族实例

  

泊松回归

泊松回归常被用于计数数据情境下的建模。计数数据不能连续取值，或只能得到0，1，2，3等自然数的数据。任何两个计数数据之间不能再细分，只能是整数。可分为记件数据和记点数据，前者指按件计数的数据，如不合格品数、冰箱数、质量检测项目数等；后者指按项计数的数据，如疾病发病数、疵点数、气泡数、产品缺陷数等。记件数据一般服从二项分布，记点数据一般服从泊松分布。

泊松回归的模型形式为：

其中表示形式为计数数据的响应变量，表示由一组相互独立的解释变量变量组成的向量，表示待估参数，可根据极大似然估计法进行估计。

根据泊松分布的概率密度函数：

我们可以求出关于的极大似然函数（其中m表示观测样本数）：

进而得到对数似然函数：

求解极大值需要解如下方程：

极大似然估计不能通过解析表达式获得解析解，由其对数似然函数为凸函数的特性，可通过Newton–Raphson或其他基于梯度下降的方法进行计算。

R语言实现泊松回归

本次建模我们使用的是由Vnabuls和Ripley（2003）提供的从1981年起比利时每年新艾滋病病例数量的数据。具体的模型表达式如下：

我们假设∼，其中是年的新病例数。假i是独立的。和是待估参数。

# 数据 & 可视化
y<- c(12,14,33,50,67,74,123,141,165,204,253,246,240)
t<-1:13
plot(t+1980,y,xlab="Year",ylab="New AIDS cases",ylim=c(0,280))

建模时我们使用R语言stats包中的glm函数。

# 建立泊松回归模型，可以看出代码非常简洁
m0 <- glm(y1~t, family=poisson) # ~左边表示响应变量，右边表示解释变量，family选取泊松分布
summary(m0)
plot(m0)

使用summary函数查看建模结果，可知和都很显著。但是残差和AIC都比较高。

从残差图可以看出违反了独立性假设，可能是由于模型中遗漏了一些重要的变量，我们尝试在模型中加入一个二次项再进行建模

# 加入二次项建模
m1 <- glm(y~t+t^2,family=poisson)
plot(m1)
summary(m1)
anova(m0,m1,test="Chisq") # 用于检验两个模型是否有显著差异

从结果可以看出，残差和AIC都得到了显著下降，独立性假设也得到了满足。但是也要注意进行变量选择，防止出现过拟合。

最后我们计算了一下拟合的置信区间，并进行了可视化展示，从下图可以看出拟合效果还是不错的。

# 计算置信区间并画图
beta.1 <- summary(m1)$coefficients[2,]
ci <- c(beta.1[1]-1.96*beta.1[2],beta.1[1]+1.96*beta.1[2])
ci # print 95% CI for beta_1
 
new.t<-seq(1,13,length=100)
fv <- predict(m1,data.frame(t=new.t),se=TRUE)
plot(t+1980,y,xlab="Year",ylab="New AIDS cases",ylim=c(0,280))
lines(new.t+1980,exp(fv$fit))
lines(new.t+1980,exp(fv$fit+2*fv$se.fit),lty=2)
lines(new.t+1980,exp(fv$fit-2*fv$se.fit),lty=2)

新冠预测练习

文章最后，小编为大家准备了全球的新冠肺炎发病数据，大家可以尝试一下使用泊松模型自己建模哦~