疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

Rでポアソン回帰分析

今回はRで学ぶデータサイエンスシリーズ『カテゴリカルデータ分析』の第7章ポアソン回帰分析のついてまとめる。
(超基本かと存じます)

はじめに

ポアソン分布というのは交通事故に代表されるように、その事象が発生する確率が極めて小さい事象に関する分布である。
このポアソン分布に従うデータの特徴とその場合のパラメータの推定法を学ぶ。

次に示す表はある市の2015年1月の脳梗塞による救急搬送の数を示している。

件数 0件 1件 2件 3件 4件 5件 6件 7件以上
日数 8 7 5 5 3 2 1 0

この結果を見ると、0件が最も多く7件以上起こった日は0である。
ポアソン分布のモデルは次のような過程で導かれるモデルである。

  1. 時間を細分化すると、各時間帯で発生しているイベントは1回だけである。
  2. 細かく分けた時間帯でのイベントの発生する確率は同じである。
  3. 他の時間帯のイベントの発生状況の影響を受けない。

ポアソン分布では、イベントの発生回数を {Y}

{ \displaystyle
P(Y = k)= \frac{\lambda^k}{k!} e^(-\lambda)       (k = 0,1,2, .... )
}

を仮定する。
λは短い時間帯でのイベントの発生確率に対するパラメータで、λの値を変えることでポアソン分布の形は変わる。
λを大きくすることで分布は右に移動し、ばらつきもそれに応じて大きくなる。


f:id:ryosukefujii0320:20160225154748p:plain
(画像出典:統計用語辞典http://aoki2.si.gunma-u.ac.jp/

パラメータの推定に関しては、最尤推定を用いる。
それに伴う、95%CIは以下のように求める。

{ \displaystyle
r - 2\sqrt{\frac{r}{m}} \leq \lambda \leq r + 2\sqrt{\frac{r}{m}}
}

ただし、この数式は近似的であることに注意する。

Rで信頼区間

t <- c(0:7)
x <- c(8,7,5,5,3,2,1,0)
m <- sum(x)
s <- t %*% x
r <- s/m
r #推定値
         [,1]
[1,] 1.935484
ci <- c(r- 2* sqrt(r/m), r+ 2* sqrt(r/m))
ci #95%信頼区間
[1] 1.435744 2.435224

よってλの推定値は1.94 (1.44 - 2.43)であることがわかる。

ここで求めたパラメータを使用して、期待値を求めてみる。

さらにその期待値と先ほどの観測値を棒グラフと折れ線グラフで描画する。

x <- c(8,7,5,5,3,2,1,0)
names(x) <- c(0:7)
y <- dpois(c(0:7), lambda=r) * sum(x)
barplot(x, ylim=c(0, 10))
lines(y)

f:id:ryosukefujii0320:20160423022134p:plain

実際にデータ解析

今回は『データ解析のための統計モデリング入門: 一般化線形モデル・階層ベイズモデル・MCMC』(いわゆる緑本)のdata3a.csv(こちらのwebページをご参照ください)を使用する。

data <- read.csv("data3a.csv", sep=",", header=T)
head(data)

#データの中身はこんな感じです
   y     x f
1  6  8.31 C
2  6  9.44 C
3  6  9.50 C
4 12  9.07 C
5 10 10.16 C
6  4  8.32 C

ここではyが種子数、xが植物の体のサイズ、fは肥料の種類ということになっています。
データを図示すると、、、

plot(data$x, data$y, pch = c(21, 19)[data$f])
legend("topleft", legend = c("C", "T"), pch = c(21, 19))

f:id:ryosukefujii0320:20160423024321p:plain

ポアソン回帰の場合はglmを指定します。(glmはもちろん一般化線形回帰モデル(Generalized Linear regression Model: GLM)ですね)

glm(y ~ x + f, data=data, family="poisson")
Call:  glm(formula = y ~ x + f, family = "poisson", data = data)

Coefficients:
(Intercept)            x           fT  
    1.26311      0.08007     -0.03200  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.81 	AIC: 476.6

この結果を見ると、このように解釈できます。
(久保先生の本より抜粋)
f:id:ryosukefujii0320:20160423024644p:plain


こんな感じで簡単に分析自体はできます。


参考文献

1. 久保拓弥著『データ解析のための統計モデリング入門: 一般化線形モデル・階層ベイズモデル・MCMC
www.amazon.co.jp



2. 藤井良宜『Rで学ぶデータサイエンス1 カテゴリカルデータ解析』
www.amazon.co.jp



3. 粕谷英一著『Rで学ぶデータサイエンス10 一般化線形モデル』
www.amazon.co.jp



20160423
RF