（2024.03.27追記）より強力なRパッケージ"mixR"が利用可能になっているようです。正規分布だけでなく、lognormal, Gamma, Weibull分布まで対応している模様（追記ここまで）
# ４部構成にしています。こちらのページはより基本形で正規混合モデルを使用しました。
# その１：正規混合モデル（オーソドックスな基本形）
# その２：ガンマ混合モデル（体長の増加に伴い分散も増加するモデル）
# その３：前後関係からユーザがグループ数を仮定して正規混合モデルを実行
# その４：体長と体幅など複数種類のデータで分ける（有効な最終手段）
*******************************************************
（注：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を混ぜたサンプル
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))
注：mclust では、BIC = 2*対数尤度 になっているようで、一般的な BIC の定義と正負が逆である。このため mclust の結果では BIC 最大のモデルがベストと判断。
E：各グループで分散は共通とするモデル、V：分散は異なるとするモデル。
グループ分け数 2 で最大 → 2 分割するのが妥当である。グループ数の推定は成功した。しかし、いずれのモデルでもグループ分けは2のモデルがベストとされたが、分散が等しいモデルが選ばれてしまうため、冒頭の仮定にしたがい分散が異なるモデルでパラメータを推定した。
cf. 標準偏差が倍くらいの違いがない限り分散が等しいモデルが選択されてしまった。
Mc <- densityMclust(Length2, modelNames="V")
# "V":分散が異なるモデル（cf. Mclustでもいいが、densityMclustの方が図示に便利）
# "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")
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 の部分は何かの数字とかに置き換えないと動作しない）。
0 件のコメント:
コメントを投稿