混合分布のパラメータ推定：その４（複数種類のデータで分ける）

（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています）

# ４部構成にしてあります、ここではグループ数の推定がうまくいかない場合の対処について。
# より基本形の例は、その１の方でまとめてあります。
# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）

*******************************************************

その１〜３では、ひとつの基準（体長とか、個体サイズを表す情報）を元にした分け方を紹介した。しかし、ここまでの方法だけではどうしても分けられないケースは出てくるだろう。

ところで、複数の基準…体長と体幅とか、体長と頭幅とか、体長と体重、のように１サンプルにつき複数の情報が得られるならば、より高い精度でコホート分離が可能である。こういうデータをお持ちの方 or 取ることが可能な方はぜひ試したらよいと思う。

そもそも、ここまでで紹介した mclust は本来はそういう複数の基準で混合分布を解くのが基本のパッケージである。基準同士の相関関係も利用して分離してくれるので強力である。
なお、多変量の混合分布なので、多変量正規混合分布（multivariate normal mixture distribution）のモデリングとなる。



### まずはサンプルデータの設定
N <- 100
mean1.2 <- c(10, 1.0)
mean2.2 <- c(13, 1.3)
var1.2 <- matrix(c(1.2, 1, 1, 1.2), 2, 2)^2 
var2.2 <- matrix(c(1.4, 1, 1, 1.4), 2, 2)^2

### 母集団からのサンプリング、N1, N2 ずつ取ってきて混ぜた
ratio <- c(0.8, 0.2) # 8 : 2 の比率でサンプリングした
# つまりトータル100で80、20個体のコホートが混じっている
 N1 <- N*ratio[1]
 N2 <- N*ratio[2]

set.seed(1) # 乱数のタネを指定
library(MASS)
Length2 <- rbind(mvrnorm(N1, mean1.2, var1.2), mvrnorm(N2, mean2.2, var2.2))  # 1, 2を混ぜたサンプルを作る

# （体長と体幅のようなサイズ関係）
par(mfrow=c(1,2))
hist(Length2[,1], col="gray")
hist(Length2[,2], col="gray")

# ２峰であることは想像がつくが、ひどい混ざり方である。このケースではどちらか一方のみを用いた解析では分離できなかった。

# 基本的なやり方はここまでと同じだが、

# データは2基準のデータがセットであり2列になるようにする（上記のLength2のように）

# BICのモデルはより複雑になる

library(mclust) 

plot(mclustBIC(Length2)) 

BICが最大となっているモデルのグループ分け数は2となり、ちゃんと分けられた。

EEVなどのモデルはコホートの二次元分布のサイズ、形状、方向の仮定の組み合わせで表されている（詳しくは、citation("mclust") で出てくるFraley et al. 2012のp.8のテーブルを参照）。

# 推定値の見方もこれまでと同じ

Mc <- densityMclust(Length2)

mc <- summary(Mc, parameters=T) 

mc

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

Density estimation via Gaussian finite mixture modeling 

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

Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components:

 log.likelihood   n df       BIC    ICL

      -309.8967 100  9 -661.2399 -664.1

Clustering table:

 1  2 

86 14    # 正解は80と20なので少しグループ１の方が多目の判定

Mixing probabilities:

        1         2 

0.8566939 0.1433061 

Means:

          [,1]      [,2]

[1,] 10.304959 14.246681  # 値の大きい方の基準は10、13が正解

[2,]  1.126679  1.721643    # 小さい方の基準は1.0、1.3が正解なので、ちょっと粗い

Variances:

[,,1]

          [,1]      [,2]

[1,] 1.7080561 0.9243498    # 1.7..と0.96...は２つの基準のそれぞれ分散

[2,] 0.9243498 0.9646183    # 0.92..の２つは共分散

[,,2]

          [,1]      [,2]

[1,] 1.1210113 0.9727447

[2,] 0.9727447 1.5516630

# 推定は多少粗いものの、このデータでこれだけ分けられれば十分だろう。

# 図示の仕方も同じ、出てくるものは標高図になる

plot(Mc, data=Length2)

横軸、縦軸に２つの基準を取り、二次元のヒストグラムを標高図にしている

# ベイズ版は省略。複数基準ならMclustで十分だと思う

混合分布のパラメータ推定：その３（グループ数をユーザ側で指定する場合のコード…）

（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています）

# ４部構成にしてあります、ここではグループ数の推定がうまくいかない場合の対処について。
# より基本形の例は、その１、２の方でまとめてあります。
# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）

*******************************************************

（その１からの続きです）
体サイズ頻度分布のヒストグラムでよくあることですが、高齢のコホートが他と比べてかなり少ない場合には、解析すると隣のグループに含められてしまうことが多々ある（想定するグループ数と異なると感じる場合…それは主観的かもしれないが）。


もし、前の時点でのグループ数などの事前情報によって、そこに確かにコホートが存在するはずという仮定ができるのならば、グループ数は既知であるとして解析してもよいだろう。完全に仮定無しでの解析はなかなか難しい…混合分布問題は完全解決には至っていないと思う。


その２でやり方を書いたような、グループ数指定（G=◯の項を加える）をしても、途中にもう一つグループを作られて、肝心の最高齢集団が無視されたりとなかなか言うことを聞いてくれない場合がある。

その２のガンマ混合モデルならば解決できるだろうが、ベイズを使うのは敷居が高いだろうし計算時間もかかる。ここでは mclust（正規混合モデル）を用いた別の解決方法を紹介する。

mclust （正規混合モデル）を使用して、グループ数をユーザ側で決めた場合の混合分布モデリングを行う。一番若い齢（体長などの値の小さい）のコホートから順に推定し、推定し終わったら除いて残りのデータで解析するというステップ・バイ・ステップ法で処理することによって、集団サイズの小さい最高齢のコホートも検出することができる。

### 母集団に用いるサンプルデータ

N <- 200 # サンプル全体のN数
mean1 <- 10
mean2 <- 15
mean3 <- 20
sd1 <- 1.0
sd2 <- 1.2
sd3 <- 1.4

### 母集団からのサンプリング、N1, N2, N3 ずつ取ってきて混ぜた
ratio <- c(0.5, 0.3, 0.2) # 5:3:2 の比率でサンプリングした
 N1 <- N*ratio[1]
 N2 <- N*ratio[2]
 N3 <- N*ratio[3]

set.seed(1) # 乱数のタネを指定
Length3 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2), rnorm(N3, mean3, sd3)) # 1, 2, 3を混ぜたサンプル

### mclust による正規混合モデリングを最若齢〜最高齢のコホートへと順に１つずつ推定するプログラム（いったん読み込ませたら、Rを終了するまで何度でも再利用できる）

step.mclust <- function(data, n.group, var="V", priors=NULL) {
require(mclust)
x.min <- 5*floor(0.2*min(data))
x.max <- 5*ceiling(0.2*max(data))
x <- seq(x.min, x.max, by=(x.max - x.min)/100)
summ <- numeric(0)
for (g in 1:n.group-1) {
mc <- summary(Mclust(data, G=n.group - g, modelNames=var), parameters=T)
n <- round(mc$pro[1]*length(data)) # グループ１の推定個体数
n2 <- round(mc$pro[2]*length(data)) # グループ２の推定個体数
m <- mc$mean[1] # グループ１の推定平均値
m2 <- mc$mean[2]
s <- sqrt(mc$variance[1]) # グループ１の推定標準偏差
s2 <- sqrt(mc$variance[2])
# １の中では確率密度が小さく、かつ２の方が確率密度が大きいものは２である
# １の中では確率密度が小さいが、２よりも確率密度が大きいものも１である
d <- dnorm(data, m, s) # dataがグループ１に属する確率密度
d2 <- dnorm(data, m2, s2) # 同、グループ２に…
dd <- d > d2 # グループ２よりもグループ１に属する確率が高い個体
p.rank <- rank(1/d) # 確率密度の高い順
gr <- data[dd ==T & p.rank <= n] # グループ１に属する可能性が高い個体（でもまだグループ１が残ってる）
# 面倒でも二段階抽出をしないとグループ１の左の裾が落ちてしまう
gr2 <- setdiff(data, gr) # dataからgrを除いた残り（この中にグループ1が残っている）
r.remain <- rank(1/dnorm(gr2, m, s)) # 残り物の中で、グループ１に属する確率が高いものの順
remain <- gr2[r.remain <= n - length(gr)] # グループ１の残党
gr <- c(gr, remain) # remainをグループ１に追加
data <- setdiff(gr2, remain) # グループ２からremainを除外したものを次のループのdataにする
summ <- rbind(summ, data.frame(n, m, s)) }
rownames(summ) <- 1:n.group
return(summ) }

### 使用法
step.mclust(Length3, 3) # データ、グループ数を指定するだけ

### 出力結果（n：各グループの推定 N数、m：同 平均値、s：同 標準偏差）

#      n        m         s
# 1  99 10.09472 0.8822265
# 2  60 14.68063 1.4391196
# 3  41 19.91919 2.3426576

# 最後のグループの標準偏差の推定が良くない。
# また、もっと最後のグループの集団サイズを減らすとけっこう推定は悪くなる…。とくに集団サイズの推定がよくない。
そういうケースにもきっちりと対応したい場合には、その２のガンマ混合モデルでがんばるなど試行錯誤するしかないだろう。

混合分布のパラメータ推定：その２（ガンマ混合モデル）

（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています）

# ４部構成にしてあります、こちらはガンマ混合モデルを使用の例。
# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）

*******************************************************

（その１からの続きです）

体サイズ頻度分布のヒストグラムでは、複数のコホートがオーバラップすることは珍しくなく、実用的には多数の重なったグループの分離はできてほしい。

まずは、その１で使用したプログラムを試したところ、平均の推定はできた。しかし、標準偏差の推定はうまくいかず、とくに真ん中のグループが狭められてしまう傾向にあった（最後のところに推定の図を載せておく）。

（注：推定の粗さが気にならなければ、"その１"の推定だけで十分だと思います）

解決策としては、目的に応じてモデルに制約を掛けることが挙げられる。ひとつの方法は分散がグループ間で共通という仮定することである。これならばその１で使用したプログラムを微修正するだけで十分に推定が可能だった（mclust()にmodelNames="E"の項を加える）。

今回はコホートの分離が目的なので、体サイズの成長に伴いコホートの分散は増加するという仮定を加えてもよいだろう。それならばガンマ分布を使用すれば、平均の増加に伴いバラツキが増加するのでグループ毎に変更するパラメータを１つ減らすことができる。

その１と同様に、ガンマ分布の母集団を作り、そこから抽出したサンプルデータから元の母集団のパラメータを推定できるか確認し、手法の有効性を確かめた。

ガンマ分布ではrate、shapeといったパラメータを使用するが、平均や標準偏差とは以下のような対応関係にある。

# 母集団の設定

N <- 200 # サンプル全体のN数

mean1 <- 10

mean2 <- 15

mean3 <- 20

sd <- 1.2 # 標準偏差は共通

ratio <- c(0.5, 0.3, 0.2) # 5 : 3 : 2 の比率で混ぜる

         N1 <- N*ratio[1]

         N2 <- N*ratio[2]

         N3  <- N*ratio[3]

# rate, shapeパラメータへ変換

        rate1 <- mean1 / sd^2

        rate2 <- mean2 / sd^2

        rate3 <- mean3 / sd^2

        shape1 <- mean1 * rate1

        shape2 <- mean2 * rate2

        shape3 <- mean3 * rate3

# shapeは平均 * rateとなっており、このため標準偏差をグループ間で共通にしても平均に比例して値が増加する。

set.seed(1) # 乱数のタネを指定

LengthG <- c(rgamma(N1, rate=rate1, shape=shape1), rgamma(N2, rate=rate2, shape=shape2), 

rgamma(N3, rate=rate3, shape=shape3)) # ガンマ分布３グループを混ぜる

嫌な感じに混ざった…標準偏差は共通にしたにもかかわらず、山は次第に幅が増加するのがガンマ分布の性質

その１同様に、Mclustの推定結果を初期値に用いた。じつはガンマ混合分布も扱えるより汎用性の高いパッケージもあるのだが推定があまりよくなかったため使用しなかった（ページの末を参照）。Mclustは正規混合分布なので分布型が違うのだが、あくまで初期値としての使用なので構わないだろう。


library(mclust)
plot(mclustBIC(LengthG)) 

# グループ分け数のモデルをBICで比較

３グループ、分散共通のモデルがベストであるという解析結果が得られた。

mc <- summary(densityMclust(LengthG), parameters=T) # 推定されたパラメータの一覧

もし、前の時点の解析結果などから、今回のグループ数推定がおかしいと判断される場合は、グループ数をユーザ側で指定してもよいだろう。G=3のようにグループ数を指定できる。
# mc <- summary(densityMclust(LengthG, G=3), parameters=T)

Mclustの結果を初期値に用いて、ベイズのガンマ混合モデル推定をした。

n.g <- 3 # グループ数の設定
data <- list(n.g=n.g, N=length(LengthG),  size=LengthG,  # 上記のデータ LengthG を使用
            alpha=mc$pro, # mclustの推定した混合比率を初期値に使用
            M=mc$mean) # mclustの推定した平均値を初期値に使用
par <- list(ratio=mc$pro, # mclustの推定した混合比率を初期値に使用
           tau=1/mean(mc$variance),F, # mclustの推定した分散の平均を初期値に使用
           m.group=mc$mean, # mclustの推定した平均値を使用
           group=sample(1:n.g, length(LengthG), replace=T),F, # 初期値不明のためランダム
           sigma=NA, delta.m=NA, shape=NA, rate=NA)

model <- function() {

    for (i in 1 : N) { # iループ

size[i] ~ dgamma(shape[group[i]], rate[group[i]]) # 混合したサイズ分布（グループ毎に変化）

group[i] ~ dcat(ratio[1:n.g]) # グループ選択の事前分布、カテゴリカル分布を使用

    } # iループ終わり

ratio[1:n.g] ~ ddirch(alpha[]) # 混合比率の事前分布：ディリクレ分布=ratioの合計が1になる

tau ~ dgamma(1, 0.01) # 分散の逆数、グループ間で共通

for (g in 1:n.g) { # gループ（n of group）

rate[g] <- m.group[g]*tau # rateパラメータ。tauは分散の逆数であることに注意

shape[g] <- m.group[g]*rate[g] # shapeパラメータ（平均 * rate）

m.group[g] ~ dnorm(M[g], 0.5)  # 各グループの平均値の事前分布

} # gループ終わり

for (gg in 2:n.g) { delta.m[gg] <- m.group[gg] - m.group[gg-1] } # ggループ

sigma <- sqrt(1/tau) # 標準偏差の事後分布の計算

} # model終わり

source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行


# 実行結果

#               mean     sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
# ratio[1]     0.515  0.037   0.444   0.489   0.514   0.540   0.590 1.001  3000
# ratio[2]     0.289  0.036   0.221   0.264   0.289   0.313   0.361 1.002  1100
# ratio[3]     0.196  0.029   0.141   0.176   0.196   0.215   0.256 1.001  2600
# m.group[1]  10.141  0.145   9.853  10.040  10.140  10.240  10.430 1.001  3000
# m.group[2]  15.021  0.232  14.590  14.870  15.020  15.170  15.500 1.001  3000
# m.group[3]  19.984  0.242  19.510  19.820  19.980  20.150  20.470 1.003   970
# sigma        1.326  0.085   1.179   1.266   1.319   1.378   1.510 1.001  3000
# delta.m[2]   4.881  0.245   4.402   4.719   4.875   5.045   5.353 1.001  3000
# delta.m[3]   4.963  0.299   4.364   4.761   4.969   5.165   5.537 1.001  3000
# shape[1]    59.152  7.140  45.460  54.287  59.025  63.972  73.270 1.001  3000
# shape[2]   129.826 15.976  98.769 119.000 129.650 140.500 161.500 1.001  3000
# shape[3]   229.892 28.998 174.195 210.200 229.550 249.500 288.800 1.001  3000
# rate[1]      5.835  0.714   4.493   5.347   5.825   6.320   7.244 1.001  3000
# rate[2]      8.644  1.066   6.618   7.918   8.631   9.370  10.750 1.001  3000
# rate[3]     11.503  1.440   8.727  10.540  11.480  12.480  14.420 1.001  3000
# deviance   680.517 15.944 655.697 669.200 678.400 689.400 717.600 1.002  1500

母集団のratioの値は0.5, 0.3, 0.2に対し、推定値（中央値）は0.514, 0.289, 0.196。平均は10, 15, 20に対し、10.140, 15.020, 19.980。標準偏差は1.2に対し、1.319。標準偏差がやや大きく推定されたものの、他の推定はかなり良い結果だった。
なお、いずれかの R.hatの値が 1 から大きく離れていたり、n.eff 値が1~2桁に落ちている場合は推定が良くない。n.iter, n.burn, n.thinの値を比例して 5~10倍に増やして再試行したほうがよいだろう。

母集団（実線）と推定値（点線）の確率密度を図示した。標準偏差がやや大きく推定されたため山が低めになったがかなり高い推定精度だ。平均のみをグループ間で変化させたにもかかわらず、ちゃんとバラツキが平均に伴い増加するパターンをモデル化することができた。

# ちなみに、上記のグラフを書くコードはこちら

hist(LengthG, breaks=seq(5*floor(0.2*min(LengthG)), 5*ceiling(0.2*max(LengthG)), by=1.25), 

ylim=c(0,40), xlab="Length (mm)", main="", lwd=1.5, col=gray(0.9))

# 棒の区画幅をものぐさするためにコードが複雑になっている…

x <- seq(5,25, 0.1) # x軸を決める

lines(x, N1*dgamma(x, rate=rate1, shape=shape1), lwd=1.5) # 実線のカーブ

lines(x, N2*dgamma(x, rate=rate2, shape=shape2), lwd=1.5)

lines(x, N3*dgamma(x, rate=rate3, shape=shape3), lwd=1.5)

n1 <- round(seg$median$ratio[1]*length(LengthG))

n2 <- round(seg$median$ratio[2]*length(LengthG))

n3 <- round(seg$median$ratio[3]*length(LengthG))

# 点線のカーブ

lines(x, n1*dgamma(x, shape=seg$median$shape[1], rate=seg$median$rate[1]), lty=2, lwd=1.5)

lines(x, n2*dgamma(x, shape=seg$median$shape[2], rate=seg$median$rate[2]), lty=2, lwd=1.5)

lines(x, n3*dgamma(x, shape=seg$median$shape[3], rate=seg$median$rate[3]), lty=2, lwd=1.5)

最後に、ちょっとトリッキーだが、３グループの正規分布からのサンプルを元に正規混合モデル（標準偏差はグループ毎に変化）で推定した結果と、同じサンプルをガンマ混合モデル（標準偏差は一定）で推定した結果を比較してみた。

平均：10, 15, 20、標準偏差：1.0, 1.2, 1.4の３つの正規分布の母集団から、計200サンプルを混合比率：0.5, 0.3, 0.2でサンプリングし、サンプルから母集団のパラメータを推定した。

### 正規混合モデルの推定結果（抜粋）

#                 mean     sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff

# ratio[1]       0.496  0.037   0.424   0.471   0.495   0.521   0.567 1.001  3000

# ratio[2]       0.272  0.040   0.194   0.246   0.272   0.298   0.349 1.002  1900

# ratio[3]       0.232  0.037   0.167   0.208   0.229   0.254   0.309 1.001  3000

# m.group[1]    10.096  0.093   9.919  10.030  10.090  10.160  10.290 1.001  3000

# m.group[2]    14.545  0.168  14.230  14.430  14.540  14.650  14.890 1.001  3000

# m.group[3]    19.748  0.413  18.840  19.540  19.790  20.010  20.400 1.002  2000

# sigma[1]       0.886  0.072   0.761   0.835   0.882   0.930   1.042 1.001  3000

# sigma[2]       0.950  0.178   0.652   0.829   0.935   1.054   1.341 1.000  3000

# sigma[3]       1.760  0.330   1.246   1.537   1.717   1.934   2.522 1.001  2300

### ガンマ混合モデルの推定結果（抜粋）

#               mean     sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff

# ratio[1]     0.506  0.037   0.433   0.481   0.505   0.530   0.579 1.001  3000

# ratio[2]     0.293  0.034   0.229   0.271   0.292   0.315   0.361 1.001  3000

# ratio[3]     0.201  0.029   0.146   0.181   0.200   0.219   0.261 1.001  3000

# m.group[1]  10.186  0.115   9.970  10.110  10.180  10.260  10.420 1.002  1200

# m.group[2]  14.834  0.166  14.510  14.720  14.840  14.940  15.160 1.001  3000

# m.group[3]  20.178  0.192  19.810  20.050  20.180  20.300  20.560 1.001  3000

# sigma        1.114  0.065   0.995   1.067   1.110   1.155   1.248 1.001  3000

比較可能な混合比率、平均で見ると、いずれもいい推定をしているものの、ガンマ混合モデルのほうがより母集団に近い値を推定をした。母集団は正規分布であるのにもかかわらずである。また正規混合モデルの標準偏差の推定はズレが大きく、とくに真ん中のグループが狭められた。なお、DICは前者が707.2、後者が666.3なので、全体としてもガンマ混合モデルのほうが当てはまりがよさそうだ（こんなケースでDICの比較は有効なのか？？）。

DICだけでなく、全体としての母集団との当てはまりの良さを図示して比較し、絵合わせでもチェックしてみた。

実線が母集団の確率密度、黒の点線が正規混合モデルによる推定、赤の点線がガンマ混合モデル。正規混合モデルはデータの重なりにより推定が乱された様子がわかる。第3グループは正規混合モデルのほうがよく当てはまっているようだが、それ以外ではガンマ混合モデルのほうが母集団のラインにより近いものとなった。推定計算にかかる時間はガンマ混合モデルのほうが長いものの、今回のように平均にともなってバラツキが大きくなるようなデータについては、元の分布型にかかわらずより無難な推定を行なってくれることが期待できる。

なお、ここまで単一のヒストグラムを想定してきたが、実際のコホート解析では当然、時系列のデータとなる。発展形として平均や標準偏差を Y[c,t] = alpha + beta*Y[c-1,t-1] のような自己回帰モデルにすると前の時点の情報を利用してもう少し推定がしやすくなるかもしれない（希望的観測？）。だんだん欲が出てきました（c：コホート、t：時間）。

cf.その１で使用したMclustよりも汎用性の高い混合モデルを扱えるパッケージを見つけた。他の分布型も扱えるという大きな利点がある。ただし、推定精度はMclustの方がよいかもしれない。大きく狂うなら、分布型が違ってもMclustの方がよいだろう。そもそもこのパッケージは複数の量的形質を組み合わせた分離を主な用途と想定していると思う。

library(flexmix) # 要インストール

mm <- FLXMRglm(でーた ~ 1, family="gaussian") # ここで分布型を指定できる

mm3 <- flexmix(. ~ ., data=data.frame(でーた), k=3, model=mm, control=list(verbose=1, nrep=5))

# k=3 は分離するグループ数、Mclustと違って一括で比較できないので、検討すべきグループ数全部について行い（k=3, 2, 1）AIC、BIC を比較する必要がある。なお、こちらでは一般的な定義通り、値の小さいモデルがよりよいモデル。

summary(mm3) # 混合比率やAIC、BIC

parameters(mm3) # 平均と標準偏差

混合分布のパラメータ推定：その１（正規混合モデル）

（2024.03.27追記）より強力なRパッケージ"mixR"が利用可能になっているようです。正規分布だけでなく、lognormal, Gamma, Weibull分布まで対応している模様（追記ここまで）

（このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものです）


# ４部構成にしています。こちらのページはより基本形で正規混合モデルを使用しました。


# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）

*******************************************************

（注：mclustでの解析の場合、citation("mclust") で手法を使用の際の引用文献が得られる。mclustは私の作成ではありません）

似通った何かを識別する手段のひとつに統計学的な分離がある。一般的な問題の中にもたくさんあるはず…例えば、外見的な何かの特徴から不良品を識別して排除するなどだろう。沢山のデータを取って比較してみるとそういった傾向が見えてくるのは興味深い。

自分の専門分野では、野外の生物の生活史研究などでよくやる体サイズ頻度のヒストグラムでコホートを分離するとか、近縁種や雌雄間で何らかの量的形質によって両者を識別するといった用途がある。そんな時に判断材料が少なく、1 つの尺度でしか比較できないようなケースでも有効な方法を考えていた。しかし、これまでに見たことのあるやり方は、統計学的にはちょっと強引な手法であったり、自分でアルゴリズムを書かないといけないような敷居が高いものだった。また、提供媒体が古くなっていたりと現代的な解析環境では使いにくくなっているようにも感じたため、Rで実行できる解析方法を紹介することにした。

以下、コホートデータを統計学的に分離し、コホートを正規分布などに当てはめた時のパラメータ推定を試みた。仮定として、サイズ頻度分布上で隣り合うコホート同士が分離されるか否かは、コホート間の平均値が統計学的に差があれば異なるコホートとみなす、コホートの分散は成長の個体差にともなって増加する（分散の差異をコホートの違いの基準には用いない）、とした。

"mclustというR パッケージでも正規混合分布を解ける"という情報提供をいただき、使用法と使い勝手を試行錯誤してみた結果、ある程度の点推定まではこれだけで行けると判断しました。当初は最初からベイズ法で解くつもりでいましたが、ベイズ法の解析についても計算負荷を抑えるためmclust の結果を初期値として計算するように書き直しました。

この手法の有効性を確認する。まず、テストデータとして、体サイズの頻度分布を想定し、予め平均値と標準偏差を指定した 2 つの正規分布を母集団を作成、その母集団から 60, 40 サンプルずつの計 100 のサンプルを抽出し混合した。それぞれ平均 10, 14、標準偏差は 1.2, 1.4 にした。

### 母集団の設定
N <- 100 # サンプル全体のN数
mean1 <- 10
mean2 <- 14
sd1 <- 1.2
sd2 <- 1.4

### 母集団からのサンプリング、N1, N2 ずつ取ってきて混ぜた
ratio <- c(0.6, 0.4) # 6 : 4 の比率でサンプリングした
N1 <- N*ratio[1]
N2 <- N*ratio[2]

set.seed(1) # 乱数のタネを指定
Length2 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2))  # 1, 2を混ぜたサンプル

嫌な感じに混じりました…。

実線のカーブはそれぞれの母集団の正規分布の確率密度。

複数の正規分布が混ざった分布を、正規混合分布（normal mixture distribution）と言う。このサンプルを解析して元の母集団のパラメータを遡り推定し、正しく母集団を再現推定できるかどうかを検証した。パラメータとして、混合比率、各正規分布の平均、標準偏差を求めた。

まず、R パッケージの mclust を使用して大まかな推定をし、その結果をベイズ推定の初期値に使用した。mclustはEMアルゴリズムによるパラメータ推定をする。

（cf. コホート解析に混合正規分布を適用するレビューがあるとの情報提供を頂きました：赤嶺 2005. 混合正規分布とEMアルゴリズム. 水産海洋研究 69: 174–183）

library(mclust) # 要インストール（cf. 追加パッケージのインストール方法）

plot(mclustBIC(Length2)) 

# グループ分け数のモデルを情報量基準 BIC で比較（図示）

注：mclust では、BIC = 2*対数尤度 になっているようで、一般的な BIC の定義と正負が逆である。このため mclust の結果では BIC 最大のモデルがベストと判断。

E：各グループで分散は共通とするモデル、V：分散は異なるとするモデル。

グループ分け数 2 で最大 → 2 分割するのが妥当である。グループ数の推定は成功した。しかし、いずれのモデルでもグループ分けは2のモデルがベストとされたが、分散が等しいモデルが選ばれてしまうため、冒頭の仮定にしたがい分散が異なるモデルでパラメータを推定した。

cf. 標準偏差が倍くらいの違いがない限り分散が等しいモデルが選択されてしまった。

Mc <- densityMclust(Length2, modelNames="V")
# "V":分散が異なるモデル（cf. Mclustでもいいが、densityMclustの方が図示に便利）

mc <- summary(Mc, parameters=T) # 推定されたパラメータの一覧を見る

mc

# ----------------------------------------------------

# Gaussian finite mixture model fitted by EM algorithm 

# ----------------------------------------------------

# Mclust V (univariate, unequal variance) model with 2 components:

#  log.likelihood   n df       BIC     ICL

#       -213.4482 100  5 -449.9223 -461.0196

# Clustering table:

#  1  2 

# 61 39 

# Mixing probabilities:

#         1         2 

# 0.6119736 0.3880264 

# Means:

#        1        2 

# 10.18958 14.18402 

# Variances: # 標準偏差でなく分散

#        1        2 

# 1.165315 1.879066 

sqrt(mc$variance) # 標準偏差に直す

#        1        2 

# 1.079497 1.370790 

推定された確率密度をヒストグラム上に一括で図示できる（見栄えは悪いが推定の確認に便利）

plot(Mc, data=Length2)

#（hist()を用いているようで、breaksの項を加えればヒストグラムの区間を操作可能）

グループ数（G= ）、推定の初期値（prior=）を指定して実行することも可能、何らかの理由で事前情報があるのならば利用していいと思う。

Mc2 <- densityMclust(Length2, modelNames="V", G=3, prior=priorControl(mean=c(15, 21, 32)))

# それでも、相対的に小さいピークは、あるはずなのに認識されないことがある。そんなときはデータの範囲を絞って推定するといい。例えば、25より大きいデータに絞るなら、Length2[Length2 > 25] 

平均値、標準偏差ともに母集団のものに近い値を推定できた。点推定でよいならばここまでで十分だろう。標準偏差の推定の粗さも許容するのならば、3グループ以上の場合でもmclustで済ませて構わないと思う。何よりもこの方法は容易で敷居が低いのがよいところだ。




**************************************************************

せっかくなのでmclustの結果を初期値として与え、さらにベイズ法で解析し、これらの値がどれくらいの範囲を取りうるのか、確率分布を求めてみた。
一般的なベイズ法のソフトウェア WinBUGS を用いた単純なコードで推定計算した。…とはいっても、多少の R と WinBUGS の経験がないと使いこなすのは難しいかもしれません。BUGSの推定がうまくいかなかった時に n.iter, n.burn, n.thin の値をどう増やしたらいいかなどは多少の勘が必要です。
コードはこちらの本を参考にした：「マルコフ連鎖モンテカルロ法（豊田秀樹、朝倉書店）」(p. 190~)。

# ベイズモデルの記述。推定計算に必要なデータ、初期値の設定。ディリクレ分布（値の合計が1になる）を2つの集団の混合比率の重み付けに利用することで正規混合モデル化した。なお重み付けをグループの選択に用いるところで、カテゴリカル分布を利用してグループ番号に変換した。

　事前分布については、今回は事前情報があるので一般的なベイズ法の場合に比べ範囲を狭めていることに注意

（注：ここでは補助ツールWBUGSを用いて実行するため通常のWinBUGS使用と比べてパラメータの記述は単純化してある）

n.g <- 2 # グループ数の指定
data <- list(n.g=n.g, N=length(Length2),  size=Length2,  # 上記のデータ Length2 を使用
            alpha=mc$pro, # mclustの推定した混合比率を使用
            M=mc$mean) # mclustの推定した平均値を使用
par <- list(ratio=mc$pro,  tau=1/mc$variance,F, # mclustの推定した混合比率、分散を使用
           m.group=mc$mean, # mclustの推定した平均値を使用
           group=sample(1:n.g, length(Length2), replace=T),F, # 初期値不明のためランダム
 sigma=NA, delta.m=NA, delta.var=NA)
           # 初期値が不要なパラメータ（他の数値計算で求まる）


### モデル式
model <- function() {
    for (i in 1 : N) { # iループ
size[i] ~ dnorm(m.group[group[i]], tau[group[i]]) # 混合したサイズ分布（グループ毎に変化）
group[i] ~ dcat(ratio[1:n.g]) # グループ選択の事前分布、カテゴリカル分布を使用
    } # iループ終わり
ratio[1:n.g] ~ ddirch(alpha[]) # 混合比率の事前分布：ディリクレ分布=ratioの合計が1になる
    for (g in 1:n.g) { # gループ（n of group）
m.group[g] ~ dnorm(M[g], 0.5)  # 各グループの平均値の事前分布
tau[g] ~ dgamma(1, 0.01) # 各グループの分散の逆数の事前分布
sigma[g] <- 1 / sqrt(tau[g]) # 同、標準偏差の事後分布の計算
} # gループ終わり
for (gg in 2:n.g) { # ggループ
delta.m[gg] <- m.group[gg] - m.group[gg-1]  # 各グループの平均値の差の事後分布の計算
delta.var[gg] <- 1/tau[gg] - 1/tau[gg-1]  # 各グループの分散値の差の事後分布の計算
} # ggループ終わり
    } # model終わり


### 計算の実行
source("WBUGS.R") # cf. WBUGS、R2WinBUGSの補助ツール
seg <- wbugs(data, par, model, n.iter=110000, n.burn=10000, n.thin=100) # 計算実行

# 実行結果

#                       mean          sd    2.5%       25%       50%       75%    97.5%  Rhat  n.eff
# ratio[1]       0.588  0.112   0.311   0.541   0.606   0.657   0.753 1.035  1500
# ratio[2]       0.412  0.112   0.247   0.343   0.394   0.459   0.689 1.042  1300
# m.group[1]    10.194  0.239   9.761  10.050  10.185  10.330  10.680 1.004  2400
# m.group[2]    14.025  0.577  12.530  13.790  14.150  14.410  14.820 1.002  3000
# sigma[1]       1.085  0.207   0.678   0.972   1.079   1.194   1.483 1.028  3000
# sigma[2]       1.488  0.380   0.977   1.222   1.402   1.673   2.383 1.014  2300
# delta.m[2]     3.831  0.509   2.431   3.636   3.932   4.153   4.526 1.004  1300
# delta.var[2]   1.137  1.544  -0.901   0.179   0.789   1.691   5.109 1.001  3000
# deviance     332.039 21.662 308.097 317.400 325.800 339.400 390.100 1.007  2000

#（注：RhatはMCMCのchain間の収束基準、n.effはMCMCの有効サンプル数で自己相関が高いと低下する。ひとつでも Rhat > 1.1 になったり、n.eff < 1000 くらいの場合は推定があやしい。とくにRhatの方はよりシビアで、Rhatの基準が満たされていない場合は、n.iter, n.burn, n.thin を5~10倍に増やして再実行する必要がある）

# 中央値も表示（要は↑の50%の値）

seg$median
$ratio
[1] 0.6058 0.3942

$m.group
[1] 10.19 14.14

$sigma
[1] 1.078 1.414

$delta.m
[1] 3.925

$delta.var
[1] 0.82195

$deviance
[1] 325.8

# 推定値の事後確率分布を図示（MCMCサンプルの頻度分布）
par(mfrow=c(3,2))
hist(seg$sims.list$ratio[,1], main="ratio1")
hist(seg$sims.list$ratio[,2], main="ratio2")
hist(seg$sims.list$m.group[,1], main="mean1")
hist(seg$sims.list$m.group[,2], main="mean2")
hist(seg$sims.list$sigma[,1], main="sd1")
hist(seg$sims.list$sigma[,2], main="sd2")

全般的に左右非対称のため median を代表値として用います。

実行結果のdelta.m（平均値差）が2.5〜97.5%まですべて同じ符号のため、2つの集団の平均値に差があると言えた。ところが、delta.var（分散差）についてはゼロをまたいでしまいました。その確率をチェック：
dv <- seg$sims.list$delta.var
length(dv[dv<0]) / length(dv) # MCMCサンプル数のうち、負の値だった個数
0.1963333 、ベイズ法でもしっかり0.05を超えていました…。ここは冒頭の仮定を踏襲します！

ratioの真の値は0.6, 0.4に対し、推定値（中央値）は0.6058, 0.3942。mclustの推定より向上した。m.group（平均値）の真の値は10, 14としたが、推定値は10.19, 14.14でmclustと同等のよい推定だった。一方、sigma（標準偏差）の真の値は1.2と1.4としたが、推定値は1.078と1.414、第2グループの推定はmclustより向上したが、前者の方は真の値より小さいままだった。

推定された混合比率、平均値、標準偏差を用いた確率密度を点線で書き加えた。元の母集団の曲線をよく再現することができた。とくに第2グループの方は母集団の曲線にかなり近いものとなった。

第1グループの方は、むしろ推定曲線のほうが当てはまりがいいように見えた。どうやらオーバラップが激しいサンプルを用いたため、第2グループの裾野が第1グループの中心付近にまで入り込み押し上げた結果、真の平均の少し右側がピークになっており、そのピークを見た目通りに真のピークと見なして推定したためだろう。それにより第2グループの裾野は少し第1グループに食われたので、第2グループは少し右にシフトしたのだろう。このあたりが推定の限界のようだ。

なお、任意の値 X がどれくらいの確率で第 1 グループに属するかは、これで求まる。
length(Length2) * 0.6058 * dnorm(X, 10.19, 1.078) # N数 x 正規分布の確率密度関数
ここで length(Length2)*0.6058 は推定された第 1 グループの N 数。X は数字列でもOK、Length2 自身を試してみてもOK（注：X の部分は何かの数字とかに置き換えないと動作しない）。

お引越し

今月で琉大でのポストが満了となり、来月からつくばの国環研に移ります。

沖縄での研究生活はとても楽しかったです。願わくばもっとフィールドに出たくはありましたが、2年間というタイムリミットを考えると、やりたいことのすべては到底できず。もしまた次に機会があるとしたら、今度はもっと腰を落ち着けていろいろなことをやれる立場で来たいものです。

沖縄ではもっぱら統計モデリングの技術を磨いていました。プロジェクトの命題の関係もありましたが、いま何が何でもしっかり身に付けたかった技術だったからというのもあります。この過程でRの使用を前提とした勉強会や研究集会など開催したのは自分自身へのフィードバックも大きく本当によき糧になりました（参加してくださった皆さま、どうもありがとうございました）。これらを通じて階層ベイズのモデリングもそこそこできるようになってきた具合なので、これからドンドン使っていく予定です。

大学学部生時代以来の二度目のつくば市ですが、昔と違って今はそれなりに人間味のあるこなれた街になりつつあるとも聞いているので楽しみではあります。近隣の皆様、どうぞよろしくお願いします。

生態学会・サンゴ–植物 企画集会やります

無事に終了しました。来場いただいた皆さま、また演者、コメンテータの皆さまありがとうございました！

########################
直前になりましたが、来る第60回日本生態学会大会（静岡）にて、標記の企画集会を開催するので、ご案内いたします。
（大会ページの本集会要旨など http://www.esj.ne.jp/meeting/abst/60/T08.html）
（今大会の案内 http://www.esj.ne.jp/meeting/60/大会案内2/）

サンゴ礁群集は生態学全般においては特異で馴染みのない存在でしょうけれども、陸上の植物群落に相対するものと位置づけると非常に興味深い類似点・相違点が見えてきます。また、サンゴならではの特性を生かした研究の数々もご紹介できると思っています。企画者である私自身もサンゴ研究の歴史が浅いため、今なお新鮮な驚きをもって本集会を楽しみにしています。陸・海いずれの研究者もご来聴を歓迎します。


企画集会 T08 3月6日 14:30-16:30 H会場

「サンゴ礁生態学：植物生態学とのアナロジー」

企画者： 熊谷直喜（琉大・熱生研）, 井口 亮（琉大・熱生研）, 向草世香（JSTさきがけ，長大・水産，琉大・熱生研）

本集会では、生態系基盤を構成する生物群として、陸上植物に相対するサンゴを捉える。それを踏まえてサンゴ礁生態学と植物生態学の類似・相違点を整理し、陸域と海域の生態系のより一般的な理解を目指す。

造礁サンゴは、動物でありながら固着性であり群体性の体制をとっている。さらに体内に共生する褐虫藻を通じて一次生産を担うことから、陸上植物と比較して議論されることが多い。また、自らが大きな物理的構造を構築することで周辺環境と異なる生態系を構成し、生物多様性と生物生産の基盤となることも植物群落との大きな共通点である。しかし、これまで両者の生態学は十分に比較対照されてはこなかった。そこで本集会では、サンゴの植物とのアナロジーに着目し、個体（群体）、集団（個体群）、群集、生態系それぞれの階層に渡って、サンゴの生態学的研究を行なっている講演者を用意した。サンゴと植物の生態学研究の比較を行うことで、陸域と海域での生態学的現象の統一的な理解を進めるための方向性について議論する。

＊コメンテータ：陸上植物の研究者お二人にお願いしています（竹中明夫・国立環境研究所、石塚航・東大総合文化）

[T08-1] はじめに：サンゴの個体・群集・生態系の基本構造について　向草世香（JSTさきがけ，長大・水産，琉大・熱生研）
[T08-2] 光合成と環境応答：光合成をおこなうことの利点と危険性　中村 崇（琉大・理）
（コメンテータ）CO2増加による生態学的応答について　井上志保里（東大・地惑）
[T08-3] 造礁サンゴの分子生態学的アプローチ：植物研究との比較　井口 亮（琉大・熱生研）
[T08-4] サンゴ群集の時空間動態モデリング：植物群集との対比　*熊谷直喜（琉大・熱生研）, 向草世香（JSTさきがけ，長大・水産，琉大・熱生研）
[T08-5] サンゴと植物：似ているが、群集と個体群の特性は大きく違う　酒井一彦（琉大・熱生研）