（このページは 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で比較
その１同様に、Mclustの推定結果を初期値に用いた。じつはガンマ混合分布も扱えるより汎用性の高いパッケージもあるのだが推定があまりよくなかったため使用しなかった（ページの末を参照）。Mclustは正規混合分布なので分布型が違うのだが、あくまで初期値としての使用なので構わないだろう。
library(mclust)
plot(mclustBIC(LengthG))
３グループ、分散共通のモデルがベストであるという解析結果が得られた。
mc <- summary(densityMclust(LengthG), parameters=T) # 推定されたパラメータの一覧
もし、前の時点の解析結果などから、今回のグループ数推定がおかしいと判断される場合は、グループ数をユーザ側で指定してもよいだろう。G=3のようにグループ数を指定できる。
# mc <- summary(densityMclust(LengthG, G=3), parameters=T)
Mclustの結果を初期値に用いて、ベイズのガンマ混合モデル推定をした。
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)
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
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倍に増やして再試行したほうがよいだろう。
なお、いずれかの 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) # 平均と標準偏差
0 件のコメント:
コメントを投稿