アールメカブ

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

Baayen Analyzing Linguistic Data の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 がデフォルトで出力するふたつのグラフの二つめのみを表示するよう指定している．

• Divisive Clustering

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

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

• 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 を指定．また各ノードは少なくとも三つの枝を持つ． パプア諸語とオセアン諸語がきれいに分類されているのが分かる．

二つまとめて 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) と実行する．

• 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))