ポアソン分布や二項分布は制約が厳しい

　これまでポアソン分布や二項分布を仮定した一般化線形モデルを考えてきた。これらの確率分布を使うことで、分散が刻々と変化するデータであっても無理ない当てはめを実現することが可能になってきた。しかしながら、これらの分布ですらまだ対応できない状況が存在する。

　ポアソン分布は平均と分散が等しい分布である。しかし、現実のデータはしばしばよくばらつくので、平均 = 分散を満たすとは限らないだろう。むしろ、大体は平均 < 分散となることが多い。二項分布も同様に確率とサイズパラメータが決まれば、分散は一意に決まるが、これも満たされるとは限らず、やはり分散のほうが大きくなりがちだったりする。このような、予想される分布の分散よりも実測値の分散が大きいことを、過分散overdispersionと呼ぶ。正規分布などを仮定するときは、通常、過分散なる状況は存在しない。なぜなら、分散はパラメータとして推定されるからだ。

　本稿では、もしデータに過分散があるときに、それを考慮せずに解析をしたらどうなるかをシミュレーションしてみよう。

過分散データのシミュレーション

　では、実際にデータを扱ってみよう。まずは、xに対するyの値を図示する。

------------------------------------------------------

library(plotn)

x <- c(23.8, 9.7, 24.7, 27.5, 22.7, 17.4, 2.6, 10.6, 27.8, 

      13.7, 11.7, 19.7, 7.8, 2.5, 27.6, 13.5, 10.2, 27.1,

      13.4, 9.9, 9.8, 14.7, 20.4, 8.3, 11.4, 19.9, 17.1, 

      15.7, 19, 4.9, 28.8, 23.6, 20.8, 26.2, 16, 20, 9, 

      12.3, 0.6, 26.9, 19.7, 26.2, 7.8, 24, 9.5, 21.7, 

       18.4, 22.7, 10.2, 11.9, 2.9, 10, 10.2, 22.1, 8.1,

      7.4, 1.6, 18.1, 12.6, 12.3)

y <- c(9, 5, 5, 6, 1, 1, 4, 18, 5, 5, 6, 11, 2, 8, 4, 11,

       6, 9, 15, 6, 8, 10, 9, 4, 5, 3, 4, 1, 4, 5, 8, 10,

       6, 4, 16, 4, 9, 4, 4, 8, 14, 5, 8, 5, 2, 8, 1, 7,

       12, 4, 7, 10, 7, 6, 6, 9, 8, 15, 5, 5)

plotn(y ~ x)#図1の描画

------------------------------------------------------

図の様子からxとyの間には明確な関係はなさそうである。ここは一度、xに対してyは関係がないと考えて、yの平均と分散を計算してみる。

------------------------------------------------------

mean(y)

## [1] 6.783333

var(y)

## [1] 14.1387

------------------------------------------------------

すると、yは平均よりも分散が大きく、ポアソン分布を仮定すると過分散の状況にある。実はこのデータ、負の二項分布と呼ばれるポアソン分布よりも大きな分散を示す確率分布から生成されている。今回は、下記のような設定で生成した。

------------------------------------------------------

b0 <- 2

b1 <- 0

size <- 6

n <- 60

x <- runif(n, 0, 30)

y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

------------------------------------------------------

上記の設定で、平均がexp(b0 + b1x)に従う負の二項分布のデータを生成できる。b0 = 0としてデータを生成しているので、xとは関係のないデータになっているわけである。とりあえず、ここでは負の二項分布でデータを生成するとポアソン分布から考えて過分散なデータができるととらえてほしい。

　では、このようなデータに関して、ポアソン分布を仮定して一般化線形モデルの下で解析したらどうなるだろうか。併せて回帰曲線も描くことにする。

------------------------------------------------------

fit1 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.7806  -1.1218  -0.3255   0.7759   3.5388  

## 

## Coefficients:

##              Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  1.939725   0.111772  17.354   <2e-16 ***

## x           -0.001640   0.006526  -0.251    0.802    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 120.94  on 59  degrees of freedom

## Residual deviance: 120.88  on 58  degrees of freedom

## AIC: 341.69

## 

## Number of Fisher Scoring iterations: 4

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y ~ x) #図2の図示

xx <- seq(min(x), max(x), length = 200)

yy <- exp(b0e + b1e * xx)

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

今回は、切片項 = 1.94、xの係数 = -0.0016でおおむね良い推定値になっている。では、これを10000回繰り返してシミュレーションする。過分散の様子を確認してもらうため、b0 = 2として上記は設定したが、ここでは一度b0 = 0としてからシミュレーションする。

------------------------------------------------------

b0 <- 0

b1 <- 0

size <- 6

n <- 60

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

  d <- data.frame(x = x, y = y)

  fit1 <- glm(y ~ x, data = d, family = poisson(link = "log")) #一般化線形モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  b0e <- c(b0e, coef(fit2)[1,1])

  b1e <- c(b1e, coef(fit2)[2,1])

}

histn(p1, xlab = "P value", ylab = "頻度") #図3の図示

histn(p2, xlab = "P value", ylab = "頻度") #図4の図示

sum(p1 < 0.05)

## [1] 681

sum(p2 < 0.05)

## [1] 670

------------------------------------------------------

設定したパラメータはともに0だから、今回は危険率を見ていることになる。すると、ともに危険率は7%弱で、5%を超過している。つまり、過分散データに無理やりポアソン分布を当てはめると第1種の過誤を犯す可能性が高まるということになる。このときの推定パラメータの期待値は以下の通り。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度") #図5の図示

histn(b1e, xlab = "b1", ylab = "頻度") #図6の図示

mean(b0e)

## [1] -0.01821896

mean(b1e)

## [1] -0.0001071337

------------------------------------------------------

平均的に見れば、ともに0に近いので、推定値自体はうまく推定されることがわかる。問題は、有意な関係がないのに過剰に検出されやすい、ということである。続いて、切片を0ではない値にしてみる。

------------------------------------------------------

b0 <- 2

b1 <- 0

size <- 6

n <- 60

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

  d <- data.frame(x = x, y = y)

  fit1 <- glm(y ~ x, data = d, family = poisson(link = "log")) #一般化線形モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  b0e <- c(b0e, coef(fit2)[1,1])

  b1e <- c(b1e, coef(fit2)[2,1])

}

histn(p1, xlab = "P value", ylab = "頻度") #図7の図示

histn(p2, xlab = "P value", ylab = "頻度") #図8の図示

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 1840

------------------------------------------------------

すると、xの係数に関する危険率は18%となり、非常に大きくなった。過分散データでは、このようにxとの有意な関係がないときでも、過剰に有意な関係を見出す原因になりうる。このときの推定パラメータの期待値は以下の通り。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度") #図9の図示

histn(b1e, xlab = "b1", ylab = "頻度") #図10の図示

mean(b0e)

## [1] 1.993428

mean(b1e)

## [1] 7.14713e-05

------------------------------------------------------

平均的に見れば、推定値自体はうまく推定されることがわかる。

　最後にxの係数が0ではない場合を考えてみる。以下のようなデータがあるとする。図示も行っておく。

------------------------------------------------------

x <- c(2.4, 29.1, 6.1, 1.7, 19.5, 16.2, 28.9, 17.2,

       2, 16.1, 4.8, 27.5, 22.8, 29.2, 3.8, 21.2, 

       16.4, 29.9, 7.9, 23.1, 8.5, 2.5, 28.9, 24.2, 

       12.2, 13.5, 4.3, 7.3, 8, 25.2, 11.9, 25.9, 21.9,

       13.6, 16.5, 12, 6.9, 19.8, 5, 12.8, 12.1, 23.2, 

       13.1, 7.8, 20, 20.9, 25.7, 9.7, 12.7, 16.7, 11,

       18.8, 28.6, 25.2, 29.7, 9.2, 29.1, 21.5, 27.4, 4.6)

y <- c(1, 5, 1, 1, 0, 2, 3, 3, 1, 1, 4, 2, 1, 14, 1,

       1, 3, 6, 5, 8, 1, 1, 6, 2, 2, 2, 1, 3, 2, 2, 

       3, 2, 5, 0, 3, 1, 0, 6, 1, 6, 1, 11, 1, 3, 5,

       8, 7, 2, 4, 3, 4, 4, 1, 4, 5, 0, 7, 8, 10, 1)

plotn(y ~ x) #図11の図示

------------------------------------------------------

実際問題としては、このデータだけをみてポアソン分布を仮定すると過分散であることを見抜くのは難しい。とりあえず、ポアソン分布を仮定して解析しよう。

------------------------------------------------------

fit1 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.6758  -0.9915  -0.2875   0.6828   2.6816  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) 0.169403   0.190490   0.889    0.374    

## x           0.056716   0.008703   6.517 7.18e-11 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 140.610  on 59  degrees of freedom

## Residual deviance:  94.947  on 58  degrees of freedom

## AIC: 261.09

## 

## Number of Fisher Scoring iterations: 5

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y ~ x) #図12の図示

xx <- seq(min(x), max(x), length = 200)

yy <- exp(b0e + b1e * xx)

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

ここで、過分散の指標として使えるものが2つあって、一つが、残差逸脱度residual decianceをその自由度で割ったものである。fit2オブジェクトの中に格納されているデータを使って計算できる。過去にならって、残差逸脱度と自由度の定義からの計算もしてみよう。

------------------------------------------------------

fit2$deviance/fit2$df.residual

## [1] 1.637019

d <- 2*(log(dpois(x = y, lambda = y)) - log(dpois(x = y, lambda = exp(b0e + b1e*x))))

sum(d)/(length(y) - 2)

## [1] 1.637019

------------------------------------------------------

もう一つ、より安定な方法がdispersion parameter Φを計算することである。Φはピアソンのχ^2と呼ばれる値を、自由度で割ったものである。

ピアソンのχ^2はポアソン回帰においては残差平方を期待値で割ったものの和であり、いわゆる検定統計量のχ^2とほとんど同じものである。定義からわかるように、要はピアソンのχ^2は期待値からのばらつきを測る指標である。

　Φの主張は以下のように式変形するとわかりやすい。つまり、不偏分散（回帰曲線からのずれの推定値）は、モデルから予測されるずれと一致する、つまりそれらの比は1になるだろう、というものである。

もしその比が1にならないのであれば、それはモデルが不適切、当てはまりが悪いものを無理やり当てはめている、と言えるわけだ。Φ > 1であれば、モデルに対して実測値のばらつきは大きい、過分散である。一方で、Φ < 1なら、それは予測されるよりも実測値のばらつきが小さい、過少分散underdispersionである。

　手計算と、関数を使って計算してみよう。

------------------------------------------------------

sum(((y - exp(b0e + b1e*x))^2)/(exp(b0e + b1e*x)))/(length(y) - 2)

## [1] 1.563703

sum(residuals(fit1, type = "pearson")^2)/(length(y) - 2)

## [1] 1.563703

------------------------------------------------------

2つの指標は、ともにポアソン分布（あるいは二項分布）を仮定すれば1に近い値となり、1を超えれば過分散である。もちろん、ポアソン分布から生成された値でも、ぴったり1になることはまずないので、1.5以上になった時は過分散を想定して解析するべきだと言われる。上記の2つの指標はともに1.5を超えてしまっているので、データはポアソン分布に対して過分散である。

　2つの指標の存在は、暗に残差逸脱度は漸近的にχ^2に一致するだろうことを言っている。実際、残差逸脱度の定義式をテイラー展開すると、漸近的に一致することを示すことができる。

　実は、今回のデータは下記のようにb0 = 0、b1 = 0.07として、負の二項分布からデータを生成していた。試しに同じ期待値の下で、ポアソン分布からデータを生成して描いてみよう。

------------------------------------------------------

b0 <- 0

b1 <- 0.07

size <- 6

n <- 60

rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

yp <- rpois(n = n, lambda = exp(b0 + b1*x))

yp <- c(2, 6, 1, 0, 6, 1, 8, 1, 0, 1, 2, 7, 6, 7, 2, 

        8, 5, 6, 5, 8, 1, 1, 9, 7, 2, 2, 0, 6, 2, 7,

        3, 7, 3, 1, 4, 2, 1, 1, 0, 1, 2, 6, 3, 2, 4, 

        3, 8, 2, 2, 1, 2, 6, 7, 2, 9, 1, 9, 7, 4, 1)

#筆者の手元で生成されたデータの一例

plotn(y ~ x) #図13の図示

overdraw(points(x, yp, col = "red"))

------------------------------------------------------

確かに、ポアソン分布に比べると今回のデータはよりばらつきが大きく見える。

　ではこの状況での危険率を調べよう。下記のようにシミュレーションする。

------------------------------------------------------

b0 <- 0

b1 <- 0.07

size <- 6

n <- 60

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

  d <- data.frame(x = x, y = y)

  fit1 <- glm(y ~ x, data = d, family = poisson(link = "log")) #一般化線形モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  b0e <- c(b0e, coef(fit2)[1,1])

  b1e <- c(b1e, coef(fit2)[2,1])

}

histn(p1, xlab = "P value", ylab = "頻度") #図14の図示

histn(p2, xlab = "P value", ylab = "頻度") #図15の図示

sum(p1 < 0.05)

## [1] 962

sum(p2 < 0.05)

## [1] 10000

------------------------------------------------------

すると、切片係数に関する危険率は10%弱となり、非常に大きくなった。このときの推定パラメータの期待値は以下の通り。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度") #図16の図示

histn(b1e, xlab = "b1", ylab = "頻度") #図17の図示

mean(b0e)

## [1] -0.01314853

mean(b1e)

## [1] 0.07024124

------------------------------------------------------

やはり、推定値自体は近い値が出てくるようだ。

おわりに

　以上のように、過分散が与える影響について議論してきた。基本的に、過分散データへの無理な当てはめによって、第1種の過誤を犯す確率が上がると考えられる。

　そもそも過分散はなぜ生じるのだろうか。いくつかの理由が考えられるが、重要な要素として、説明変数の中に、さらにその下の階層構造が存在する場合が考えられる。その場合は、その階層構造を考慮したうえでの解析、一般化線形混合モデルの出番となる。あるいは、一般化線形混合モデルから自然に導かれる発想として、過分散に対応できる確率分布を使う、という手段も考えられる。今後は、一般化線形混合モデルに着目しながら、解析を行っていこう。