Intro

We investigate how phylogeny drives the clustering of mammals species in the Wallacea region.

library(ape)
library(phytools)
## Loading required package: maps
library(methClust)
library(rasterVis)
## Loading required package: raster
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following objects are masked from 'package:ape':
## 
##     rotate, zoom
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
library(gtools)
library(sp)
library(rgdal)
## rgdal: version: 1.2-20, (SVN revision 725)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-7
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
## 
##     layer
library(maps)
library(mapdata)
library(mapplots)
library(scales)
library(ggthemes)

Phylo function

phylo.counts = function(counts, tree, collapse_at){
  root_node <- length(tree$tip.label) + 1
  root_age <- ape::branching.times(tree)[names(ape::branching.times(tree)) == root_node]
  trees_at_slice <- phytools::treeSlice(tree, root_age - collapse_at)
  counts_at_slice <- base::as.data.frame(counts)
  for( i in 1:length(trees_at_slice)){
   # print(i)
    new.column <- base::as.data.frame(base::rowSums(counts[,trees_at_slice[[i]]$tip.label]))
    colnames(new.column) <- trees_at_slice[[i]]$tip.label[1]
    drops <- trees_at_slice[[i]]$tip.label
    counts_at_slice <- counts_at_slice[,!(names(counts_at_slice) %in% drops)]
    counts_at_slice <- cbind(counts_at_slice, new.column)
  }
  counts_at_slice <- base::as.matrix(counts_at_slice)
  return(counts_at_slice)
}

Data

birds <- get(load("../data/wallacea_birds.rda"))
latlong_chars <- rownames(birds)
latlong <- cbind.data.frame(as.numeric(sapply(latlong_chars, 
                              function(x) return(strsplit(x, "_")[[1]][1]))),
         as.numeric(sapply(latlong_chars, 
                           function(x) return(strsplit(x, "_")[[1]][2]))))
phylo_names <- read.csv("../data/names_matched_to_phylogeny.csv", row.names = 1)
species_to_consider <- intersect(phylo_names[,1], colnames(birds))
birds2 <- birds[,species_to_consider]
colnames(birds2) <- phylo_names[match(species_to_consider, phylo_names[,1]),2]

Birds

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
seq2 <- seq(5, 50, length.out = 10)

Model

phylo_gom <- list()
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 <- 10
  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=2, 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]]
}
save(phylo_gom, file = "../output/phylo_gom_birds/gom_2.rda")

Phylogenet cluster plots

phylo_gom <- get(load("../output/phylo_gom_birds/gom_2.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_2/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()
}

This R Markdown site was created with workflowr