分布の集中度を調べる

ランダムか、集中か

水田の調査圃場において、正方形に区切った調査区画毎にニカメイガの卵塊が次のように見つかったとする。

egg <- c(2,4,1,1,
         4,0,3,3,
         4,3,0,1,
         3,2,1,1,
         2,1,4,3,
         2,2,1,0,
         1,0,5,1,
         2,1,3,2,
         2,2,2,2,
         4,1,1,2,
         3,2,0,1,
         1,2,1,0)

もし、ニカメイガの産卵が先に卵塊があるか否かに影響されず、ランダムに行われるのであれば、その分布は二項分布に従うことが予想される。
また、今回の例のように、任意の個体の特定の調査区への飛来確率pが非常に小さく*1、総個体数Nが大きい場合、二項分布をポアソン分布で近似できる。
よって、得られた分布のポアソン分布への適合度を確認することで分布がランダムか否かを確認できる。
また、ポアソン分布は分散と平均値が等しいという性質を持つことから、得られた分布の分散と平均の比\sigma^2/m(実際にはs^2/\bar{x})が1に等しいか、大きいか、小さいかを確認することでもポアソン分布かどうかの判断ができる。

ポアソン分布への適合度を確認する

まず、得られた分布から度数分布表を作成する。そして、分布がポアソン分布に従っていた場合の各度数の理論値を求める。ポアソン分布はパラメータとして母平均mのみを持つ。ここでは平均値\bar{x}を母平均の不偏推定量として用いる。

> ## ランダム分布か、集中分布か(ポアソン分布への当てはめ)
> # 度数分布の計算
> (e.table <- table(egg))
egg
 0  1  2  3  4  5 
 6 15 14  7  5  1 
> # ポアソン分布における各度数の比率計算
> (e.table <- data.frame(count = c(e.table, 0, 0),
+ p = dpois(0:7, mean(egg))))
  count           p
1     6 0.156583374
2    15 0.290331673
3    14 0.269161656
4     7 0.166356857
5     5 0.077113335
6     1 0.028596195
7     0 0.008837019
8     0 0.002340758
> # 最後の項はまとめる
> e.table$p[8] <- 1 - sum(e.table$p[1:7])

比率は合計が1になるように、分布に現れなかった項まで計算し、余分な分はまとめておく。

> # 区画数の理論値計算
> (e.table <- cbind(e.table, phi = e.table[,2] * length(egg)))
  count           p        phi
1     6 0.156583374  7.5160020
2    15 0.290331673 13.9359203
3    14 0.269161656 12.9197595
4     7 0.166356857  7.9851291
5     5 0.077113335  3.7014401
6     1 0.028596195  1.3726174
7     0 0.008837019  0.4241769
8     0 0.003019892  0.1449548
> # 理論値5以下の区画はまとめる
> (e.table <- rbind(e.table[1:4,], lapply(e.table[5:8,], sum)))
  count         p       phi
1     6 0.1565834  7.516002
2    15 0.2903317 13.935920
3    14 0.2691617 12.919759
4     7 0.1663569  7.985129
5     6 0.1175664  5.643189

ここで理論値が小さな項目をまとめるのは、観測値は0と1の間をとらないため、後にカイ二乗統計量を計算する際にこれが大きくなりすぎるのを防ぐため、ということのようである。
次にカイ二乗統計量=観測値と理論値の差の自乗の総和を求める。

> # 理論値との差
> (e.table <- cbind(e.table, d = abs(e.table$count - e.table$phi)))
  count         p       phi         d
1     6 0.1565834  7.516002 1.5160020
2    15 0.2903317 13.935920 1.0640797
3    14 0.2691617 12.919759 1.0802405
4     7 0.1663569  7.985129 0.9851291
5     6 0.1175664  5.643189 0.3568109
> # 差の自乗
> (e.table <- cbind(e.table, d2 = e.table$d^2))
  count         p       phi         d        d2
1     6 0.1565834  7.516002 1.5160020 2.2982620
2    15 0.2903317 13.935920 1.0640797 1.1322655
3    14 0.2691617 12.919759 1.0802405 1.1669196
4     7 0.1663569  7.985129 0.9851291 0.9704794
5     6 0.1175664  5.643189 0.3568109 0.1273140
> # カイ二乗統計量の計算
> (e.table <- cbind(e.table, "d2/phi" = e.table$d2/e.table$phi))
  count         p       phi         d        d2     d2/phi
1     6 0.1565834  7.516002 1.5160020 2.2982620 0.30578251
2    15 0.2903317 13.935920 1.0640797 1.1322655 0.08124799
3    14 0.2691617 12.919759 1.0802405 1.1669196 0.09032054
4     7 0.1663569  7.985129 0.9851291 0.9704794 0.12153584
5     6 0.1175664  5.643189 0.3568109 0.1273140 0.02256065
> (chi <- sum(e.table$"d2/phi"))
[1] 0.6214475

そして得られたカイ二乗統計量を検定にかける。自由度は階級数-2=5となる。

> # 検定(自由度=階級数-2)
> pchisq(chi, 3, lower.tail = FALSE)
[1] 0.8915054

以上より分布がポアソン分布に従わないとは言えない。
ニカメイガの卵塊はポアソン分布に近い分布をするだろう(実際にはプロットしてみたほうが良い)。

分布の集中度を調べる

上述のように分布がポアソン分布に従うのであれば、s^2/\bar{x}を調べることで分布の集中度、一様度が判定出来る。s^2/\bar{x}

  • 1より小さい = 一様分布
  • 1である = ランダム分布(ポアソン分布)
  • 1より大きい = 集中分布

である。
さらに、s^2/\bar{x}はF分布するので、F検定によって分布が集中分布であるか否かを検定できる。自由度はn_1=n-1, n_2=\inftyを用いる。
s^2/\bar{x}が1より小さい場合は、\bar{x}/s^2について検定を行うことで一様分布であるか否かの検定をする。

> ## バリアンスと平均値の関係を確認する
> var(egg)/mean(egg)
[1] 0.8489123
> qf(0.95, length(egg)-1, Inf)
[1] 1.361726
> pf(var(egg)/mean(egg), length(egg)-1, Inf, lower.tail=FALSE)
[1] 0.7590387
> pf(mean(egg)/var(egg), length(egg)-1, Inf, lower.tail=FALSE)
[1] 0.188295

ここで注意が必要なのは、s^2/\bar{x}は平均値によって値が変わってしまうため、s^2/\bar{x}が大きいほど分布が集中している、とは言えないという点である。
この点を改良した指数として、I_\delta指数(Morisita's indexとも)がある。この指数はs^2/\bar{x}と同様にポアソン分布で1、一様分布で1より小、集中分布で1より大となるが、平均値に影響されず*2、平均値の異なる集団同士の集中度の比較に用いる事ができる。
I_\delta=\frac{\sum x_i(x_i-1)}{N(N-1)}
ここでx_ii番目の区画における個体数、Nは総個体数である。

> ## I_δ示数の利用
> i.delta <- function(x){
+ length(x) * (sum(x * (x - 1))/(sum(x)*(sum(x)-1)))
+ }
> i.delta(egg)
[1] 0.9193054

*1:二項分布の平均値はNpに等しいことから求めると、約2%

*2:一部例外はある