ポアソン回帰(比率データ)

ポアソン回帰はカウントデータあるいはイベントの発生率をモデル化する際に用いられる。このページでは、馬に蹴られて死亡したプロシア軍の死亡率について考える。

データは R の pscl パッケージに入っている。パッケージがない場合、インストールしておいてからロードする。

install.packages("pscl")
library(pscl)

実際のはデータは次のようにロードできる。14 軍団(corp)の 20 年間(year)の死亡者数(y)が記載されている。

data(prussian)
head(prussian)
##   y year corp
## 1 0   75    G
## 2 2   76    G
## 3 2   77    G
## 4 1   78    G
## 5 0   79    G
## 6 0   80    G

まず、20 年間で見た時、各軍団における死亡数(死亡率)に違いがあるかどうかについて解析する。

res <- glm(y ~ corp, family = poisson, offset = log(year), data = prussian)
summary(res)
## Call:
## glm(formula = y ~ corp, family = poisson, data = prussian, offset = log(year))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6340  -1.1051  -0.8166   0.5278   2.0496  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.660e+00  2.500e-01 -18.640   <2e-16 ***
## corpI        3.224e-09  3.535e-01   0.000   1.0000    
## corpII      -2.877e-01  3.819e-01  -0.753   0.4512    
## corpIII     -2.877e-01  3.819e-01  -0.753   0.4512    
## corpIV      -6.931e-01  4.330e-01  -1.601   0.1094    
## corpIX      -2.076e-01  3.734e-01  -0.556   0.5781    
## corpV       -3.747e-01  3.917e-01  -0.957   0.3387    
## corpVI       6.062e-02  3.483e-01   0.174   0.8618    
## corpVII     -2.877e-01  3.819e-01  -0.753   0.4512    
## corpVIII    -8.267e-01  4.532e-01  -1.824   0.0681 .  
## corpX       -6.454e-02  3.594e-01  -0.180   0.8575    
## corpXI       4.463e-01  3.202e-01   1.394   0.1633    
## corpXIV      4.055e-01  3.227e-01   1.256   0.2090    
## corpXV      -6.931e-01  4.330e-01  -1.601   0.1094    
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 321.10  on 279  degrees of freedom
## Residual deviance: 294.97  on 266  degrees of freedom
## AIC: 628.05
## 
## Number of Fisher Scoring iterations: 5

(Intercept) は基準となる軍団 G の解析結果である。ここで、軍団 G について、実際の観測値とポアソン回帰により予測した値をプロットすると以下のようになる。

library(MASS)
par(oma = c(0, 0, 0, 2))
predicted <- dpois(0:5, lambda = 2.231e-01)
observed <- prussian[prussian[, "corp"] == "G", 1]
truehist(observed, xlim = c(-1, 6), h = 1, prob = F, ylab = "observed", xlab = "", x0 = -0.5, col = "grey", axes = F)
axis(4)
par(new = T)
plot(0:5, predicted, xlim = c(-1, 6), ylab = "", xlab = "fatalities", pch = 20, cex = 2, axes = F)
axis(2)
mtext("fitted",side = 4, line = 3)
axis(1)

次に、もう 1 例を取り上げる。corpXI の結果について見てみる。

library(MASS)
par(oma = c(0, 0, 0, 2))
predicted <- dpois(0:5, lambda = 2.745975)
observed <- prussian[prussian[, "corp"] == "XI", 1]
truehist(observed, xlim = c(-1, 6), h = 1, prob = F, ylab = "observed", xlab = "", x0 = -0.5, col = "grey", axes = F)
axis(4)
par(new = T)
plot(0:5, predicted, xlim = c(-1, 6), ylab = "", xlab = "fatalities", pch = 20, cex = 2, axes = F)
axis(2)
mtext("fitted",side = 4, line = 3)
axis(1)

References

  1. 久保 拓弥 確率と情報の科学 データ解析のための統計モデリング入門 ― 一般化線形モデル・階層ベイズモデル・MCMC. 2012.
  2. Dobson AJ. An Introduction to Generalized Linear Models. Second Edition. 2002.
  3. Carolyn AJ. Poisson Regression or Regression of Counts (& Rates). PDF
  4. AIDS/STI -related database Japan. HIV患者とAIDS患者報告数の年次推移. AIDS/STI