6_phyloglm.rmd

---
title: "Niche Traits"
author: "Rutger Vos (@rvosa)"
date: "27-2-2019"
output: html_document
---

Load the required libraries:

```{r load_libs}
library(dismo, quietly=T)
library(phytools, quietly = T)
library(phylolm, quietly = T)
library(dplyr, quietly = T)
source(Data.R)
```

Initialize global variables:

```{r globals}
REPO_HOME <- paste(getwd(),'/../',sep='')
```

Load the list of taxa:

```{r taxa}
taxa.names <- scan(paste(REPO_HOME, "/data/filtered/taxa.txt", sep = ""), sep = "\n", what = character())
```

Define the list of GIS layers, i.e. the traits:

```{r layers}
gis.layers <- raster::getData(
  "worldclim",
  var = "bio",
  res = 5,
  path = paste(REPO_HOME, "/data/GIS", sep = ""),
  download = T
)

# The location of the TIFF file with the stacked layers described directly above
files.names <- list.files(paste(REPO_HOME, "/data/GIS/5_deg", sep = ""))

# Turn the file names into layer names: strip the prefix (which might include
# the resolution) and strip the file extension
gis.layers.names <- files.names
gis.layers.names <- gsub('current_5arcmin_','',gis.layers.names)
gis.layers.names <- gsub('.tif','',gis.layers.names)

# Combine the layer names with those we've already read from BIOCLIM
gis.layers.names <- c(names(gis.layers),gis.layers.names)

# Iterate over files
for (i in 1:length(files.names)) {
  
  # Stack with previously read layers
  gis.layers <- stack(
    gis.layers,
    
    # Read as raster
    raster(
      
      # Construct file name
      paste(REPO_HOME, "/data/GIS/5_deg/", files.names[i], sep = "")
    )
  )
}

# Apply all names
names(gis.layers) <- gis.layers.names
```

Phylogenetic generalized linear modeling
----------------------------------------

We are going to look at the trait means as predictor variables for domestication, so 
we create a data frame with domestication as a state:

```{r domestication}
traits.means.file <- sprintf('%s/results/OMI/normalized_MaxEnt_values.csv', REPO_HOME)

# do the averaging if we haven't already
if (!file.exists(traits.means.file) ) {

  # instantiate data frame from a matrix of taxa x layers
  traits.means <- data.frame(matrix(
    nrow = length(taxa.names), 
    ncol = length(layers.names) + 1, 
    dimnames = list(taxa.names,c("IsDomesticated",layers.names)))
  )
  
  # vector of all domesticated taxa
  dom.taxa <- c(
    "Bos_frontalis_gaurus", "Bos_grunniens_mutus", "Bos_javanicus", "Bos_taurus_primigenius",
   	"Bubalus_bubalis_arnee", "Camelus_bactrianus", "Camelus_dromedarius", "Capra_hircus_aegagrus",
   	"Equus_przewalskii", "Lama_glama_guanicoe", "Ovis_aries_orientalis",
   	"Rangifer_tarandus", "Sus_scrofa", "Vicugna_vicugna"
  )
  
  # iterate over taxa
  for ( i in 1:length(taxa.names) ) {
    taxon.name <- taxa.names[i]
    message(taxon.name)
    
    # store whether this taxon is domesticated
    if ( taxon.name %in% dom.taxa ) {
      traits.means[i,1] <- T
    } else {
      traits.means[i,1] <- F
    }
    
    # read occurrences, create df with numeric coordinates
    taxon.occ <- get_occurrences(REPO_HOME,taxon.name)
    taxon.coord <- as.data.frame(cbind(taxon.occ[, c("decimal_longitude", "decimal_latitude")]))
    taxon.coord$decimal_longitude <- as.numeric(as.character(taxon.coord$decimal_longitude))
    taxon.coord$decimal_latitude <- as.numeric(as.character(taxon.coord$decimal_latitude))  
    
    # extract data, store in mean values
    taxon.values <- extract(data.layers, taxon.coord, df = T, fun = 'mean')
    for ( j in 1:length(layers.names) ) {
      layer.name <- layers.names[j] 
      traits.means[taxon.name,layer.name] <- mean(taxon.values[,layer.name], na.rm = T)
    }  
  }
  
  # write result
  write.csv(
    traits.means, 
    file = traits.means.file, 
    quote = F, 
    row.names = T
  )
} else {
  
  # read means from file
  traits.means <- read.csv(
    traits.means.file, 
    header = T, 
    sep = ",", 
    row.names = 1
  )
}
```

And load a phylogenetic tree:

```{r phylo}
# reconcile tree and traits
tree.file <- sprintf('%s/data/phylogeny/ungulates.tree',REPO_HOME)
tree <- read.tree(file = tree.file)
tip.idx <- match(row.names(traits.means),tree$tip.label)
tip.idx <- tip.idx[!is.na(tip.idx)]
tree <- keep.tip(tree,tip.idx)
traits.means.tree <- traits.means[tree$tip.label,]

# run variable selection
pglm.formula.tmpl <- 'IsDomesticated  ~ bio4 + %s'
aic <- vector(mode = "numeric", length = length(gis.layers.names))
lnl <- vector(mode = "numeric", length = length(gis.layers.names))
results <- data.frame( aic = aic, lnl = lnl, row.names = gis.layers.names)

for ( i in 1:length(gis.layers.names) ) {
  l <- gis.layers.names[i]
  f <- sprintf(pglm.formula.tmpl, l)
  message(l)
  r <- phyloglm(f, traits.means.tree, tree, method = "logistic_MPLE", btol = 30)
  results[l,'lnl'] <- r$logLik
  results[l,'aic'] <- r$aic
}

```

After iterating over all layers, we find the layer that produces the lowest AIC
in the results df. We add that layer to the formula, and then iterate again to
find the layer that minimizes the AIC in combination. If we keep doing this we
find the combination that overall minimizes the AIC:

MaxEnt 
bio4 67.32846

raw occurences 

bio3 69.77260

bio3 + OrganicCarbon_5min 68.11076
bio3 + OrganicCarbon_5min + bio2  67.96198
bio3 + OrganicCarbon_5min + bio2 + bio17 63.80596
bio3 + OrganicCarbon_5min + bio2 + bio17 + BulkDensity_5min 61.71878




```{r plot3d}
f <- IsDomesticated  ~ bio3 + OrganicCarbon_5min + bio2 + bio17 + BulkDensity_5min
r <- phyloglm(f, traits.means.tree, tree, method = "logistic_MPLE", btol = 30)
summary(r)


p3d <- as.matrix(dplyr::select(traits.means.tree, bio3, OrganicCarbon_5min, bio2, bio17, BulkDensity_5min))
phylomorphospace3d(tree,p3d)
```

Let's plot some of these values on the niche clustering dendrogram:

```{r trace}
# load tree
# CSV file location of all pairwise comparisons of niche overlap in Gowers's D
overlap.file <- paste(REPO_HOME, '/results/maxent/phylo_tree_gower.tree', sep = '')
niche.tree <- read.tree(overlap.file)


dom.tip.idx <- match(dom.taxa,niche.tree$tip.label)
dom.tip.idx <- dom.tip.idx[!is.na(dom.tip.idx)]

# iterate over layers
for ( i in 1:length(layers.names) ) {
  layer.name <- layers.names[i]
  layer.file <- sprintf('%s/probeersel/%s_contMap.pdf',REPO_HOME,layer.name)
  
  # write continuous mapping if not exists
  if (!file.exists(layer.file)) {
  
    # extract named vector of local layer values, subset to match tips in trees
    trait.vector <- traits.means[,layer.name]
    names(trait.vector) <- row.names(traits.means)
    trait.vector <- trait.vector[niche.tree$tip.label]
    
    # do the mapping
    cm <- contMap(
      niche.tree,
      trait.vector,
      plot = F
    )
    
    
    # write result
    pdf(file = layer.file)
    plot(
      cm,
      type = "fan",
      fsize = 0.2,
      lwd = 2,
      outline = F,
      ftype = 'i'
    )  
    tiplabels(
      tip = dom.tip.idx, 
      frame = "none", 
      pch = 20
    )
    dev.off()
  }

}
```


Filtering on trait values
-------------------------

Let's see which other ungulates fall inside the value range of the domesticated ones for
all traits. We are both going to consider all the domesticates as well as the set that
leaves out the reindeer (which are often considered semi-domesticated).

```{r filter}
# split data in domesticated and non, extract trait names vector
traits.means$taxon.name <- row.names(traits.means)
dom.df <- filter(traits.means,IsDomesticated)
undom.df <- filter(traits.means,!IsDomesticated)
layers.names <- names(dplyr::select(dom.df,-IsDomesticated,-taxon.name))

# to exclude the reindeer from the value ranges:
#dom.df <- filter(dom.df, taxon.name != 'Rangifer_tarandus')

# iterate over layers
for ( i in 1:length(gis.layers.names) ) {
  layer <- gis.layers.names[i]
  dom.values <- dom.df[,layer]
  min <- min(dom.values)
  max <- max(dom.values)
  undom.df <- filter(undom.df, UQ(as.name(layer)) >= min & UQ(as.name(layer)) <= max)
}

# clean up result
row.names(undom.df) <- undom.df$taxon.name
undom.df <- dplyr::select(undom.df, -IsDomesticated, -taxon.name)
row.names(undom.df)
```