R_Divisive_Clustering_diana関数とnj関数 の変更点 - アールメカブ

アールメカブ


R_Divisive_Clustering_diana関数とnj関数 の変更点


パプアニューギニアとオセアニアの諸言語の系統発生を,それらの分布的特徴を使って分類するという課題である.

Baayen [[Analyzing Linguistic Data:http://www.amazon.co.jp/Analyzing-Linguistic-Data-Introduction-Statistics/dp/0521882591/]] のpp.154--

 phylogeny.dist <- dist(phylogeny[, 3:ncol(phylogeny)], 
          method = "binary")

# 一方の言語のラベルを大文字に変更し,他方は小文字のままにして対照する
 plotname <- as.character(phylogeny$Language)
 plotname[phylogeny$Family == "Papuan"] <-  
     toupper(plotname[phylogeny$Family == "Papuan"])

# plot する.ここで引数 which.plot = 2 の意味は,print.diana がデフォルトで出力するふたつのグラフの二つめのみを表示するよう指定している.

- &size(18){&color(red){Divisive Clustering};};

 plot(diana(phylogeny.dist), labels = plotname, 
 cex = 0.8,  main = " ", xlab = " ",  which.plot = 2)

なお Baayen では  col = c("black","white") と言う引数が与えられているが,これは print.diana メソッドに固有の定義はないようである.


- &size(18){&color(red){neighbor-joining による unrooted tree  を作成する};};
 library(ape)
 
 phylogeny.dist.tr <- nj(phylogeny.dist)
 
# nj  オブジェクトのラベルはデフォルトでは数値なので変更する 
 families <- 
     as.character(phylogeny$Family[as.numeric
         (phylogeny.dist.tr$tip.label)])
 
 languages <- 
     as.character(phylogeny$Language[as.numeric
         (phylogeny.dist.tr$tip.label)])
 
 phylogeny.dist.tr$tip.label <- languages
 
 plot(phylogeny.dist.tr, type = "u", 
       font = as.numeric(as.factor(families)))

引数 type = "u"  は unrooted tree を指定.また各ノードは少なくとも三つの枝を持つ.
パプア諸語とオセアン諸語がきれいに分類されているのが分かる.

#ref(baayen156-2.png,left,nowrap,Divisive_Clustering))

二つまとめて png 画像にする.十分なマージンが必要なので 
par(mfrow = c(1,2), oma = c(1,1,1,1), mar = c(1,1,1,1), bg = "white"); dev.print(device = png, file = "baayen156.png",width = 480, height = 480)
と実行する.

- &color(red){CONSENSUS TREE}; の作り方
- &size(18){&color(red){CONSENSUS TREE}; の作り方};

これはブートストラップによるTreeに観察されないサブグループを collapse させる手法である.その結果としてmultichotomy であるTreeが得られる.今の Baayen の例を使って実行してみる.

# p.158 # Bootstrap によるクラスタリング結果の評価

 B <- 200
 btr <- list()
 length(btr) <- B
 
 for(i in 1:B){
   trB <- nj(dist(papuan.mat[, sample(ncol(papuan.mat),
          replace = T)], method = "binary"))
  trB$tip.label <- 
 as.character(papuan.meta$Language[as.numeric
            (trB$tip.label)])
  btr[[i]] <- trB
 }
 
 props <- prop.clades(papuan.dist.tr, btr) / B
 plot(papuan.dist.tr, type = "u", 
     font = as.numeric(as.factor(fonts)))
 nodelabels(thermo = props, plecol = c("black",
  , "grey"))
 btr.consensus <- consensus(btr, p = 0.5) # 時間がかかる
 x <- btr.consensus$tip.label
 x <- data.frame(Language = x, Node = 1:length(x))
 x <- merge(x, papuan.meta, 
        by.x = "Language", by.y = "Language")
 head(x)
 x <- x[order(x$Node),]
 x$Geography <- as.factor(x$Geography)
 plot(btr.consensus, type = "u", font = 
  as.numeric(x$Geography))

#ref(baayen159-2.png,center,nowrap,Consensus Clustering)