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

（このページは 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で十分だと思う