phylogenetic cluster birds : K = 10

Create phylogenic counts data

tree_file <-  read.tree("../data/birds.nwk")
colnames(birds2) <- gsub(" ","_",colnames(birds2))
tips_to_drop <- setdiff(tree_file$tip.label, colnames(birds2))
new_tree <- drop.tip(tree_file, tips_to_drop, 
                     trim.internal = TRUE, subtree = FALSE, root.edge = 0,
                     rooted = is.rooted(tree_file),
                     collapse.singles = TRUE, interactive = FALSE)
if(is.null(new_tree)){
  new_tree <- tree_file
}
birds_data <- birds2

Model

phylo_gom <- list()
seq2 <- seq(5, 50, length.out = 10)
for(num in 1:length(seq2)){
  birds_phylo <- phylo.counts(birds_data, new_tree, collapse_at = seq2[num])
  cat("Processing over", dim(birds_phylo)[2], "species \n")
  birds_phylo[birds_phylo > 0] <- 1
  M <- 1
  fits_list <- list()
  L_array <- c()
  for(m in 1:M){
    counter = 0
    while(counter != 1){
      tmp <- try(meth_topics(birds_phylo, 1 - birds_phylo,
                             K=10, tol = 10,
                             use_squarem = FALSE), TRUE)
      if(!inherits(tmp, "try-error")){
        counter = 1
      }else{
        counter = 0
      }
    }
    fits_list[[m]] <- tmp
    cat("We are at iteration", m, "\n")
  }
  
  loglik <- unlist(lapply(fits_list, function(x) return(x$L)))
  ids <- which.max(loglik)
  phylo_gom[[num]] <- fits_list[[ids]]
}

## Processing over 1531 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 177050.8, 1749.2, 2152.9, 1227, 1009.1, 550.2, 895.9, 319.3, 257, 129.3, 181.4, 39.4, 72.7, 28.9, 45.3, 39, 42.5, 50.5, 38.5, 33, 33.2, 58.5, 55.4, 75.6, 50.9, 36, 37.8, 15.1, 35, 21.7, 26.6, 74.6, 27.7, done.
## We are at iteration 1 
## Processing over 936 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 57666.5, 1582.2, 740.9, 480.7, 262.8, 183.2, 175.6, 121.7, 63.3, 52.9, 80.7, 31.4, 32.5, 20.3, 27.9, 32.2, 56.5, 23.5, 31.5, 21.6, 14.1, 10.4, 43.2, 42, 72.4, 5.6, done.
## We are at iteration 1 
## Processing over 590 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 44574.9, 1277.1, 344.3, 570, 250.5, 170.5, 165.3, 94.3, 71.5, 81.4, 37.5, 37.8, 44.4, 39.3, 21, done.
## We are at iteration 1 
## Processing over 390 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 53671.6, 900.9, 346.3, 506.7, 1205.6, 186, 138.9, 108.3, 32.3, 51.5, 44.7, 10.7, done.
## We are at iteration 1 
## Processing over 271 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 23785.5, 678.6, 787.1, 510.7, 163.5, 94.3, 100, 71.6, 41.8, 17.5, 14.4, done.
## We are at iteration 1 
## Processing over 185 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 8683, 147.3, 188.2, 108.3, 145.3, 59.4, 124.9, 112.1, 62.1, 52.9, 17.1, 32.3, done.
## We are at iteration 1 
## Processing over 122 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 9890.4, 363.8, 141.7, 104.4, 153.1, 39.3, 99.1, 41.9, 31.2, done.
## We are at iteration 1 
## Processing over 90 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 4814.1, 70.8, 10.9, done.
## We are at iteration 1 
## Processing over 73 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 6719.6, 69.7, 151.1, 20.4, done.
## We are at iteration 1 
## Processing over 59 species 
## 
## Estimating on a 664 samples collection.
## log posterior increase: 7844.4, 95, done.
## We are at iteration 1

save(phylo_gom, file = "../output/phylo_gom_birds/gom_10.rda")

Phylogenet cluster plots

phylo_gom <- get(load("../output/phylo_gom_birds/gom_10.rda"))

color = c("red", "cornflowerblue", "cyan", "brown4", "burlywood", "darkgoldenrod1",
          "azure4", "green","deepskyblue","yellow", "azure1")
intensity <- 0.8
for(num in 1:length(seq2)){
  png(filename=paste0("../docs/phylo_birds_K_10/geostructure_phylo_", seq2[num], ".png"),width = 1000, height = 800)
map("worldHires",
    ylim=c(-18,20), xlim=c(90,160), # Re-defines the latitude and longitude range
    col = "gray", fill=TRUE, mar=c(0.1,0.1,0.1,0.1))
lapply(1:dim(phylo_gom[[num]]$omega)[1], function(r)
  add.pie(z=as.integer(100*phylo_gom[[num]]$omega[r,]),
          x=latlong[r,1], y=latlong[r,2], labels=c("","",""),
          radius = 0.5,
          col=c(alpha(color[1],intensity),alpha(color[2],intensity),
                alpha(color[3], intensity), alpha(color[4], intensity),
                alpha(color[5], intensity), alpha(color[6], intensity),
                alpha(color[7], intensity), alpha(color[8], intensity),
                alpha(color[9], intensity), alpha(color[10], intensity),
                alpha(color[11], intensity))));
dev.off()
}