過小分散なカウントデータを扱いたい
背景
カウントデータをモデリングしようと思ったとき、まず思い浮かべる分布といえばポアソン分布だと思います。しかしポアソン分布は期待値と分散が等しいという性質があるため、実際のデータに当てはめようとすると、(ポアソン分布から期待されるものよりも)分散が大きい(過分散)または小さい(過小分散)という問題に直面することがあります。
そのようなとき、過分散ならともかく1、過小分散の場合にどんな分布を当てはめれば良いのか知らなかったのですが、先日某所でダブルポアソン分布というものを教えてもらったので試してみます。
doublepoisson
を触ってみる
いきなりですが、ダブルポアソン分布で調べてみたところRではrmutil
パッケージでdoublepoisson
関数が使えるようなので早速インストールして使ってみます。
install.packages("rmutil") library(rmutil)
rmutil
では、runif
やdnorm
などと同じようにr/d/p/qを頭に付けたdoublepois
関数が使えます。ダブルポアソン分布に従う乱数を生成してみましょう。Overdispersion parameterとしてs
が指定できるのでひとまず小さめの値にしてみます。
set.seed(123) n <- 500 m <- 5 # 平均 s <- 2 # Overdispersion parameter d <- data.frame(y = rdoublepois(n, m, s))
> head(d) y 1 4 2 6 3 5 4 7 5 8 6 3
ぱっと見た感じでは通常のポアソンと変わりありません。このときの平均と分散は:
> mean(d$y) [1] 5.016 > var(d$y) [1] 2.452649
でした。うーん、分散がs
と結構違ってますが、これ合っているんですかね。。。
ひとまず次に進めます。rpois
と比べ分散が小さいデータを生成できているか確認してみましょう。カウントごとの頻度を描き2、その上に(通常の)ポアソン分布の確率関数を引いてみます。
## プロット用に最大値を得る max_count <- max(d$y) ## 平均がmのポアソン分布における確率密度を得る nrm_poi <- dpois(0:max_count, m) ## データの分布を得る dens <- table(d$y) / nrow(d) ## プロット用に最大密度を得る。まるめの単位はよしなに。 max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 ## 欠損に対応するために0からmax_countまでのベクトルを用意して値を代入する dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens ## ヒストグラムではなく棒グラフで作図するので、座標を指定する xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) ## 棒グラフに確率密度関数を載せる plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "") rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue")
通常のポアソン分布で期待される密度(青線)と比較するとダブルポアソン分布では平均値周りの密度が高くなっており、全体としてバラツキが小さくなっているようです。
ではs
の値を変えるとどうなるでしょうか?簡単に変更できるように関数化しておきます。
plot_double_poisson <- function(n, m, s) { d <- data.frame(y = rdoublepois(n, m, s)) max_count <- max(d$y) nrm_poi <- dpois(0:max_count, m) dens <- table(d$y) / nrow(d) max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "", main = paste0("s = ", s), cex.main = 0.7) rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue") }
s
を1, 3, 5, 10とした場合のグラフを描いてみます。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 5, 1) plot_double_poisson(500, 5, 3) plot_double_poisson(500, 5, 5) plot_double_poisson(500, 5, 10)
おや? s
が小さい方がデータのバラツキが大きくなっていますね。s
が10のときにはデータの半数以上が5となっているようです。またs
が1のときにはほぼポアソン分布と一致しているように見えます3。この傾向はm
を変更しても変わらないようです。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 10, 1) plot_double_poisson(500, 10, 3) plot_double_poisson(500, 10, 5) plot_double_poisson(500, 10, 10)
想像していた結果とは異なるのですが、このOverdispersion parameterを変更することでポアソン分布では扱えなかった過小分散なカウントデータに対応できそうです。
切片を推定してみる
glm
で推定してみる
それでは続いて、この過小分散なデータに対してGLMを当てはめたときの結果を確認してみましょう。改めてサンプルデータを作成しますが、十分に過小分散となるようs
は10としておきましょう。
set.seed(123) n <- 500 m <- 5 s <- 10 d <- data.frame(y = rdoublepois(n, m, s))
glm
を使ってフィッティングしますが、ここで切片として5が推定されれば問題ありません。またリンク関数はlog
とします。
res_glm <- glm(y ~ 1, d, family = poisson("log"))
> coef(res_glm) (Intercept) 1.609838
切片は1.609
と表示されていますが、逆リンク関数exp
で戻してみると、
> exp(coef(res_glm)) (Intercept) 5.002 > mean(d$y) [1] 5.002
一致しています。過小分散であっても切片は正しく推定されるようです。しかしポアソン分布であると仮定した場合、その期待値と分散は等しい(E(Y) = Var(Y))はずなのですが、実際の分散を見てみると:
> var(d$y) [1] 0.4709379
大きく異なります。
dglm
で推定してみる
それではこのようなデータに対してどのようにモデリングすれば良いかということなのですが、調べてみたところdglm
というパッケージ・関数があるようなので使ってみましょう。
install.packages("dglm") library(dglm)
> res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log")) stats::glm.fit(x = Z, y = d, weights = rep(1/2, N), offset = doffset, でエラー: 'x' の中に NA/NaN/Inf があります
あらら、何らかのエラーで止まってしまいました。原因を探りたいところですが、ひとまずs
を小さくして先に進みます。
set.seed(123) n <- 500 m <- 5 s <- 3 d <- data.frame(y = rdoublepois(n, m, s))
res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log"))
今度は大丈夫なようですので結果を見てみましょう。
> exp(res_dglm$coefficients[1]) (Intercept) 5.01 > mean(d$y) [1] 5.01
今度もm
は推定できているようです。分散はどうでしょうか。
> exp(res_dglm$dispersion.fit$coefficients) (Intercept) 0.3372774 > var(d$y) [1] 1.685271
なんか、やっぱり分散は違いますね。。。
実はここまで触れてこなかったのですが、ダブルポアソン分布についてはこちらのブログを参考にさせて頂きました:
この記事ではダブルポアソン分布の分散として
と紹介されていますので、もしかしたら分散はm
との積を求める必要があるかもしれません。
> sqrt(exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients)) (Intercept) 1.299908 > var(d$y) [1] 1.685271
さっきよりもだいぶ近い値が得られるようになりましたが、右辺が平方根になっているのが気になります(分散なので)。取ってみましょう。
> exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients) (Intercept) 1.68976 > var(d$y) [1] 1.685271
だいぶ近い!…のですが、分散の定義…。
回帰係数を推定してみる
glm
、dglm
で推定してみる
なかなか良さそうな感触が得られたので、切片だけでなく回帰係数も推定してみましょう。以下のようにサンプルデータを作成します。s
は3としました。
set.seed(123) n <- 500 b <- 1.5 x <- runif(n, 3, 5) m <- b * x s <- 3 y <- rdoublepois(n, m, s) d <- data.frame(y, x)
平均m
を与えたときのダブルポアソン分布を見てみましょう。
hist(d$y, col = "gray", ylim = c(0, 0.5), probability = T, main = "", xlab = "", breaks = "scott")
また平均と分散を確認すると:
> mean(d$y) [1] 5.982 > var(d$y) [1] 2.751178
このようになりました。
それぞれのデータを、glm
とdglm
それぞれでフィッティングしてみます。まずはglm
から。
res_glm <- glm(y ~ x, d, family = poisson("log"))
> exp(coef(res_glm)) (Intercept) x 2.057136 1.302974
回帰係数は真の値(1.5)にまぁまぁ近いものが推定されています。続いてdglm
でフィッティングすると:
res_dglm <- dglm(y ~ x, ~1, data = d, family = poisson("log"))
> exp(res_dglm$coefficients) (Intercept) x 2.057136 1.302974
glm
と同じ結果が得られました。
では分散を見てみましょう。m
にはd$x
の平均値を用いることにします。
m_hat <- exp(res_dglm$coefficients %*% c(1, mean(d$x)))
> m_hat [,1] [1,] 5.914556 > mean(m) [1] 5.985851
> exp(res_dglm$dispersion.fit$coefficients) * m_hat [,1] [1,] 1.947131
うーん、合いません(真の値は3)。。。
ちなみにこのデータをlm
でフィッティングするとどうなるかと言うと:
res_lm <- lm(y ~ x, d)
> var(res_lm$residuals) [1] 1.935673
あまり変わりませんね。
optim
で推定してみる
もう少しチャレンジしてみましょう。dglm
の代わりにoptim
で推定してみます。目的関数を以下のように定義します:
my_obj_fun <- function(b) { b1 <- b[1] b2 <- b[2] b3 <- b[3] ret <- sum(ddoublepois(d$y, exp(b1 + d$x * b2), b3, log = TRUE)) ret }
パラメータのベクトルを適当に用意して、optim
にかけてみましょう。
b <- c(0, 1, 3) res_opt <- optim(b, my_obj_fun, control = list(fnscale = -1))
結果を確認します:
> exp(res_opt$par)[-3] # 回帰係数 [1] 2.050437 1.303756 > res_opt$par[3] # Overdispersion [1] 3.064928
おぉ!回帰係数はglm
、dglm
の結果と大体一致しています。またs
の推定値は3.06となっており、これも真の値と近しいものになっています。optim
いいですね。
終わりに
今回は過小分散なカウントデータを扱うためにrmutil
とdglm
というパッケージを使ってみました。optim
による推定は思いの外うまくいったようですが、dglm
の結果は少しパッとしませんでした。ただしs
の求め方については本当にこれで合っているのか自信がないため、推定できているかについての判断は保留です。ご存知の方がいらっしゃれば教えてください。