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

（このページは 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

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