R で混合正規分布の分解 ～自分のデータを解析～

前回はパッケージ mixtools の使い方を整理するために、テストデータを作成して、そこから正規分布を単離する方法をご紹介しました。ですが、やはり実際の全長のデータを解析してこそコホート解析と言えるのではないでしょうか？それに、このような新しいパッケージや関数に出会ったとき「自分のデータでも解析できるのか？」という点も気になりますよね。

ということで、今回は自分のデータを使って解析できるよう、コードを書き換えてみました。といっても、R で混合正規分布の分解 ～テストデータを解析～ のテストデータの準備のところを、tcltk の OpenFile を使って csv を読み込む の内容に変えただけです。

解析の手順は R で混合正規分布の分解 ～テストデータを解析～ と一緒なのでそちらを御覧ください。

まぁ、「解析」なんて大それた言い方してますけど、ただ単に「ピークと思わしきところをクリックするだけ」なんですよねぇ…。これを「解析」とか言ってしまうと、いろんな人から怒られそうですが…。

ともあれ、以下がそのコードです。

#作成日時 : 14:52 2015/02/20
#作成者名：T.Nakano
 rm(list=ls())
 gc();gc()

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

# データの読み込み -------------------------------------------------------------
# データが整理されている csv ファイルを選択する。
 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)
 DT <- read.csv(fileName, header=T)

 # DT に全長のデータが入っている。
 # View(DT)

 # 混合正規分布の分解 -----------------------------------------------------------
 # ヒストグラムを書く -----
 pitch=(max(DT$TL)-min(DT$TL))/15     # だいたいバーが45本
 BKS <- seq(-pitch, max(DT$TL)+pitch*2, 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(コホート番号=1:length(int$x), 平均体長=round(NM$mu), 標準偏差=round(NM$sigma, 2), 個体数=round(NM$lambda*length(DT$TL)))

​

テストデータをこちらに置いておきます。

今回のコードの実行用に使用するなり、自分のデータを整形する参考にするなり、ご自由にお使い下さい。

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