R で混合正規分布の分解 ～連続でコホート解析する～

前回の 2 回までで、混合正規分布から正規分布を単離する方法をご紹介しました。
これらのコードは調査 1 回分のデータを解析するものでした。

しかし、通常、個体群の調査は 1 回で終わることなく、同じフィールドで繰り返し行われます。ということは、データも調査 1 回分ではなく、複数回分あるはずです。

前回のコードで解析しようとすると、データをわざわざ Excel で 1 回分だけ取り出して、csv で保存して、R で解析して、また別の 1 回分のデータを Excel 取り出して、csv で保存して、R で解析して、また……、csv で……、解析して……。

面倒くさいこと極まりないですよね。

しかも、1 回 1 回の解析を別々で行うということは、コホートに ID を付けずに解析することになるので、成長や生残を追跡することができません。成長や生残が追跡できないということは、もはやコホート解析とは呼べないんじゃないかとも思います。



ということで、今回は、前回のコードを拡張し、データを全部ひっくるめて読み込んで、1 回分のデータを抽出するところから R にやってもらい、そしてコホートに ID を付けていくことで、後々の成長・生残の解析ができる形に持って行こうというコードを作りました。


さらに、このコードではデータファイルがあるフォルダに、正規分布が割り当てられたグラフが PDF 形式で書き出され、それぞれのパラメーターの値が csv 形式で書き出されるようにしておきました。


テストデータをこちらに置いておきます。
今回のコードの実行用に使用するなり、自分のデータを整形する参考にするなり、ご自由にお使い下さい。

#作成日時 : 15:33 2015/05/12
#作成者名：T.Nakano
rm(list=ls())
gc();gc()

library(mixtools)
library(tcltk)
# /////////////////////////////////////////////////////////////////////////////////

#ダイアログボックスから文字を記入する関数 ----------------------------------------
library(tcltk)
modalDialog <- function(title,question,entryInit,entryWidth=30,returnValOnCancel="ID_CANCEL"){
  dlg <- tktoplevel()
  tkwm.deiconify(dlg)
  tkgrab.set(dlg)
  tkfocus(dlg)
  tkwm.title(dlg,title)
  textEntryVarTcl <- tclVar(paste(entryInit))
  textEntryWidget <- tkentry(dlg,width=paste(entryWidth),textvariable=textEntryVarTcl)
  tkgrid(tklabel(dlg,text="       "))
  tkgrid(tklabel(dlg,text=question),textEntryWidget)
  tkgrid(tklabel(dlg,text="       "))
  ReturnVal <- returnValOnCancel
  onOK <- function()
  {
    ReturnVal <<- tclvalue(textEntryVarTcl)
    tkgrab.release(dlg)
    tkdestroy(dlg)
  }
  OK.but     <-tkbutton(dlg,text="   OK   ",command=onOK)
  tkgrid(OK.but)
  tkgrid(tklabel(dlg,text="    "))
  tkfocus(dlg)
  tkbind(textEntryWidget, "<Return>", onOK)
  tkwait.window(dlg)

  return(ReturnVal)

}

# データの読み込み --------------------------------------------------------------
# データが整理されている csv ファイルを選択する。
library(tcltk)
fileName <- tclvalue(tkgetOpenFile())
act1 <- unlist(strsplit(fileName,"/"))
act1 <- act1[-length(act1)]

directory <- NULL
if(length(act1)!=1){
for(i in 1:(length(act1)-1)){
directory <- paste(directory,act1[i],"/",sep="")
}
}
directory <- paste(directory,act1[length(act1)],sep="")
setwd(directory)
WR <- read.csv(fileName, header=T)
WR$TL <- as.numeric(as.character(WR$TL))

# コホート解析開始 ---------------------------------------------------------------
PARA <- NULL
unidate <- unique(WR$date)
for(i in 1:length(unidate)){
DT <- WR[which(WR$date==unidate[i]),]

for(j in 1:200){
# ヒストグラムを書く -------------------------------------------------------------
win.graph()
#pitch=(max(WR$TL)-min(WR$TL))/15     # だいたいバーが15本
pitch=1
BKS <- seq(-pitch, max(WR$TL)+pitch*2, by=pitch)

TLhist <- hist(DT$TL, breaks=BKS, plot=F)
hist(DT$TL, breaks=BKS, main=unidate[i], ylim=c(0,max(TLhist$counts)*1.5))

# 初期値を設定する ---------------------------------------------------------------
# ピークと思われるところを全てクリックする。
# 右クリック → 停止 で locator() を止める。
int <- locator()

# 正規分布を単離する -------------------------------------------------------------
if(length(int$x)==1){
NM <- data.frame(mu=mean(DT$TL), sigma=sd(DT$TL), lambda=1)
xdt <- seq(min(DT$TL), max(DT$TL), length=200)
lines(xdt, NM$lambda*dnorm(xdt, mean=NM$mu, sd=NM$sigma)*length(DT$TL)*pitch, col=2)

}else if(length(int$x) > 1){
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)
}
}    # if(length(int$x)==1)

cond <- tkmessageBox(message="コホート分けをやり直しますか？",icon="question",type="yesno",default="yes")
if(as.character(cond)=="no") break
}     # for(j in 1:200)

PARA1 <- data.frame(mu=NM$mu, sigma=NM$sigma, N=NM$lambda*length(DT$TL), lambda=NM$lambda)

maxdt <- PARA1$lambda*dnorm(PARA1$mu, mean=PARA1$mu, sd=PARA1$sigma)*length(DT$TL)*pitch
text(PARA1$mu, maxdt, labels=paste("age", 1:length(maxdt), sep=""), pos=3)

CN <- NULL
for(j in 1:length(maxdt)){
ReturnVal <- modalDialog("cohort name", paste("age", j, " → ") ,"")
CN <- c(CN, ReturnVal)
}

hist(DT$TL, breaks=BKS, main=unidate[i], ylim=c(0,max(TLhist$counts)*1.5))
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)
}
text(PARA1$mu, maxdt, labels=CN, pos=3)

PARA1 <- data.frame(cohort=CN, date=unidate[i], PARA1)

PARA <- rbind(PARA, PARA1)
}     # for(i in 1:length(unidate))

View(PARA)

# 結果を保存する ----------------------------------------------------
write.csv(PARA, file="混合正規分布の分解_結果_パラメーター.csv", row.names=F)

pdf("混合正規分布の分解_結果_ヒストグラム.pdf")
mxv <- NULL
for(i in 1:length(unidate)){
DT <- WR[which(WR$date==unidate[i]),]
BKS <- seq(-pitch, max(WR$TL)+pitch*2, by=pitch)
TLhist <- hist(DT$TL, breaks=BKS, plot=F)
mxv <- c(mxv, max(TLhist$counts))
}

BKS <- seq(-pitch, max(WR$TL)+pitch*2, by=pitch)

for(i in 1:length(unidate)){
DT <- WR[which(WR$date==unidate[i]),]
hist(DT$TL, breaks=BKS, main="", ylim=c(0,max(mxv)*1.5))
legtxt <- paste(unidate[i], "\n", "N = ", nrow(DT) ,sep="")
legend("topleft", legend=legtxt, bty="n", cex=0.8)

PARA1 <- PARA[which(PARA$date==unidate[i]),]
xdt <- seq(-pitch, ceiling(max(WR$TL))+pitch*2+0.5, length=200)
for(j in 1:nrow(PARA)){
lines(xdt, PARA1$lambda[j]*dnorm(xdt, mean=PARA1$mu[j], sd=PARA1$sigma[j])*length(DT$TL)*pitch, col=2)
}

maxdt2 <- NULL
for(j in 1:nrow(PARA)){
act1DT <- max(PARA1$lambda[j]*dnorm(xdt, mean=PARA1$mu[j], sd=PARA1$sigma[j])*length(DT$TL)*pitch)
maxdt2 <- c(maxdt2, act1DT)
}
text(PARA1$mu, maxdt2, labels=PARA1$cohort, pos=3)

}

dev.off()

​

このコードを実行すると、まずデータファイルを指定するよう求められます。

データファイルをしてすると、ファイルが読み込まれ、1 回目の調査のデータが抽出されます。

そして、前回同様、locator が発動し、ピークをクリックして初期値を与えるよう求められます。

全てクリックし終えたところで、右クリック → 停止 で locator() を止めます。

すると今回は、「コホート分けをやり直しますか？」と R に質問されます。

「はい」か「いいえ」で答えてあげて下さい。

「はい」の場合、もう一度ピークをクリックし直せます。

「いいえ」の場合、それ以降のコードが実行されます。

すると、コホートの名前を入力するよう求められます。

コホート名を入力するウィンドウに書いてある age1 というのは、図中の age1 に対応しています。コホートの数だけコホート名を入力するよう求めてくるので、全てに名前を付けてあげて下さい。(これがコホートの ID になります。)

ここまでが調査 1 回分の解析の手順です。

次からは 2 回目の調査のデータが抽出され、上記内容が繰り返されます。

2 回目以降は、前回の単離結果を見ながらピークをクリックしていくのが良いと思ったので、前回の結果を残しておくことにしました。

以下のようにグラフをずらすと前回の結果を見ることができます。

ちなみに前のグラフには、命名したコホート名がグラフ内に記されています。

前回の結果を確認しながら今回のピークを判定しクリックして下さい。

全てクリックし終えたら locator を止めて下さい。

（ここでも、やり直すかと聞かれるので答えてあげて下さい。）

次に、前回のコホート名を確認しながら今回のコホート名を入力して下さい。
ここで、2 回目以降は、前回の判定を確認しながら、同じコホートだと判断できるものには同じコホート名を入力して下さい。この時、大文字と小文字は区別されるので注意してください。

ここまでくれば、あとは同じことを繰り返していくだけです。

そして、全ての作業が終了すると、PDF と csv ファイルが作成されます。


[修正]
2015/10/07
先日、ヒストグラムが上手く書けないというご指摘をいただいたので、その部分を変更しました。

2015/10/15
日付と個体数をグラフ左上に表示できるように変更しました。

2016/02/10
データを読み込んだとき、TL が Factor という型になってしまい、正常に解析できないという問題が見つかりました。
これについて、
WR <- read.csv(fileName, header=T)
の下に、
WR$TL <- as.numeric(as.character(WR$TL))
を入れることで対処しました。

2016/02/10
混合正規分布でない場合、すなわち単一の正規分布で構成されている場合でも解析できるように変更しました。