柚子快報(bào)邀請(qǐng)碼778899分享：泊松回歸統(tǒng)計(jì)模型

http://yzkb.51969.com/

背景知識(shí)

廣義線性模型（GLM）是線性模型的擴(kuò)展，通過(guò)聯(lián)結(jié)函數(shù)建立響應(yīng)變量的數(shù)學(xué)期望值與線性組合的預(yù)測(cè)變量之間的關(guān)系。其特點(diǎn)是不強(qiáng)行改變數(shù)據(jù)的自然度量，數(shù)據(jù)可以具有非線性和非恒定方差結(jié)構(gòu)。是線性模型在研究響應(yīng)值的非正態(tài)分布以及非線性模型簡(jiǎn)潔直接的線性轉(zhuǎn)化時(shí)的一種發(fā)展。

廣義線性模型定義 指數(shù)分布族定義 指數(shù)分布族實(shí)例

泊松回歸

泊松回歸常被用于計(jì)數(shù)數(shù)據(jù)情境下的建模。計(jì)數(shù)數(shù)據(jù)不能連續(xù)取值，或只能得到0，1，2，3等自然數(shù)的數(shù)據(jù)。任何兩個(gè)計(jì)數(shù)數(shù)據(jù)之間不能再細(xì)分，只能是整數(shù)?？煞譃橛浖?shù)據(jù)和記點(diǎn)數(shù)據(jù)，前者指按件計(jì)數(shù)的數(shù)據(jù)，如不合格品數(shù)、冰箱數(shù)、質(zhì)量檢測(cè)項(xiàng)目數(shù)等；后者指按項(xiàng)計(jì)數(shù)的數(shù)據(jù)，如疾病發(fā)病數(shù)、疵點(diǎn)數(shù)、氣泡數(shù)、產(chǎn)品缺陷數(shù)等。記件數(shù)據(jù)一般服從二項(xiàng)分布，記點(diǎn)數(shù)據(jù)一般服從泊松分布。

泊松回歸的模型形式為：

其中表示形式為計(jì)數(shù)數(shù)據(jù)的響應(yīng)變量，表示由一組相互獨(dú)立的解釋變量變量組成的向量，表示待估參數(shù)，可根據(jù)極大似然估計(jì)法進(jìn)行估計(jì)。

根據(jù)泊松分布的概率密度函數(shù)：

我們可以求出關(guān)于的極大似然函數(shù)（其中m表示觀測(cè)樣本數(shù)）：

進(jìn)而得到對(duì)數(shù)似然函數(shù)：

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

極大似然估計(jì)不能通過(guò)解析表達(dá)式獲得解析解，由其對(duì)數(shù)似然函數(shù)為凸函數(shù)的特性，可通過(guò)Newton–Raphson或其他基于梯度下降的方法進(jìn)行計(jì)算。

R語(yǔ)言實(shí)現(xiàn)泊松回歸

本次建模我們使用的是由Vnabuls和Ripley（2003）提供的從1981年起比利時(shí)每年新艾滋病病例數(shù)量的數(shù)據(jù)。具體的模型表達(dá)式如下：

我們假設(shè)～，其中是年的新病例數(shù)。假i是獨(dú)立的。和是待估參數(shù)。

#?數(shù)據(jù)?&?可視化

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))

建模時(shí)我們使用R語(yǔ)言stats包中的glm函數(shù)。

#?建立泊松回歸模型，可以看出代碼非常簡(jiǎn)潔

m0?<-?glm(y1~t,?family=poisson)?#?~左邊表示響應(yīng)變量，右邊表示解釋變量，family選取泊松分布

summary(m0)

plot(m0)

使用summary函數(shù)查看建模結(jié)果，可知和都很顯著。但是殘差和AIC都比較高。

從殘差圖可以看出違反了獨(dú)立性假設(shè)，可能是由于模型中遺漏了一些重要的變量，我們嘗試在模型中加入一個(gè)二次項(xiàng)再進(jìn)行建模

#?加入二次項(xiàng)建模

m1?<-?glm(y~t+t^2,family=poisson)

plot(m1)

summary(m1)

anova(m0,m1,test="Chisq")?#?用于檢驗(yàn)兩個(gè)模型是否有顯著差異

從結(jié)果可以看出，殘差和AIC都得到了顯著下降，獨(dú)立性假設(shè)也得到了滿足。但是也要注意進(jìn)行變量選擇，防止出現(xiàn)過(guò)擬合。

最后我們計(jì)算了一下擬合的置信區(qū)間，并進(jìn)行了可視化展示，從下圖可以看出擬合效果還是不錯(cuò)的。

#?計(jì)算置信區(qū)間并畫圖

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)

新冠預(yù)測(cè)練習(xí)

文章最后，小編為大家準(zhǔn)備了全球的新冠肺炎發(fā)病數(shù)據(jù)，大家可以嘗試一下使用泊松模型自己建模哦~

柚子快報(bào)邀請(qǐng)碼778899分享：泊松回歸統(tǒng)計(jì)模型

http://yzkb.51969.com/

文章鏈接