個体群生態学の分野では、採取してきた生物の全長を使って、年級群というグループに分けるコホート解析を行うことが定石となっています。コホート解析では、主に混合正規分布から複数の正規分布を単離して、一つの正規分布を一つのコホートと見なします。
今回は、このコホート解析の肝となる混合正規分布の分解を行う R のパッケージ mixtools を見つけたので、その使い方を整理します (詳細は CRAN の PDF へ)。
以下のコードでは、一旦自分で混合正規分布になるデータを作成し、次に関数 normalmixEM を使って単一の正規分布に分解します。分解したあと、精度確認ってほどじゃないけど、確認のために、自分で作った混合正規分布データの各正規分布のパラメーターと推定されたパラメーターを並べて表示します。
なお、本来ならパラメーター推定のために初期値を入力する必要があるのですが、今回は関数 locator を使ってクリックで指定するようにしました。
#作成日時 : 14:52 2015/02/20
#作成者名：T.Nakano
rm(list=ls())
gc();gc()
library(mixtools)
# /////////////////////////////////////////////////////////////////////////////////
# テストデータの準備 --------------------------------------------------------------
# 混合正規分布のデータを作成する関数 ------------------------------------------
fun_CNDdata <- function(para){
dtwr <- NULL
for(i in 1:nrow(para)){
paraDT <- rnorm(para[i,1], para[i,2], para[i,3])
dtwr <- c(dtwr, paraDT)
}
dtwr
}
# 混合正規分布データの準備 -------------------------------------------------------
# パラメーターの設定 -----
DT <- PARA <- NULL
Ni <- 1
for(i in 1:Ni){
# 個体数 -----
mn <- 20
mx <- 500
exdt <- mn+runif(20)*(mx-mn) # mn から mx までの乱数作成
actn <- as.vector(summary(exdt)[c(1,2,3,5,6)]) # 5つのデータを抽出
actn <- round(actn[rev(order(actn))])
# 平均殻長 -----
mn <- 5
mx <- 50
exdt <- mn+runif(20)*(mx-mn)
actμ <- as.vector(summary(exdt)[c(1,2,3,5,6)])
actμ <- round(actμ[order(actμ)])
# 標準偏差 -----
mn <- 1.5
mx <- 3
exdt <- mn+runif(20)*(mx-mn)
actσ <- as.vector(summary(exdt)[c(1,2,3,5,6)])
actσ <- round(actσ[order(actσ)], 2)
# データフレームを作成 -----
eval(parse(text=paste("para", i, " <- data.frame(n=actn, μ=actμ, σ=actσ)", sep="")))
eval(parse(text=paste("dtwr", i, " <- fun_CNDdata(para", i, ")", sep="")))
eval(parse(text=paste("DT <- rbind(DT, data.frame(date=", i, ", TL=dtwr", i, "))", sep="")))
eval(parse(text=paste("PARA <- rbind(PARA, data.frame(date='Day", i, "', para", i, "))", sep="")))
}
# DT に混合正規分布のデータが入っている。
# View(DT)
# 混合正規分布の分解 -------------------------------------------------------------
# ヒストグラムを書く -----
pitch=ceiling((max(DT$TL)-min(DT$TL))/45) # だいたいバーが45本
BKS <- seq(-0.5, ceiling(max(DT$TL))+pitch*2+0.5, by=pitch)
TLhist <- hist(DT$TL, breaks=BKS, plot=F)
hist(DT$TL, breaks=BKS, ylim=c(0, max(TLhist$counts)*1.5))
# 初期値を設定する -----
# ピークと思われるところを全てクリックする。
# 右クリック → 停止 で locator() を止める。
int <- locator()
# 正規分布を単離する -----
NM <- normalmixEM(DT$TL, mu=int$x)
xdt <- seq(min(DT$TL), max(DT$TL), length=200)
for(j in 1:length(int$x)){
lines(xdt, NM$lambda[j]*dnorm(xdt, mean=NM$mu[j], sd=NM$sigma[j])*length(DT$TL)*pitch, col=2)
}
# 精度確認 ----------------------------------------------------------------------
# 平均 -----
data.frame(設定値=PARA$μ, 推定値=round(NM$mu))
# 標準偏差 -----
data.frame(設定値=PARA$σ, 推定値=round(NM$sigma, 2))
# 個体数 -----
data.frame(設定値=PARA$n, 推定値=round(NM$lambda*length(DT$TL)))
上記プログラムを実行すると、以下の様なグラフが表示されます。
この時点で既に locator が発動しているので、(windows では) マウスカーソルが + になっています。
次に、この状態でピークと思われるところをクリックしていきます。クリックしたところの値を locator が取得し、関数 normalmixEM の初期値とします。今回は 5 つの正規分布を混ぜてあるので、ピークと思われるところは 5 箇所あると思います。
全てクリックし終えたところで、右クリック → 停止 で locator() を止めます。すると、以降のコードが実行されて、単離された正規分布の図と、そのパラメーターの値が表示されます。
※ただし、今回のプログラムでは、混合正規分布のデータを作成する際に乱数を発生させているので、この画像の通り再現されるとは限りません。
