HPsubfield_eco.Rmd

---
title: "Ecological Correlates of Hippocampal Subfield Volumes Across Primates" 
author: |
  | Brian M. Schilder, Bioinformatician II
  | Raj Lab
  | Department of Neuroscience
  | Icahn School of Medicine at Mount Sinai
  | NYC, New York
date: "4/29/2017"
output:
  rmarkdown::html_document: 
    theme: yeti
    highlight: zenburn
    code_folding: show 
    toc: true 
    toc_float: true
    smooth_scroll: true
    number_sections: false  
editor_options: 
  chunk_output_type: console
---

# Setup

Read comments in the setup section below for directions on install additional dependencies.

```{r setup, warning=F, message=F, include=T, dpi = 600}
setwd("~/Desktop/Hippo_Eco/")
library(readxl)
library(dplyr)  
library(pander) 
library(reshape2)
library(shiny)
library(DT)

library(ggplot2) #devtools::install_github('cran/ggplot2')    
library(gplots)
library(corrplot) 
library(ggdendro) 
library(ggrepel) 
library(ggbiplot) #devtools::install_github("vqv/ggbiplot")
library(gridExtra)
library(grid) 
library(cowplot)
library(plotly)
library(heatmaply)
library(tidyverse)

# library(xlss) # For xlss, must first download JDK: https://www.oracle.com/technetwork/java/javase/downloads/jdk11-downloads-5066655.html
library(ggtree)# BiocManager::install("ggtree")
library(ggimage)
library(ggstance)


library(phytools) 
library(caper)
library(ape) 
library(dispRity)
# install.packages("https://cran.r-project.org/src/contrib/Archive/genlasso/genlasso_1.3.tar.gz", repos=NULL, method="libcurl")
# devtools::install_github("khabbazian/l1ou")
library(l1ou)

# For interactive 3D plotting
library(knitr)
library(rgl) # Macs need to download X11 first. Link here: http://xquartz.macosforge.org/trac/wiki
knitr::knit_hooks$set(webgl = hook_webgl)
knitr::opts_chunk$set(echo=T, warning=F, message=F)

# thm <- knitr::knit_theme$get("solarized-dark")
# knitr::knit_theme$set(thm)

sessionInfo()
```


# Import Data

```{r Import Data}
#### LOAD PREVIOUSLY COMPUTED RESULTS #### #### #### #### #### #### 
# load("~/Desktop/Hippo_Eco/Data/HippoEco_AllResults.RData")
#### #### #### #### #### #### #### #### #### #### #### #### #### 

# Frahm HP-subfield data (1982) + Barger whole-HP & brain (2015) + Pantheria eco + ATWP eco
HP_orig <- as.data.frame(read_excel("Data/Frahm_HPsub_eco.xlsx", sheet = "Final_Data", na = "NA"))

## Designate groups
HP_orig$Genus <- gsub("_.*","",HP_orig$Species_10K)

HP_orig$Clade <- 
  ifelse(grepl("Homo", HP_orig$Genus), "Homo sapiens", 
         ifelse(grepl("Gorilla|Hylobates|Pan", HP_orig$Genus), "Apes", 
         ifelse(grepl("Cercopithecus|Erythrocebus|Lophocebus|Miopithecus|Nasalis|Papio|Piliocolobus|Pygathrix", HP_orig$Genus), "Old World Monkeys", 
         ifelse(grepl("Alouatta|Aotus|Ateles|Callicebus|Callimico|Callithrix|Cebus|Lagothrix|Pithecia|Saguinus|Saimiri", HP_orig$Genus), "New World Monkeys",
         ifelse(grepl("Tarsius", HP_orig$Genus), "Tarsiiformes",
         ifelse(grepl("Avahi|Cheirogaleus|Eulemur|Indri|Lepilemur|Microcebus|Propithecus|Varecia|Galago|Galagoides|Loris|Nycticebus|Otolemur|Perodicticus|Daubentonia", HP_orig$Genus), "Strepsirrhines",
         "NA"))))))
# Order Clade for plots
orderedClades <- c("Homo sapiens","Apes","Old World Monkeys","New World Monkeys","Tarsiiformes","Strepsirrhines")
HP_orig$Clade <- factor(HP_orig$Clade, labels=orderedClades, levels=orderedClades, ordered=T)


# Rename some variables
HP_orig <- HP_orig %>% dplyr::rename(retroHP=HippocampusRetro,
                                      FD=FasciaDentata,
                                      Sub=Subiculum,
                                      Fibers=HP_HS_fibers,
                                      Body_Size=Body,
                                      # "EC*"=Schizocortex,
                                      EC=Schizocortex, 
                                      Diet_Breadth = DietBreadth,
                                      Group_Size = GroupSize_filled,
                                      Home_Range = HR_filled,)
 

# Log Neuroanatomical vars
#cols <- c("Body_Size","BrainVol","retroHP","FD","Hilus","CA3","CA2","CA1","Sub")
#HP_orig[cols] <- log(HP_orig[cols])

# Log Eco vars
# HP_orig$PopulationDensity_n_km2 <- scales::rescale(log(HP_orig$PopulationDensity_n_km2),c(1,10))
# HP_orig$GroupSize_filled <-  scales::rescale(log(HP_orig$GroupSize_filled),c(1,10))
# HP_orig$Home_Range <-  scales::rescale(log(HP_orig$Home_Range),c(1,10))


# Add combined regions as variables
HP_orig$retroHP.EC <- HP_orig$retroHP + HP_orig$EC
HP_orig$DG <- HP_orig$FD + HP_orig$Hilus
HP_orig$CA2.3 <- HP_orig$CA2 + HP_orig$CA3

# Correct units of Barger (2015) data
HP_orig[,c("BrainVol","HPvol")] <- HP_orig[,c("BrainVol","HPvol")] * 1000


# Create retoHP minus each subregion
  HP_orig$Body__Brain <- HP_orig$Body_Size - HP_orig$BrainMass
  HP_orig$Brain__HP <- HP_orig$BrainVol - HP_orig$retroHP
  HP_orig$Brain__EC <- HP_orig$BrainVol - HP_orig$EC
  HP_orig$Brain__Fib <- HP_orig$BrainVol - HP_orig$Fibers

  HP_orig$retroHP__DG <- HP_orig$retroHP - HP_orig$DG
  HP_orig$retroHP__FD <- HP_orig$retroHP - HP_orig$FD
  HP_orig$retroHP__Hil <- HP_orig$retroHP - HP_orig$Hilus
  HP_orig$retroHP__CA2.3 <- HP_orig$retroHP - HP_orig$CA2.3
  HP_orig$retroHP__CA3 <- HP_orig$retroHP - HP_orig$CA3
  HP_orig$retroHP__CA2 <- HP_orig$retroHP - HP_orig$CA2
  HP_orig$retroHP__CA1 <- HP_orig$retroHP - HP_orig$CA1
  HP_orig$retroHP__Sub <- HP_orig$retroHP - HP_orig$Sub
colnames(HP_orig)


# Convert ACtivityCycle to categorical
# ActCycle_dict <- setNames(c("A","B","C",NA), unique(HP_orig$ActivityCycle))
# HP_orig$ActivityCycle <- ActCycle_dict[HP_orig$ActivityCycle]
## Make shapes and colors consistent
shapes <- c(17,16,18,15,7,9)
names(shapes) <-levels(HP_orig$Clade) 
#gg_color_hue <- function(n) { # ggplot default colors
  #hues = seq(15, 375, length=n+1)
  #hcl(h=hues, l=65, c=100)[1:n]
  #}
#colors <- gg_color_hue(6)
colors <- c("#C99800","#F8766D","#00BCD8","#00BA38","blue2","#FD61D1")
# Remove human shape/color if humans were removed
if(!("Homo" %in% HP_orig$Genus)){colors<-colors[-1]; shapes<-shapes[-1]}
names(colors) <- names(shapes)
```

## Filter & Normalize

```{r Filter & Normalize}
# log normalize neuro and/or eco data
## Don't need to normalize the neuro vars bc they're already normal and on the same scale.
## Besides, I'm using the residuals for the selective regime analyses.
var_list <- c(
              # "Diet_Breadth", # Can't normalize bc it's ordinal data, not numeric like the others
              "Population_Density", # "PopulationDensity_n_km2",
              "Group_Size",
              "Home_Range")

# vec <- HP$data$CA1_res

scale_ <- function(vec, lower=1, upper=2){
  scaled <- as.vector(scales::rescale(vec, to=c(lower:upper))) 
  return(scaled)
}

log_scale <- function(vec){
   logged <- log10(vec)
   scaled <- scale_(logged)
   # scaled <- scaled + 0.00000000000000001
   # scaled <- zscore(logged, rescale_01 = T)
   return(scaled)
}

scale_log <- function(vec){
   scaled <- scale_(vec)# + .Machine$double.xmin
   logged <- log10(scaled)
   return(logged)
}


scale_minmax <- function(vec){ 
  scaled <- vec + abs(min(vec)) + 0.1
  return(scaled)
}

logTrans <- function(vec, constant=1){
  logged <- log10(vec + constant - min(vec))
  return(logged)
}

logTrans_scale <- function(vec, constant=1){
   logged <- log10(vec + constant - min(vec))
   scaled <- scale_(logged)
   return(scaled)
}

scale_logTrans <- function(vec, constant=1){
  scaled <- scale_(vec)
  logged <- logTrans(scaled)
  return(logged)
}

scaleMM_logTrans <- function(vec, constant=1){
  scaled <- scale_minmax(vec)
  logged <- logTrans(scaled, constant = constant)
  return(logged)
}

logTrans_missing <- function(vec){ 
  indices <- match(NA, vec)
  if(is.na(indices)){
     completed <- scale_logTrans(vec)
  } else {
     non.missing <- vec[!is.na(vec)] %>% scale_logTrans()
     completed <- R.utils::insert(x=non.missing, ats=indices, values=rep(NA, length(indices)))
  } 
  return(completed)
}

logTrans_missing_scale <- function(vec){
  trans <- logTrans_missing(vec)
  scaled <- scale_(trans)
  return(scaled)
}

rm_outliers <- function(vec, n_outliers=1){
  # vec <- HP$data$Population_Density
  outliers <- boxplot(vec, plot=F)$out %>% sort(decreasing = T)
  if(length(outliers)>0){
    indices <- match(outliers[n_outliers], vec)
    vec[indices] <- NA
  } 
  return(vec)
}



log_missing <- function(vec){ 
  indices <- match(NA, vec)
  if(is.na(indices)){
     completed <- scale_logTrans(vec)
  } else {
     non.missing <- vec[!is.na(vec)] %>% log10()
     completed <- R.utils::insert(x=non.missing, ats=indices, values=rep(NA, length(indices)))
  } 
  return(completed)
}

log_missing_scale <- function(vec){
  logged <- log_missing(vec)
  scaled <- scale_(logged)
  return(scaled) 
}

# HP_scaled <- HP_orig
# for(v in var_list){
#   print(v)
#   scaled <- log_scale(HP_orig[[v]] )
#   res <- shapiro.test(scaled)
#   print(res) 
#   HP_scaled[v] <- scaled
#   # hist(HP_orig[[v]])
#   # hist(scaled)
#   # From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.
# }
# HP_orig <- HP_scaled


# Apply filter to find variables with more than ##% of species missing that data, so that they can be excluded from the model (or have data added to them)
# HP_orig <- HP_orig %>% dplyr::select(which(colMeans(is.na(.)) < 0.40))
# ************** Exclude strepsirrhines # ************** 
# HP_orig <- subset(HP_orig, Clade != "Strepsirrhines")
```

## Tree & Comparative Files

```{r Tree & Comparative Files}
# Primate tree from 10K Trees 
tree10K <- read.nexus(file = "Data/Phylo/10k_Primates.nex")

## Rescale tree
# Need to rescale the tree due to errors with pgls described here: http://blog.phytools.org/2011/12/error-message-from-brownielite-and.html#
#scale <- 100
#tree10K$edge.length <- tree10K$edge.length/max(nodeHeights(tree10K)[,2])*scale

# Create comparative.data
HP <- comparative.data(phy=tree10K, data=HP_orig, names.col=Species_10K, 
                        vcv=T, na.omit=F, warn.dropped=T)
# Put Species_10K back as col too
HP$data["Species_10K"] <- row.names(HP$data)
```

# Create Functions


## Assign norm function
```{r}
# res_norm_func=function(x){x}
# resHR_norm_func=scale_
# norm_func=scale_
```


## PGLS residuals & plot function

```{r PGLS residuals & plot function}

pgls.plot_curve <- function(comp_data, x, y, model, xy_plot, x_norm_func, y_norm_func){
  x.values <- x_norm_func(as.numeric(comp_data$data[[x]]))
  x.smooth<- seq(min(x.values),max(x.values),.01)
  
  exp.predict <- exp(predict(model, setNames(list(x.smooth), x)))
  exp.line <- data.frame(x=x.smooth, y=exp.predict)
  return( xy_plot + geom_line(data = exp.line, aes(x, y), inherit.aes = F) )
}


pgls.xy_plot <- function(comp_data, pgls_summary, 
                           response_var, predictor_var, 
                           norm_func=function(x){x},
                           x_prefix="", x_suffix="", 
                           y_prefix="", y_suffix="", interact=F){
    xy_plot_legend <- ggplot(comp_data$data, aes(y=norm_func(eval(parse(text=response_var))), 
                                            x=norm_func(eval(parse(text=predictor_var))),
                                            fill=Clade, color=Clade, shape=Clade) )  +
    geom_abline(slope = pgls_summary$coefficients[2,1], intercept = pgls_summary$coefficients[1,1], color="black") +
    geom_point(size=3) + 
    labs(y=paste0(y_prefix,response_var,y_suffix), x=paste0(x_prefix, predictor_var,x_suffix)) + 
    theme(plot.title = element_text(hjust=0.5, size=8)) + 
    scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + scale_fill_manual(values=colors)   
  if(interact){
    print(xy_plot_legend %>% ggplotly())
  } else{
    print(xy_plot_legend)
  }
  xy_plot <- xy_plot_legend + theme(legend.position = "none")
  return(xy_plot)
}

pgls_function <- function(comp_data, 
                          response_var, 
                          predictor_var,  
                          lambda='ML', 
                          interact=T, 
                          norm_func=log10,#log_missing,    
                          y_prefix="",
                          x_prefix="", 
                          y_suffix=" (mm3)",
                          x_suffix=" (mm3)"){   
  model <- caper::pgls(data=comp_data, 
                        norm_func(eval(parse(text=response_var))) ~ 
                        norm_func(eval(parse(text=predictor_var))), lambda=lambda) 
  pgls_summary <- summary(model)
  print(pgls_summary)
  ## Get res
  res <- resid(model)
  res_name <- paste(response_var,"res",sep="_")
  ## Put back into comparative.data
  comp_data$data[res_name] <- as.vector(res)# as.vector(scales::rescale(res, c(1:2))) #
  new_comp_data <- comp_data
  # Put into actual comparative.data in global environment
  # comp_data$data[res_name] <- as.vector(scales::rescale(res, c(1:2)))
  assign('HP',comp_data,envir=.GlobalEnv)
  
  #  Test if residuals are normally distributed
  #shapiro.test <- shapiro.test(x = comp_data$data[paste(response_var,"res",sep="_")])
  
  # Residuals XY Plot
    # radj <- as.numeric(pgls_summary$adj.r.squared)
    # p <- as.numeric(pgls_summary$coefficients[2,4])
    #title = title=paste(response_var,"vs.", predictor_var, paste("\n adj.R^2=", round(radj,4),", p= ",round(p,4)))
  
   xy_plot <- pgls.xy_plot(comp_data, pgls_summary, 
                           response_var, predictor_var )
    
   # xy_plot <- pgls.plot_curve(comp_data, 
   #                            x=predictor_var, y=response_var, 
   #                            model=model, xy_plot=xy_plot, 
   #                            x_norm_func = norm_func,
   #                            y_norm_func = function(x){x})
    
 
  # + annotate("text",x=5, y=10, hjust=0, label =  paste("adj.R^2 =", round(radj,4), "; p = ", round(p,4)))
  

# Plot fancy contMaps

## Absolute vols plot (not that informative since it's basically just tracking brain size)
  #HP_data <- as.matrix(comp_data$data[response_var])
  #names(HP_data) <- comp_data$phy$tip.label
  # phylo_plot <- contMap(comp_data$phy, HP_data, type="phylogram", plot=TRUE)
## Residuals plot
  HPres_data <- as.matrix(comp_data$data[res_name])
  names(HPres_data) <- comp_data$phy$tip.label
  res_phylo_plot <- contMap(comp_data$phy, HPres_data, type="phylogram", 
                            plot=FALSE, fsize=c(0.5, 1))
    plot(res_phylo_plot, leg.txt=paste(res_name),lwd=3, fsize=c(0.5, 1))

# Use plot() to edit further
  #plot(res_phylo_plot)
  #par(cex.lab=.5)

# Top residuals
  #res <- DG.HP_res$residuals
  res_data <- data.frame(species=row.names(res),res=res)
 
  
  print(paste("Top 5 Positive Residuals:", res_name))
  sort_res <- res_data[with(res_data, order(-res)), ]
  row.names(sort_res) <- NULL
  print(sort_res[1:5,])
  print("=====================================================")
  print(paste("Top 5 Negative Residuals:", res_name))
  sort_res <- res_data[with(res_data, order(res)), ]
  row.names(sort_res) <- NULL
  print(sort_res[1:5,])
# Output  
  output <- list(response_var=response_var, predictor_var=predictor_var,
                 pgls_summary=pgls_summary, residuals=res, sort_res=sort_res,
                 xy_plot= xy_plot, 
                 res_phylo_plot=res_phylo_plot,
                 new_comp_data=new_comp_data)
  return(output)
}
```

## Eco analyses function

**NOTE**: Some variables may not be appropriate to put in model together due to colinearity. To solve this, could:
* Remove certain variables
* Perform a stepwise regression: http://blog.phytools.org/2014/06/performing-stepwise-phylogenetic.html


### Phylostep

Optimize the mode through systematic variable selection.

```{r Phylostep}
run_eco.variable_selection <- function(pgls_function.out,
                                       pgls_eco_raw,
                                       eco_vars,
                                       res_norm_func=function(x){x},
                                       resHR_norm_func=scale_,
                                       norm_func=scale_,
                                       bootstraps=100){
  # pgls_function.out <- CA1.HP_res
  message(">>>>> PHYLOSTEP:: <<<<<<<")
  library(future); plan(multiprocess)
  # Filtering
  dat <- pgls_function.out$new_comp_data
  dat <- subset(dat, Genus!="Homo")
  keep.species <- row.names(dat$data[,c(pgls_function.out$response_var, eco_vars)] %>% na.omit() %>% droplevels())
  dat <- subset(dat, Species %in% keep.species) 
  # Run variable selection
  # text.formula <- paste(pgls_function.out$response_var, "~", paste(eco_vars,collapse="+"))
  text.formula <- paste0(pgls_function.out$response_var,"_res ~ ", as.character(pgls_eco_raw$formula)[3])  
  best.model <- phylolm::phylostep(data = dat$data, 
                                   phy = dat$phy, 
                                   direction = "both", 
                                   model = "lambda",
                                   boot = bootstraps,
                                   formula = eval(parse(text=text.formula)))
  message(">>>>> PHYLOSTEP:: Optimized model results <<<<<<<")
  print(summary(best.model) ) 
  # Re-run with pgls function just because the summary output is nicer 
  message(">>>>> PHYLOSTEP:: New Model results <<<<<<<")
  print(best.model$formula)
  new.model <- caper::pgls(data = dat, formula = eval(parse(text=best.model$formula)))
  summary(new.model)
  
  return(best.model)
} 
```



```{r Eco analyses function} 
nonNA_subset <- function(HP, 
                         response_var,
                         no_homo = T,
                         remove_outliers = F,
                         eco_vars =  c("Diet_Breadth",
                                        "Population_Density",
                                        "Group_Size",
                                        "Home_Range",
                                        "Residual_HomeRange")){
  dat <- HP
  if(no_homo){
     print("Removing Homo sapiens...")
     # Phylo object understand subset command
    dat <- subset(dat, Genus!="Homo")
   } 
  keep.species <- row.names(dat$data[,c(response_var, eco_vars)] %>% 
                              na.omit() %>% droplevels())
  dat <- subset(dat, Species %in% keep.species) 
  
  
 if(remove_outliers!=F){
   print("Removing outliers...")
   for(v in eco_vars){
      x <- as.numeric(dat$data[[v]])
      dat$data[v] <- rm_outliers(x, remove_outliers)
   }
 }  
  return(dat) 
}
# MODEL HERE 
# Analyze all eco variables in one model
run_eco <- function(pgls_function.out,  
                    res_norm_func=function(x){x},
                    resHR_norm_func=scale_,
                    norm_func=scale_,#log_missing_scale,#scale_log, 
                    no_homo=T, 
                    remove_outliers=0,
                    variable_selection=F){
  dat <- pgls_function.out$new_comp_data 
  eco_vars <- c("Diet_Breadth",
                "Population_Density",
                "Group_Size",
                "Home_Range",
                "Residual_HomeRange")
  dat <- nonNA_subset(HP = dat, 
                      response_var = pgls_function.out$response_var,
                       no_homo = no_homo,
                       eco_vars = eco_vars)
 
   pgls_eco_raw <- pgls(data=dat, formula=
    res_norm_func(eval(parse(text=paste(pgls_function.out$response_var,"res",sep="_"))))~  
      
    resHR_norm_func(Diet_Breadth) +  
    norm_func(Population_Density) +
    norm_func(Group_Size) +
    norm_func(Home_Range) + 
    resHR_norm_func(Residual_HomeRange),# Residual Home range has to be normalized differently bc it's residual (can't log negatives)
    # log10(ActivityCycle) +  
    # log10(HabitatBreadth) +
    # log10(ATWP_perForaging) + 
    # log10(ATWP_DR),
  lambda='ML'
  )  
  # + HR_filled + GroupSize_filled + ATWP_DR + ATWP_perForaging + ATWP_perTravel
  pgls_eco_summary <- summary(pgls_eco_raw)
  
  
  # Variable selection  using phylostep
  # best.model <- phylolm::phylostep(data = dat$data, 
  #                                  phy = dat$phy, 
  #                                  direction = "both", 
  #                                  model = "BM",
  #                                  formula = pgls_eco_raw$formula)
  # new.model <- pgls(data = dat, formula = eval(parse(text=best.model$formula)))
  # summary(new.model)
  if(variable_selection){
     optimized.model <- run_eco.variable_selection(pgls_function.out, 
                                                pgls_eco_raw, 
                                                eco_vars,
                                                res_norm_func, 
                                                resHR_norm_func,
                                                norm_func)
  } else {optimized.model <- NULL }
 
  print(pgls_eco_summary)
  output <- list(pgls_eco_summary=pgls_eco_summary, 
                 pgls_eco_raw=pgls_eco_raw,
                 optimized.model=optimized.model)
  return(output)
}
```


## Eco plot functions

```{r Eco plot functions} 
# Plot eco variable function
eco_plots <- function(pgls_function.out, 
                      eco_var, 
                      run_eco.out, 
                      figLabel="", 
                      res_norm_func=function(x){x},
                      norm_func_plot=log10,
                      x_prefix=""){
  
    pgls_eco_sum <- run_eco.out$pgls_eco_summary
    response_var <- pgls_function.out$response_var
    data <- pgls_function.out$new_comp_data$data
    brain <- data[,paste(response_var,"res",sep="_")] 
    eco_x <-data[,paste(eco_var)] 
    eco.row <- pgls_eco_sum$coefficients[grep(eco_var, rownames(pgls_eco_sum$coefficients)), ]
      slope <- as.numeric(eco.row[["Estimate"]])
      p <- as.numeric(eco.row[["Pr(>|t|)"]])
  #  Initialize plot 
  # if(eco_var=="Residual_HomeRange"){norm_func <- res_norm_func}
   eco_plot <- ggplot(data, aes(y=res_norm_func(brain), x=norm_func_plot(eco_x), 
                                fill=Clade, color=Clade, shape=Clade))  + 
     theme_classic() + ggtitle(figLabel)
  # Conditionally add regression line
  if ( p<=0.05){ 
      eco_plot <- eco_plot +  
      geom_smooth(na.rm=F, inherit.aes=F, method="lm", alpha = .15, 
                  data=data, aes(y=res_norm_func(brain), x=norm_func_plot(eco_x))) 
  } 
  # Conditionally add legend
  if (eco_var == "Residual_HomeRange"){
    eco_plot =  eco_plot + theme(legend.position=c(1.75, 0.35), legend.background = element_rect(fill="whitesmoke"),
                                 plot.title = element_text(hjust= -.1))
  } else{
    eco_plot = eco_plot + theme( legend.position="none",
                                 plot.title = element_text(hjust= -.1))
  }    
  # Add remaining featuress
  eco_plot <- eco_plot + geom_point(size=3, alpha=.7) +  
         xlab(paste(x_prefix,eco_var,sep="")) + ylab(paste(response_var,"res",sep="_")) +   
         scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + scale_fill_manual(values=colors) 

 
  #geom_abline(slope=SLOPE, intercept=INTERCEPT)
    #assign(paste("eco_plot",paste(response_var,"res",sep="_"), "vs",eco_var,sep="."),eco_plot) 
      return(eco_plot) # %>% ggplotly()
  }
  
# MAKE MORE FLEXIBLE (dependent on which model variables actually used) 
all_eco_plots <- function(pgls_function.out, 
                          norm_func_plot = log_scale,
                          resHR_norm_func_plot=scale_
                          ){
  eco=run_eco(pgls_function.out = pgls_function.out)
  print(eco$pgls_eco_summary)
  p_f.o <-pgls_function.out
  # e0 <- eco_plots(eco_var = "Nocturnality", run_eco.out=eco)
  e1 <-eco_plots(pgls_function.out = p_f.o, eco_var = "Diet_Breadth", run_eco.out=eco, figLabel="(a)", 
                 norm_func_plot = resHR_norm_func_plot)
  #e2 <-eco_plots(p_f.o, eco_var = "HomeRange_km2", run_eco.out=eco)
  #e3 <-eco_plots(p_f.o, eco_var = "HomeRange_Indiv_km2", run_eco.out=eco)
  e4 <-eco_plots(p_f.o, eco_var = "Population_Density", run_eco.out=eco, figLabel="(b)", norm_func_plot = norm_func_plot)
  # e5 <-eco_plots(p_f.o, eco_var = "SocialGrpSize", run_eco.out=eco)
  # e6 <-eco_plots(p_f.o, eco_var = "ATWP_DR", run_eco.out=eco)
  # e7 <-eco_plots(p_f.o, eco_var = "ATWP_perForaging", run_eco.out=eco)
  #e8 <-eco_plots(p_f.o, eco_var = "ATWP_GroupSize", run_eco.out=eco)
  #e9 <-eco_plots(p_f.o, eco_var = "ATWP_HR", run_eco.out=eco)
  #e10 <-eco_plots(p_f.o, eco_var = "ATWP_perTravel", run_eco.out=eco)
  e11 <-eco_plots(p_f.o, eco_var = "Group_Size", run_eco.out=eco, figLabel="(c)", norm_func_plot = norm_func_plot)
  e12 <-eco_plots(p_f.o, eco_var = "Home_Range", run_eco.out=eco, figLabel="(d)", norm_func_plot = norm_func_plot)
  e13 <-eco_plots(p_f.o, eco_var = "Residual_HomeRange", run_eco.out=eco, figLabel="(e)", 
                  norm_func_plot = resHR_norm_func_plot)
  
  # Additional variables
  # e14 <-eco_plots(p_f.o, eco_var = "ActivityCycle", run_eco.out=eco, figLabel="(f)")
  # e15 <-eco_plots(p_f.o, eco_var = "HabitatBreadth", run_eco.out=eco, figLabel="(g)") 
  # e16 <-eco_plots(p_f.o, eco_var = "ATWP_perForaging", run_eco.out=eco, figLabel="(h)") 
  # e17 <-eco_plots(p_f.o, eco_var = "ATWP_DR", run_eco.out=eco, figLabel="(i)") 
  
    
  mp <- gridExtra::grid.arrange(e1,e4,e11,e12,e13, 
                          # e14,e15,e16,e17, 
                          ncol=3, 
                           top = textGrob(paste("Ecological Plots:",
                                                paste(pgls_function.out$response_var,"res",sep="_")),
                                          gp=gpar(fontsize=18)))
  return(mp)
                         #e6,e7,e10
    # mod_vars <- list(eco_vars=c("DietBreadth","PopulationDensity_n_km2","GroupSize_filled", "Home_Range", "Residual_HomeRange")) 
  # subplot(e1,e4,e11,e12,e13, nrows = 2, shareY=F, shareX=F, titleX = T, titleY=T, margin = 0.1) 

  # fluidPage(
  #   fluidRow( column(4, e1), column(4, e4), column(4, e11)),
  #   fluidRow( column(4, e12), column(8, e13) )
  # )
  #plot_grid( plotlist = list(e1,e4,e11,e12,e13), labels = list()) 
}
```


# Plots & Analyses

## Calculate Residual_HomeRange

```{r Calculate Residual_HomeRange} 
suffix = " (log mm3)"#expression(" (log"~mm^3~')')
 
# HR_phylo <- contMap(HP$phy, HP$data$Home_Range, type="phylogram", plot=F)
# plot(HR_phylo, leg.txt="Home Range",lwd=3, fsize=c(0.5, 1)) 
  
# We know HomeRange scales with BodySize, so perhaps using "extra" HR relative to body weight is a more relavent measure.
  ## Alternatively could include BodySize as covariate in eco model.
# HP$data$Body.Size <- log10(HP$data$Body_Size)
HR.body_res <- pgls_function(HP, "Home_Range", "Body_Size", 
                             norm_func = log_scale,
                             x_prefix = "", x_suffix = " (g)", 
                             y_prefix = "", y_suffix = " (km2)") 

# Put back into Comp.data
HP$data$Residual_HomeRange <- HR.body_res$residuals
```

## Retrocommissural Hippocampus

### retroHP vs. Total Brain Vol

- For eco analysis, use ``scale_log`` (or else PGLS will throw a "reciprocal" error).  
- For eco plot, just plot ``log10`` of each ecological varialbe (except *Residual_HomeRange*, which is already normal), and the untransformed y-axis residuals (which are also already normal).

```{r retroHP vs. Total Brain Vol}
# HP vs. Brain
HP.Brain_res <- pgls_function(HP, "retroHP", "Brain__HP", y_suffix = suffix, x_suffix = suffix) 
run_eco(pgls_function.out = HP.Brain_res)
all_eco_plots(pgls_function.out = HP.Brain_res)
```

### __***retroHP vs. Medulla***__

```{r retroHP vs. Medulla}
# HP vs. Medulla
HP.Med_res <- pgls_function(HP, "retroHP", "Medulla", y_suffix = suffix, x_suffix = suffix)
retroHP.Medulla <- run_eco(HP.Med_res)
all_eco_plots(HP.Med_res)

# Hyp: Relationship with overall brain vol?
#pgls_function(HP, "retroHP_res", "BrainVol")
```

#### retroHP.vs.Medulla Phenogram

```{r retroHP.vs.Medulla Phenogram, fig.height=8,fig.width=7}
HPres <- HP.Med_res$new_comp_data$data$retroHP_res
names(HPres) <- HP.Med_res$new_comp_data$phy$tip.label
  ## Get colors from ggplot defaults
  gg_color_hue <- function(n) {
    hues = seq(15, 375, length = n + 1)
    hcl(h = hues, l = 65, c = 100)[1:n]
  }
  colz = gg_color_hue(6)

phenogram(HP.Med_res$new_comp_data$phy, x=HPres, spread.labels=T, 
          fsize=.75, colors=colz, axes=c(time, HPres),
          main="Phenogram:\n retroHP vs. medulla PGLS residuals", 
          xlab="Mya", ylab="retroHP residuals") #spread.cost = c(1, 0)
## Color branches by clades??
```

#### retroHP_res vs. BrainVol

```{r retroHP_res vs. BrainVol}
#Hyp: As brain size scales up, does retroHP_res becomes proportionally smaller?
HPres.Brain_res <- pgls_function(HP, "retroHP_res", "log10(Brain__HP)", y_suffix = "", 
                                 norm_func = function(x){x})
retroHP.Brain <- run_eco(pgls_function.out = HPres.Brain_res)
all_eco_plots(HPres.Brain_res) 
#Result: No relationship.
```

#### retroHP_res vs. Encephalization

```{r retroHP_res vs. Encephalization}
#Hyp: As Encephalization scales up, does retroHP_res become proportionally smaller?
Enceph_res <- pgls_function(HP, "BrainVol", "Body__Brain", y_suffix = suffix, x_suffix = suffix)
#Don't normalize residuals 
pgls_function(HP, "retroHP_res", "BrainVol_res", x_suffix = "", y_suffix = "")
run_eco(HPres.Enceph_res)
all_eco_plots(HPres.Enceph_res)
#Result: No relationship.
```


## Entorhinal Cortex*

* Entorhinal Cortex*: "Included are entorhinal, perirhinal, presubicular and parasubicular cortices and the underlying white matter. These cortices are characterized by the presence of one or several, almost cell-free layers of sublayers."
+ Stephan, H. et al. (1981) New and revised data on volume of brain structures in insectivores and primates. Folia primatol. 35, 1–29

### Entorhinal Cortex* vs. Medulla

NOTE: Using EC.med instead of EC.HP 1) flips the axes of the pPCA, and 2) gets rid of selective regime shift in human lineage. This is because humans are more divergent in EC.med (due to our massively enlarged neocotex, which EC is part of).


```{r Entorhinal Cortex* vs. Medulla}
EC.Med_res <-pgls_function(HP, "EC", "Medulla",  y_suffix = suffix, x_suffix = suffix)
run_eco(pgls_function.out = EC.Med_res)
all_eco_plots(EC.Med_res)
```

### Entorhinal Cortex* vs. BrainVol

- Happ-only eco model approaching sig (p=~0.06) after normalizing and using phylostep for variable selection.

```{r Entorhinal Cortex* vs. BrainVol, eval=F}
EC.Brain_res <-pgls_function(HP, "EC", "Brain__EC", y_suffix = suffix, x_suffix = suffix)
run_eco(EC.Brain_res)
all_eco_plots(EC.Brain_res)
```

### Entorhinal Cortex* vs. retroHP

**Note:: Used as final EC**

- Happ-only eco model sig (p=~0.01) after normalizing and using phylostep for variable selection.

```{r Entorhinal Cortex* vs. retroHP}
EC.HP_res <-pgls_function(HP, "EC", "retroHP", norm_func = log10,
                           y_suffix = suffix, x_suffix = suffix)
EC.retroHP <- run_eco(EC.HP_res)
all_eco_plots(EC.HP_res)
#pgls_function(HP, "EC_res", "BrainVol")
```

## HP Fibers

### HP Fibers vs. Medulla

```{r HP Fibers vs. Medulla}
fibers.Med_res <-pgls_function(HP, "Fibers", "Medulla", norm_func = log10,
                                y_suffix = suffix, x_suffix = suffix)
run_eco(fibers.Med_res)
all_eco_plots(fibers.Med_res)
```

### HP Fibers vs. BrainVol

```{r HP Fibers vs. BrainVol}
fibers.brain_res <-pgls_function(HP, "Fibers", "Brain__Fib", norm_func = log10,
                                  y_suffix = suffix, x_suffix = suffix)
run_eco(fibers.brain_res)
all_eco_plots(fibers.brain_res)
```

### HP Fibers vs. retroHP + EC * * *

Result: No relationship
```{r HP Fibers vs. retroHP + EC}
fibers.HP_res <-pgls_function(HP, "Fibers", "retroHP.EC", lambda =  0.001, norm_func = log10,
                               y_suffix = suffix, x_suffix = suffix)
Fibers.retroHP.EC <- run_eco(fibers.HP_res)
all_eco_plots(fibers.HP_res)
#pgls_function(HP, "Fibers_res", "BrainVol") 
```

## Subfields

### DentateGyrus vs. retroHP

Result: No relationship
```{r DentateGyrus vs. retroHP}
DG.HP_res <- pgls_function(HP, "DG", "retroHP__DG", norm_func = log10,
                            y_suffix = suffix, x_suffix = suffix)
DG.retroHP <- run_eco(DG.HP_res)
all_eco_plots(DG.HP_res)
#pgls_function(HP, "DG_res", "BrainVol")
```

### FasciaDentata vs. retroHP
 
Result: No relationship
```{r FasciaDentata vs. retroHP}
FD.HP_res <- pgls_function(HP, "FD", "retroHP__FD", norm_func = log10,
                            y_suffix = suffix, x_suffix = suffix)
FD.retroHP <- run_eco(FD.HP_res, variable_selection = T)
mp <- all_eco_plots(FD.HP_res)
ggsave("./Figures/FINAL/Fig4.png", mp, dpi=400, height=6, width=9)


# improve model
dat <- nonNA_subset(HP = HP, eco_vars = "Residual_HomeRange", response_var = "FD_res")
model <- caper::pgls(data = dat, formula = FD_res ~ Residual_HomeRange, lambda = "ML") 
summary(model)
#pgls_function(HP, "FD_res", "BrainVol")
```

### Hilus vs. retroHP
 
Result: No relationship
```{r Hilus vs. retroHP}
Hilus.HP_res <- pgls_function(HP, "Hilus", "retroHP__Hil", norm_func = log10,
                               y_suffix = suffix, x_suffix = suffix)
Hilus.retroHP <- run_eco(Hilus.HP_res, norm_func = logTrans_missing)
all_eco_plots(Hilus.HP_res)
#pgls_function(HP, "Hilus_res", "BrainVol")
```

### CA2/3 vs. retroHP

Result: Positive relationship between CA2.3_res and retroHP_res.
```{r CA2/3 vs. retroHP}
CA2.3.HP_res <- pgls_function(HP, "CA2.3", "retroHP__CA2.3", norm_func = log10,
                               y_suffix = suffix, x_suffix = suffix)
CA2.3.retroHP <- run_eco(CA2.3.HP_res)
all_eco_plots(CA2.3.HP_res)
#pgls_function(HP, "CA2.3_res", "BrainVol") # No relationship
#pgls_function(HP, "CA2.3_res", "retroHP") # No relationship
#pgls_function(HP, "CA2.3_res", "retroHP_res") # Signficant***
```

### CA3 vs. retroHP

Result: CA3_res gets smaller as retroHP_res gets larger. Humans deviate drastically from this trend by having an unusually large CA3_res.
```{r CA3 vs. retroHP}
CA3.HP_res <- pgls_function(HP, "CA3", "retroHP__CA3", norm_func = log10,
                             y_suffix = suffix, x_suffix = suffix)
CA3.retroHP <- run_eco(CA3.HP_res)
all_eco_plots(CA3.HP_res)

#pgls_function(HP, "CA3_res", "BrainVol")# No relationship
#pgls_function(HP, "CA3_res", "retroHP")# No relationship
#pgls_function(HP, "CA3_res", "retroHP_res") # Approaching sig **
```

### CA2 vs. retroHP

Result: CA2_res gets smaller has retroHP_res gets larger. Humans deviate from this by having somewhat larger CA2_res than expected.
```{r CA2 vs. retroHP}
CA2.HP_res <- pgls_function(HP, "CA2", "retroHP__CA2", norm_func = log10,
                             y_suffix = suffix, x_suffix = suffix)
CA2.retroHP <- run_eco(CA2.HP_res)
all_eco_plots(CA2.HP_res) 
#pgls_function(HP, "CA2_res", "BrainVol")# No relationship
#pgls_function(HP, "CA2_res", "retroHP")# No relationship
#pgls_function(HP, "CA2_res", "retroHP_res") # Significant *** 
```

### CA1 vs. retroHP

Result: No relationship.
```{r CA1 vs. retroHP}
CA1.HP_res <- pgls_function(HP, "CA1", "retroHP__CA1", y_suffix = suffix, x_suffix = suffix)
CA1.retroHP <- run_eco(pgls_function.out = CA1.HP_res, variable_selection = T)
ca1.plot <- all_eco_plots(CA1.HP_res)
print(ca1.plot)
ggsave("./Figures/FINAL/Fig5.png", ca1.plot, dpi=400, height=6, width=9)

 
dat <- nonNA_subset(HP = HP, eco_vars = "Residual_HomeRange", response_var = "CA1_res")
model <- caper::pgls(data = dat, formula = CA1_res ~ Residual_HomeRange, lambda = "ML") 
summary(model)

CA1.HP_res <- pgls_function(dat, "CA1_res", "Residual_HomeRange", 
                            y_suffix = suffix, x_suffix = suffix, norm_func = function(x){x})

CA1.HP_res$xy_plot + 
  geom_point(data = HP$data["Homo_sapiens",], aes(x=Residual_HomeRange, y=CA1_res), 
             shape="2", color="goldenrod2", size=4)
 
# Phylogram
res_name <- "Residual_HomeRange"

cont_map <- function(comp_data, res_name){
  # comp_data <- HP
  HPres_data <- as.matrix(comp_data$data[res_name])
  names(HPres_data) <- comp_data$phy$tip.label
  res_phylo_plot <- contMap(comp_data$phy, 
                            HPres_data, 
                            type="phylogram", 
                            plot=F, fsize=c(0.5, 1)) 
  
 # fancyTree(tree = comp_data$phy, 
  #           x=HPres_data, 
  #           type="contmap",
  #           fsize=c(0.5, 1))
    ## simulate correlated trait data
  # http://blog.phytools.org/2013/09/plotting-facing-trees-using-phytools.html 
  return(res_phylo_plot)
}

fsize <- c(.5, 1)

plot.new()
par(mfrow=c(1,1)) 

res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name) 
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize) 
png("./Figures/pheno1.png", res = 100, height=10)

dev.off()

res_name <- "FD_res"
cp1 <- cont_map(HP, res_name)
plot(cp1, leg.txt=paste(res_name),lwd=3, fsize=fsize)

res_name <- "CA1_res"
cp2 <- cont_map(HP, res_name)
plot(cp2, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()


 
 

# HR.HRres <- pgls_function(HP, "Home_Range", "Residual_HomeRange", norm_func = function(x){x}, y_suffix = " (km2)", x_suffix = "")
# HR.HRres.log <- pgls_function(HP, "log10(Home_Range)", "Residual_HomeRange", norm_func = function(x){x}, y_suffix = " (km2)", x_suffix = "")
# 
# HR.HRres$xy_plot + labs(title="Untransformed") +
# HR.HRres.log$xy_plot + labs(title="Log transformed")
# 
# 

# ALternative methods for selective regime shift analysis 
# shifts <- phylolm::OUshifts(y=dat$data$CA1_res, phy = dat$phy, nmax=5)
# phylolm::plot.OUshifts(shifts)

# logTrans_missing(CA1.HP_res$new_comp_data$data$CA1_res)
# Reduced model with only the sig eco vars
pgls.plot_predictions <- function(pgls.object, x="Residual_HomeRange", y="Home_Range"){
  # pgls.object <- CA1.HP_res
  comp_data <- pgls.object$new_comp_data
  model <- caper::pgls(data=comp_data, formula = eval(parse(text=paste(y,"~",x))))
  pgls_summary <- summary(model) 
  pgls.xy_plot(comp_data, pgls_summary, response_var = y, predictor_var = x)
 
  
  x.values <- as.numeric(comp_data$data[[x]])
  x.smooth<- seq(min(x.values),max(x.values),.01)
  setNames(list(x.smooth), x)
  exp.predict <- exp(predict(model, list(Residual_HomeRange = x.smooth)))
  exp.line <- data.frame(x=x.smooth, y=exp.predict)
  pgls.object$xy_plot + geom_line(data = exp.line, aes(x,y), inherit.aes = F) 
}


```

### Subiculum vs. retroHP

Result: No relationship
```{r Subiculum vs. retroHP}
Sub.HP_res <- pgls_function(HP, "Sub", "retroHP__Sub", y_suffix = suffix, x_suffix = suffix)
Sub.retroHP <- run_eco(Sub.HP_res)
all_eco_plots(Sub.HP_res) 
#pgls_function(HP, "Sub_res", "BrainVol")
```


# Summaries
## All neuroanatomical regression plots
### Scaling diagnostic multiplots
```{r Scaling diagnostic multiplots}
grid.arrange(HP.Med_res$xy_plot + ggtitle("(a)") + 
               theme(plot.title = element_text(hjust= -.13, size = 18)), 
             HP.Brain_res$xy_plot+ ggtitle("(b)") + 
               theme(plot.title = element_text(hjust= -.13, size = 18)), 
             HPres.Enceph_res$xy_plot + ggtitle("(c)") + 
               theme(plot.title = element_text(hjust= -.13, size = 18)),
             HPres.Brain_res$xy_plot + ggtitle("(d)") + 
               theme(plot.title = element_text(hjust= -.13, size = 18)), 
             HR.body_res$xy_plot + ggtitle("(e)") +
                theme(legend.position=c(1.3, 0.35), legend.background = element_rect(fill="whitesmoke"),
                      plot.title = element_text(hjust= -.13, size = 18)),
             ncol = 3, 
             top = textGrob("Scaling Diagnostics", gp=gpar(fontsize=18)))

# HR.body_res$pgls_summary
```

### EC & Fibers multiplots
```{r EC & Fibers multiplots}
grid.arrange(EC.HP_res$xy_plot, 
             EC.Med_res$xy_plot, 
             EC.Brain_res$xy_plot,
             fibers.HP_res$xy_plot, 
             fibers.Med_res$xy_plot, 
             fibers.brain_res$xy_plot, 
             ncol=3, top=textGrob("EC & Fiber Regressions", gp=gpar(fontsize=18)))
```

### retroHP subfield multiplots
```{r retroHP subfield multiplots} 
mp <- cowplot::plot_grid(# Scaling Diagnostics
                   HP.Med_res$xy_plot + theme_classic() + 
                     theme(legend.position="none", plot.background = element_rect(fill="whitesmoke")), 
                   HP.Brain_res$xy_plot + theme_classic() + 
                     theme(legend.position="none", plot.background = element_rect(fill="whitesmoke")), 
                   HPres.Brain_res$xy_plot + theme_classic() + 
                     theme(legend.position="none", plot.background = element_rect(fill="whitesmoke")), 
                   HPres.Enceph_res$xy_plot + theme_classic() +
                     theme(legend.position="none", plot.background = element_rect(fill="whitesmoke")), 
                   HR.body_res$xy_plot + theme_classic() + 
                     theme(legend.position="none", plot.background = element_rect(fill="whitesmoke")), 
                   # Subregional Plots
                   EC.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   DG.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   FD.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   Hilus.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   CA2.3.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   CA3.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   CA2.HP_res$xy_plot + theme_classic() + theme(legend.position="none"),
                   CA1.HP_res$xy_plot + theme_classic() + theme(legend.position="none"), 
                   Sub.HP_res$xy_plot + theme_classic() +  theme(legend.position=c(1.5, 0.5), 
                                                            legend.background = element_rect(fill="whitesmoke"),
                                                            plot.title = element_text(hjust= -.13, size = 18)),
                   labels = paste0("(",letters,")"),
                   label_x = .15,
                   ncol=5)
print(mp)
ggsave("./Figures/FINAL/Fig3.png", mp, dpi = 400, height=8, width=15)
```

## Subfield Proportions {.tabset .tabset-fade .tabset-pills}

Plot the volumetric proportions of each subregion of the hippocampal complex.

### Subfield Proportions

```{r Subfield Proportions} 
subregion_proportions <- function(HP_data, subregions){ 
  subfield_props <- HP_data %>% 
    subset(select=c("Clade", subregions)) %>%   
    mutate(HP_complex = rowSums( .[subregions])) %>% 
    mutate_at(vars(-Clade), funs(./ HP_complex*100)) %>%
    # group_by(Clade) %>% 
    summarise_at(vars(-Clade), funs(mean(., na.rm=TRUE))) %>% 
    subset(select=-HP_complex) %>%
    melt(id.vars = c("Clade"),  variable.name = "Subregion", value.name = "% Total Volume") %>% arrange(Clade)
   subfield_props
  
  bp <- ggplot(subfield_props, aes(x=Clade, y = `% Total Volume`, fill = Subregion)) + 
    geom_bar(position = "fill", stat="identity") + labs(x="Clade", y="% Total Volume") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  bp %>% ggplotly()
} 
# Absolute
subregion_proportions(HP_data = HP$data, 
                      subregions = c("Sub", "CA1", "CA2", "CA3", "Hilus", "FD", "EC","Fibers")) 
```

### Residual Proportions 

```{r Residual Proportions}
residual_proportions <- function (HP_data, residuals){
  df <- HP_data %>%  subset(select=c("Clade", residuals)) %>% 
             group_by(Clade) %>% summarise_at(vars(-Clade), funs(mean(., na.rm=TRUE))) %>% 
    melt(id.vars = c("Clade"),  variable.name = "Subregion", value.name = "Subregion Residuals") %>% arrange(Clade)
  
  bp <- ggplot(df, aes(x=Clade, y = `Subregion Residuals`, fill = Subregion)) + 
    geom_bar(position = "fill", stat="identity") + labs(x="Clade", y="Subregion Residuals") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  bp %>% ggplotly()
}   

# # Residuals
residual_proportions(HP$data, residuals=c("retroHP__Sub", "retroHP__CA1", "retroHP__CA2",
                                           "retroHP__CA3", "retroHP__Hil", "retroHP__FD") )
```


## Summary stats

### PGLS: Neuroanatomy regression stats

```{r PGLS: Neuroanatomy regression stats}
#pgls_function.out <- Sub.HP_res
extract_neuro <- function(pgls_function.out){
  # Get comparison name
  var_name <- paste(pgls_function.out$response_var,"res",sep="_")
  pgls_summary <- pgls_function.out$pgls_summary
  # Get model p.value
    f <- pgls_summary$fstatistic
    model.p <- round( pf(f[1],f[2],f[3],lower.tail=F) ,3)
  # Get other measures
  Multiple.R2 <- round(pgls_summary$r.squared,3)
  Adj.R2 <- round(pgls_summary$adj.r.squared[1],3)
  F.stat <- round(pgls_summary$fstatistic[1],3)
  Lambda <- round(pgls_summary$param[2],3)

  output <- list(Multiple.R2=Multiple.R2, Adj.R2=Adj.R2, F.stat=F.stat, Lambda=Lambda,
                 model.pval=model.p)
  return(output)
}

neuro.summ <- rbind(
      data.frame(Residual="retroHP.med_res", extract_neuro(HP.Med_res)),
      data.frame(Residual="retroHP.brain_res", extract_neuro(HP.Brain_res)),
      data.frame(Residual="EC.HP_res", extract_neuro(EC.HP_res)),
      data.frame(Residual="Fiber.HP_res", extract_neuro(fibers.HP_res)),
      data.frame(Residual="DG.HP_res", extract_neuro(DG.HP_res)),
      data.frame(Residual="FD.HP_res", extract_neuro(FD.HP_res)),
      data.frame(Residual="Hilus.HP_res", extract_neuro(Hilus.HP_res)),
      data.frame(Residual="CA2.3.HP_res", extract_neuro(CA2.3.HP_res)),
      data.frame(Residual="CA3.HP_res", extract_neuro(CA3.HP_res)),
      data.frame(Residual="CA2.HP_res", extract_neuro(CA2.HP_res)),
      data.frame(Residual="CA1.HP_res", extract_neuro(CA1.HP_res)),
      data.frame(Residual="Sub.HP_res", extract_neuro(Sub.HP_res))
      )
# Bonferonni correction for multiple tests
alpha <- round(0.06/11, 4)
paste("Alpha = ", alpha)
neuro.summ$Sig <-  ifelse(neuro.summ$model.pval<0.05 & neuro.summ$model.pval>0.01, "*",
              ifelse(neuro.summ$model.pval<0.01 & neuro.summ$model.pval>alpha,"**", 
              ifelse(neuro.summ$model.pval<alpha,"***", "")))
row.names(neuro.summ)<-NULL
colnames(neuro.summ) <- c("Residual Name",  "Multiple R<sup>2<sup>","Adj. R<sup>2<sup>", 
                          "F stat.", "λ", "Model p-value", "Sig.")

createDT <- function(dt, caption, save=F, scrollY=700){
  if (save!=F){ 
    dir.create("./Tables", showWarnings=F)
    write.csv(dt, file=file.path("Tables",paste(save,"csv",sep=".")), quote = F, row.names = F)  
  }
    datatable( dt, caption=caption,
      extensions = 'Buttons', escape = F,
      options = list( dom = 'Bfrtip', 
                      buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), 
                      scrollY = scrollY, paging = F

    )
  ) 
  
}

createDT(neuro.summ, caption="Neuroanatomical PGLS Regressions Summary", save="neuro.results") 

```

### PGLS: HP-eco stats

```{r PGLS: Optimized HP-eco stats}
optimized.models <- list(retroHP.Medulla, 
  retroHP.Brain, 
  EC.retroHP, 
  Fibers.retroHP.EC, 
  DG.retroHP, 
  FD.retroHP, 
  Hilus.retroHP,
  CA2.3.retroHP,
  CA3.retroHP,
  CA2.retroHP, 
  CA1.retroHP,
  Sub.retroHP)

# mod <- optimized.models[[1]]
optimized.df <- lapply(optimized.models, function(mod){  
  model <- mod$optimized.model
  pgls_eco_summary <- summary(model)
  
  pgls_eco_summary
  digits <- 3
  sjstats::boot_p(model$bootstrap) 
  sjstats::p_value(model) 
   
  AIC <- model$aic %>% round(digits)
  LogLik <- phylolm::logLik.phylolm(model)$logLik %>% round(digits)
  lambda <- model$optpar  %>% round(digits)
  
  
  
  # Get model p.value 
  model.p <- round( pf(f[1],f[2],f[3],lower.tail=F) ,3)
  # Get other measures
  Multiple.R2 <- round(pgls_eco_summary$r.squared,3)
  Adj.R2 <- round(pgls_eco_summary$adj.r.squared[1],3)
  F.stat <- round(pgls_eco_summary$fstatistic[1],3)
  Lambda <- round(pgls_eco_summary$param[2],3) 
  predictors <- gsub("norm_func|resHR_norm_func|res_norm_func|[)]|[(]","",
                     paste(mod$optimized.model$varNames[-1], collapse=", ")) 
  
  output <- list(`Response`=mod$optimized.model$varNames[[1]],
                 `Optimized Predictors`=predictors,
                 Multiple.R2=Multiple.R2, 
                 Adj.R2=Adj.R2, 
                 F.stat=F.stat, 
                 Lambda=Lambda,
                 model.pval=model.p) 
  return(output)
  }) %>% data.table::rbindlist()

alpha <- round(0.05/12, 4)
paste("Alpha = ", alpha)optimized.df
optimized.df$Sig <-  ifelse(optimized.df$model.pval<0.05 & optimized.df$model.pval>0.01, "*",
              ifelse(optimized.df$model.pval<0.01 & optimized.df$model.pval>alpha,"**", 
              ifelse(optimized.df$model.pval<alpha,"***", "")))
```


```{r PGLS: HP-eco stats}
#pgls_function.out <- Sub.HP_res
extract_eco <- function(pgls_function.out){
  # Get comparison name
  var_name <- paste(pgls_function.out$response_var,"res",sep="_")
  run_eco.out <-run_eco(pgls_function.out)
  pgls_eco_summary <-run_eco.out$pgls_eco_summary
  # Get model p.value
    f <- pgls_eco_summary$fstatistic
    model.p <- round( pf(f[1],f[2],f[3],lower.tail=F) ,3)
  # Get other measures
  Multiple.R2 <- round(pgls_eco_summary$r.squared,3)
  Adj.R2 <- round(pgls_eco_summary$adj.r.squared[1],3)
  F.stat <- round(pgls_eco_summary$fstatistic[1],3)
  Lambda <- round(pgls_eco_summary$param[2],3)
  output <- list(Multiple.R2=Multiple.R2, Adj.R2=Adj.R2, F.stat=F.stat, Lambda=Lambda,
                 model.pval=model.p)
  return(output)
}
eco.summ <- rbind( 
      data.frame(Residual="retroHP.med_res", extract_eco(HP.Med_res)),
      data.frame(Residual="retroHP.brain_res", extract_eco(HP.Brain_res)),
      data.frame(Residual="EC.HP_res", extract_eco(EC.HP_res)),
      data.frame(Residual="Fiber.HP_res", extract_eco(fibers.HP_res)),
      data.frame(Residual="DG.HP_res", extract_eco(DG.HP_res)),
      data.frame(Residual="FD.HP_res", extract_eco(FD.HP_res)),
      data.frame(Residual="Hilus.HP_res", extract_eco(Hilus.HP_res)),
      data.frame(Residual="CA2.3.HP_res", extract_eco(CA2.3.HP_res)),
      data.frame(Residual="CA3.HP_res", extract_eco(CA3.HP_res)),
      data.frame(Residual="CA2.HP_res", extract_eco(CA2.HP_res)),
      data.frame(Residual="CA1.HP_res", extract_eco(CA1.HP_res)),
      data.frame(Residual="Sub.HP_res", extract_eco(Sub.HP_res))
      )
# Bonferonni correction for multiple tests
alpha <- round(0.05/11, 4)
paste("Alpha = ", alpha)
eco.summ$Sig <-  ifelse(eco.summ$model.pval<0.05 & eco.summ$model.pval>0.01, "*",
              ifelse(eco.summ$model.pval<0.01 & eco.summ$model.pval>alpha,"**", 
              ifelse(eco.summ$model.pval<alpha,"***", "")))
row.names(eco.summ)<-NULL
colnames(eco.summ) <- c("Residual Name",  "Multiple R<sup>2<sup>","Adj. R<sup>2<sup>", 
                          "F stat.", "λ", "Model p-value", "Sig.")
createDT(eco.summ, caption="Ecological PGLS Regressions Summary", save = "eco.results")
```

# Multidimensional plots

## Set up interactive 3D plots & PCA data

```{r Set up interactive 3D plots & PCA data}
# Select which vars to put in PCA
## All residuals
all.vars <- c("retroHP_res", "EC_res", "Fibers_res", "DG_res", "FD_res", "Hilus_res", "CA2.3_res","CA3_res", "CA2_res", "CA1_res", "Sub_res")
## Simplified
vars <- c("EC_res", "Fibers_res", "FD_res", "Hilus_res","CA3_res", "CA2_res", "CA1_res", "Sub_res") 
Clade <- factor(HP$data$Clade, c("Homo sapiens","Apes","Old World Monkeys","New World Monkeys","Tarsiiformes","Strepsirrhines"), ordered=T)
plot.alldata <- dplyr::select(HP$data, one_of(all.vars)) %>% na.omit()
plot.data <- dplyr::select(HP$data, one_of(vars)) %>% na.omit()
```


## Non-phylogenetic PCA

```{r Non-phylogenetic PCA, fig.height=8, fig.width=10}
PCA <- prcomp(plot.data, scale=TRUE, center=TRUE)
pander::pander(PCA$rotation[,1:4]) # Show PC1-4
summPCA <- summary(PCA)

# Simple biplot
##biplot(PCA)
# Fancy biplot

ggbiplot(PCA, obs.scale=1, var.scale=1, groups=Clade, ellipse=T, circle=F, alpha=1) +
  scale_color_discrete(name = '') + theme_classic()


# Polygonal PCA
## Prep PCA data
grouped_omit <- na.omit(all)
Species <- rownames(HP$data)
  levels(Species) <- Species
DFforPlot <- data.frame(PC1=PCA$x[,1], PC2=PCA$x[,2], Species=Species)
DFforPlot["Clade"] <- HP$data$Clade
  Ape <- subset(DFforPlot, Clade=="Apes")
  NWM <- subset(DFforPlot, Clade=="New World Monkeys")
  OWM <- subset(DFforPlot, Clade=="Old World Monkeys")
  Stp <- subset(DFforPlot, Clade=="Strepsirrhines")
  Trs <- subset(DFforPlot, Clade=="Tarsiiformes")
## Plot PCA
ggplot(data=DFforPlot, aes(x=PC1, y=PC2, color=Clade, fill=Clade, shape=Clade, label=Species)) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  geom_polygon(aes(x=PC1, y=PC2), alpha=.3, data=Ape[chull(Ape$PC1, Ape$PC2),]) +  
  geom_polygon(aes(x=PC1, y=PC2), alpha=.3, data=NWM[chull(NWM$PC1, NWM$PC2),]) + 
  geom_polygon(aes(x=PC1, y=PC2), alpha=.3, data=OWM[chull(OWM$PC1, OWM$PC2),]) + 
  geom_polygon(aes(x=PC1, y=PC2), alpha=.3, data=Stp[chull(Stp$PC1, Stp$PC2),]) + 
  geom_polygon(aes(x=PC1, y=PC2), alpha=.3, data=Trs[chull(Trs$PC1, Trs$PC2),]) + 
  geom_point(size=3) + labs(title = "Residual Hippocampal Subregions PCA", x = paste("PC1","(", round(summPCA$importance[2,1]*100, 2), "%)"), y = paste("PC2","(", round(summPCA$importance[2,2]*100, 2), "%)")) +
  geom_text(aes(label=Species),vjust=-1,hjust=.5) +    
  scale_x_continuous(expand = c(.1, .1)) + scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + scale_fill_manual(values = colors) 
```

### 3D PCA
```{r 3D PCA, webgl=TRUE} 
plot3d(PCA$x[,1:3], col=as.numeric(as.factor(Clade)), size=1, type='s', main="3D PCA: Residual Hippocampal Subregions")
  # Add 3D vectors
  text3d(PCA$x[,1:3], texts=rownames(PCA$x), col=as.numeric(as.factor(Clade))) # points
  text3d(PCA$rotation[,1:3]*5, texts=rownames(PCA$rotation), col="purple3") # vectors
  coords <- NULL
  for (i in 1:nrow(PCA$rotation)) {
   coords <- rbind(coords, rbind(c(0,0,0),PCA$rotation[i,1:3]))
    }
  lines3d(coords*5, col="purple3", lwd=2)
  #legend3d("topright", legend=unique(Clade), col=as.numeric(unique(Clade)))
```



## Phylogenetic PCA

NOTE: By including DG_res and CA2.3_res, you're inflating the contribution of those regions (which, as it so happens, are extremely divergent in humans). So adding combinations of these measures multiple times overinflates how unique humans are.
```{r Phylogenetic PCA, fig.height=8, fig.width=10}
phy.PCA <- phyl.pca(HP$phy, plot.data, method="lambda")
pander::pander(phy.PCA$L[,1:4]) # Show PC1-4
phy.PCA$lambda # Show lambda
# Select which vars to put in PCA
summPCA <- summary(phy.PCA)
 
# Simple biplot
biplot(phy.PCA$S, 
       phy.PCA$L, 
       xlabs=rep("o", nrow(phy.PCA$S)),
       main = "Phylogenetic PCA: PC Axes", 
       col=c("blue","purple"))
# (ggbiplot doesn't recognize phyl.pca format)

# Polygonal PCA
## Get Propotion of Variance for each PC
POV <- diag(phy.PCA$Eval)/sum(phy.PCA$Eval)*100
## PC1 vs. PC2
grouped_omit <- na.omit(all)
Species <- rownames(HP$data)
  levels(Species) <- Species
DFforPlot <- data.frame(PC1=phy.PCA$S[,1], PC2=phy.PCA$S[,2], Species=Species)
DFforPlot["Clade"] <- HP$data$Clade
  Ape <- subset(DFforPlot, Clade=="Apes")
  NWM <- subset(DFforPlot, Clade=="New World Monkeys")
  OWM <- subset(DFforPlot, Clade=="Old World Monkeys")
  Stp <- subset(DFforPlot, Clade=="Strepsirrhines")
  Trs <- subset(DFforPlot, Clade=="Tarsiiformes")
## Plot PCA
PC12 <- ggplot(data=DFforPlot, aes(x=PC1, y=PC2, color=Clade, fill=Clade, shape=Clade, label=Species)) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle =element_text(hjust = 0.5)) +
  stat_ellipse(geom="polygon", level=0.95, alpha=0.2) +
  geom_point(size=3) + labs(title = "Phylogenetic PCA", 
                            subtitle = "Residual Hippocampal Subregions",
                            x = paste("PC1","(", round(POV[1], 2), "%)"), 
                            y = paste("PC2","(", round(POV[2], 2), "%)")) +
  scale_x_continuous(expand = c(.1, .1)) + 
  scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + 
  scale_fill_manual(values = colors) +
  geom_text_repel(aes(label=Species), box.padding=unit(0.35, "lines"), 
                  point.padding = unit(0.1, "lines"))
PC12# %>% ggplotly()

## PC3 vs. PC4
DFforPlot <- data.frame(PC3=phy.PCA$S[,3], PC4=phy.PCA$S[,4], Species=Species)
DFforPlot["Clade"] <- HP$data$Clade
  Ape <- subset(DFforPlot, Clade=="Apes")
  NWM <- subset(DFforPlot, Clade=="New World Monkeys")
  OWM <- subset(DFforPlot, Clade=="Old World Monkeys")
  Stp <- subset(DFforPlot, Clade=="Strepsirrhines")
  Trs <- subset(DFforPlot, Clade=="Tarsiiformes")
## Plot PCA 
PC34 <- ggplot(data=DFforPlot, aes(x=PC3, y=PC4, color=Clade, fill=Clade, shape=Clade, label=Species)) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  stat_ellipse(geom="polygon", level=0.95, alpha=0.2) +
  geom_point(size=3) + labs(title = "Phylogenetic PCA:", 
                            subtitle = "Residual Hippocampal Subregions", 
                            x = paste("PC3","(", round(POV[3], 2), "%)"), 
                            y = paste("PC4","(", round(POV[4], 2), "%)")) +
  scale_x_continuous(expand = c(.1, .1)) + 
  scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + 
  scale_fill_manual(values = colors) +
  geom_text_repel(aes(label=Species), box.padding=unit(0.35, "lines"), 
                  point.padding = unit(0.1, "lines"))
PC34 %>% ggplotly()
```


### 3D Phylogenetifc PCA
```{r 3D Phylogenetifc PCA, webgl=TRUE} 
plot3d(phy.PCA$S[,1:3], col=as.numeric(as.factor(Clade)), size=1, type='s', main="Phylogenetic 3D PCA: Residual Hippocampal Subregions")
  # Add 3D vectors
  #text3d(phy.PCA$S[,1:3],texts=rownames(phy.PCA$S),col=as.numeric(as.factor(Clade))) #points
  text3d(phy.PCA$L[,1:3], texts=rownames(phy.PCA$L), col="purple3") # vectors
  coords <- NULL
  for (i in 1:nrow(phy.PCA$L)) {
   coords <- rbind(coords, rbind(c(0,0,0),phy.PCA$L[i,1:3]))
    }
  lines3d(coords, col="purple3", lwd=2)
  #legend3d("topright", legend=unique(Clade), col=c("blue","turquoise","green","black","red","magenta"))
```

## Hierarchical Clustering + Heatmap {.tabset .tabset-fade .tabset-pills}

### Heatmaply

```{r Heatmaply, out.width="120%"}
Spectral <- grDevices::colorRampPalette(RColorBrewer::brewer.pal( length(orderedClades), "Spectral"))
heatmaply(plot.alldata, k_row = length(orderedClades), plot_method= "plotly", fontsize_row = 8, 
          row_side_colors =  list(Clade=Clade), row_side_palette = Spectral, 
          colorbar_xpos = 1.02,  key.title="Scaled Residuals") %>% 
          colorbar( len = 0.3,  tickfont = list(size = 8), titlefont = list(size = 12), which = 2)  
```

### Phylo.heatmap

```{r Phylo.heatmap, fig.height=7, fig.width=9.8}
# # Method 1: ggdendrogram
# hc <- hclust(dist(plot.alldata), "ave")
# ggdendrogram(hc, rotate=T) + labs(title="Multifactorial Cluster Dendrogram: \nResidual Hippocampal Subregions") + theme(plot.title = element_text(hjust = 0.5)) 
#+ geom_text(data=plot.data, aes(label=Clade, x=x, y=0, colour=Clade))

# Heatmaps


## Default
#heatmap( as.matrix(plot.alldata), main="Hierarchical Clustering Heatmap")
#heatmap.2
# heatmap.2(as.matrix(plot.alldata), srtCol=45) # ColSideColors=colz
## NMF
#NMF::aheatmap(as.matrix(plot.alldata))
## Ape
phylo.heatmap(HP$phy, plot.alldata, standardize = T) 
 
```

## Principal Components: Trees

### PC1 ContMap

```{r PC1 ContMap}
phy.PC1 <- phy.PCA$S[,1]
res_phylo_plot <- contMap(HP$phy, phy.PC1, type="phylogram", plot=FALSE)
  plot(res_phylo_plot, leg.txt="PC1",lwd=3, fsize=c(0.5, 1))
```

### PC1 Phenogram
```{r PC1 Phenogram, fig.height=8,fig.width=7}
phenogram(HP$phy, x=phy.PC1, spread.labels=T, fsize=.75, colors=colz, axes=c(time, phy.PC1), main="Phenogram:\n PC1", xlab="Mya", ylab="PC1") #spread.cost = c(1, 0), 
```

### PC2 ContMap
```{r PC2 ContMap}
phy.PC2 <- phy.PCA$S[,2]
res_phylo_plot <- contMap(HP$phy, phy.PC2, type="phylogram", plot=FALSE)
  plot(res_phylo_plot, leg.txt="PC2",lwd=3, fsize=c(0.5, 1))
```

### PC2 Phenogram
```{r PC2 Phenogram, fig.height=8,fig.width=7}
phenogram(HP$phy, x=phy.PC2, spread.labels=T, fsize=.75, colors=colz, axes=c(time, phy.PC2), main="Phenogram:\n PC2", xlab="Mya", ylab="PC2") #spread.cost = c(1, 0), 
```

## Principal Components: Scaling

### PC1 vs. retroHP_res

```{r PC1 vs. retroHP_res}
#Hyp: The hippocampus becomes increasingly reorganized as residual retroHP becomes larger.
HP$data$PC1 <- scales::rescale(phy.PC1,c(1,2))
#PC1.HP_res <- pgls_function(HP, "PC1", "retroHP") # No relationship
#PC1.Brain <- pgls_function(HP, "PC1", "BrainVol") # No relationship
PC1.HPres_res <- pgls_function(HP, "PC1", "retroHP_res", interact = F)
CA3.PC1_res <- pgls_function(HP, "CA3_res", "PC1") # Strong relationship, duh.

#Hyp: Since PC1 is largely driven by CA2.3_res, and PC1 correlates with retroHP_res, is there a relationship between retroHP_res and CA3_res?
CA3.HPres_res <- pgls_function(HP, "CA3_res", "retroHP_res") #Result: Approaching sig.
#Result: Significant positive relationship! But somewhat weak R2
```

#### PC2 vs. retroHP_res

```{r PC2 vs. retroHP_res}
#Hyp: The hippocampus becomes increasingly reorganized as residual retroHP becomes larger.
HP$data$PC2 <- scales::rescale(phy.PC2,c(1,2))
#PC2.HPres_res <- pgls_function(HP, "PC2", "retroHP") # No relationship
#PC2.HPres_res <- pgls_function(HP, "PC2", "BrainVol") # No relationship
#PC2.HPres_res <- pgls_function(HP, "PC2", "retroHP_res") # No relationship
PC2.HPres_res <- pgls_function(HP, "PC2", "Hilus_res") # Strong correlation, duh.
#Result: No relationship.
```

## Principal Component: Eco

### PC-eco functions

```{r PC-eco functions}
# Run eco PC
PCrun_eco <- function(comp_data, response_var){ 
  HP.noHomo <- subset(comp_data, subset=! Genus=="Homo")
  pgls_eco_raw <- pgls(data=HP.noHomo, formula=
    eval(parse(text=paste(response_var))) ~ 
                     DietBreadth + PopulationDensity_n_km2 +
                    GroupSize_filled + Home_Range + Residual_HomeRange,
                   lambda='ML') 
  pgls_eco_summary <- summary(pgls_eco_raw)

  output <- list(pgls_eco_summary=pgls_eco_summary, pgls_eco_raw=pgls_eco_raw, response_var=response_var, comp_data=comp_data)
  return(output)
}


# Plot PC1-eco
PCeco_plots <- function(eco_var, run_eco.out, figLabel){
  
    pgls_eco_summary <- run_eco.out$pgls_eco_summary
    response_var <- run_eco.out$response_var
    data <- run_eco.out$comp_data$data
    brain <-data[,paste(run_eco.out$response_var)]
    eco_x <- data[,paste(eco_var)] 
      slope <- as.numeric(pgls_eco_summary$coefficients[eco_var,1])
      p <- as.numeric(pgls_eco_summary$coefficients[eco_var,4])
      
  #  Initialize plot
  eco_plot <- ggplot(data, aes(y=brain, x=eco_x, fill=Clade, color=Clade, shape=Clade))  + 
    theme_classic() + ggtitle(figLabel)
  
  # Conditionally add regression line
  if ( p<=0.05){ 
    eco_plot <- eco_plot +  
      geom_smooth(na.rm=F, inherit.aes=F, method="lm", alpha = .15, data=data, aes(y=brain, x=eco_x))
  } 
  # Conditionally add legend
  if (eco_var == "Residual_HomeRange"){
    eco_plot =  eco_plot + theme(legend.position=c(1.75, 0.35), legend.background = element_rect(fill="whitesmoke"),
                                 plot.title = element_text(hjust= -.3))
  } else{
    eco_plot = eco_plot + theme( legend.position="none",
                                 plot.title = element_text(hjust= -.3))
  }    
  # Add remaining featuress
  eco_plot <- eco_plot + geom_point(size=3, alpha=.7) +  
         xlab(paste(eco_var,sep="")) + ylab(paste(response_var,"res",sep="_")) +   
         scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + scale_fill_manual(values=colors) 

 
  #geom_abline(slope=SLOPE, intercept=INTERCEPT)
    #assign(paste("eco_plot",paste(response_var,"res",sep="_"), "vs",eco_var,sep="."),eco_plot) 
      return(eco_plot) # %>% ggplotly()
  }
  
# MAKE MORE FLEXIBLE (dependent on which model variables actually used)
PCall_eco_plots <- function(comp_data, response_var){
  eco=PCrun_eco(comp_data, response_var)
  e1 <-PCeco_plots(eco_var = "DietBreadth", run_eco.out=eco, figLabel="A")
  e4 <-PCeco_plots(eco_var = "PopulationDensity_n_km2", run_eco.out=eco, figLabel="B")
  e11 <-PCeco_plots(eco_var = "GroupSize_filled", run_eco.out=eco, figLabel="C")
  e12 <-PCeco_plots(eco_var = "Home_Range", run_eco.out=eco, figLabel="D")
  e13 <-PCeco_plots(eco_var = "Residual_HomeRange", run_eco.out=eco, figLabel="E")

  # mod_vars <- list(eco_vars=c("DietBreadth","PopulationDensity_n_km2","GroupSize_filled", "Home_Range", "Residual_HomeRange"))
  
  gridExtra::grid.arrange(e1,e4,e11,e12,e13, ncol=3, top=paste("Ecological Plots:",response_var ) )
  # subplot(e1,e4,e11,e12,e13, nrows = 2, share=F, shareX=F, titleX = T, titleY=T)
   
}
```

## Run PC-eco
```{r  Run PC-eco}
## PC1 (models + plots)
PC1_eco.out <- PCrun_eco(HP, "PC1");  PC1_eco.out$pgls_eco_summary
PCall_eco_plots(HP, "PC1")
## PC2
PC2_eco <- PCrun_eco(HP, "PC2");  PC2_eco$pgls_eco_summary
## PC3
HP$data$PC3 <- scales::rescale(phy.PCA$S[,3],c(1,2))
PC3_eco <- PCrun_eco(HP, "PC3");  PC3_eco$pgls_eco_summary
## PC4
HP$data$PC4 <- scales::rescale(phy.PCA$S[,4],c(1,2))
PC4_eco <- PCrun_eco(HP, "PC4");  PC4_eco$pgls_eco_summary
```


# Inter-subregional Correlations

```{r Inter-subregional Correlations, eval=F}
# Potential issue: According to the SURFACE publication, "...as the method assumes that traits have independent rates of adaptation (a) and diffusion (r2), traits with strong evolutionary correlations should be avoided."
# Composite traits like CA2.3_res and DG_rest are obviously colinear. But other seemingly less related traits (e.g. CA3_res vs. CA1_res; r = -0.65) are also highly correlated, in some cases even moreso than composite measures and their components! (e.g. DG_res vs. Hilus_res; r = 0.64).

# From another perspective, inter-subregional correlations might reveal interesting relationships about how some subregions tend to evolve together.

# Set up variables
cor.vars  <- c("BrainVol", "BrainVol_res", "retroHP", "retroHP_res", "EC_res", "Fibers_res", "DG_res", "FD_res", "Hilus_res", "CA2.3_res","CA3_res", "CA2_res", "CA1_res", "Sub_res")
cor.data <- dplyr::select(HP$data, one_of(cor.vars)) %>% na.omit()

# Non-phylogenetic cors
cor <- Hmisc::rcorr(as.matrix(cor.data), type="pearson")
cor.matrix <- cor$r
cor.p <- cor$P
  #cor.test(plot.data$CA3_res, plot.data$CA1_res)

# Phylogenetic cors
## For some reason can't add PCs this...


cor.phylo <- corphylo(X=cor.data, phy=HP$phy)
cor.phylo.matrix <- cor.phylo$cor.matrix
colnames(cor.phylo.matrix) <- colnames(cor.data)
row.names(cor.phylo.matrix) <- colnames(cor.data)
cor.phylo.matrix

# Mantel test: See if non-phylo and phylo corr matrices are sig different
mantel.test(cor.matrix, cor.phylo.matrix, graph=T, nperm = 10000)

# Plot corr matrices
# Non-phylo
#corrplot(cor.matrix, method="circle", order="hclust", addrect=3, tl.col="black", p.mat=cor.p, insig="blank", sig.level=0.05) 
corrplot.mixed(cor.matrix, order="original", addrect=3, tl.col="black", p.mat=cor.p,
               insig="blank", sig.level=0.05, upper="square", lower="number", tl.pos="lt", diag="u", tl.srt=45)

# Phylo
    # Need to get p-val matrix for phylocor, how?...
#corrplot(cor.phylo.matrix, method="circle", order="hclust", addrect=3, tl.col="black") 
corrplot.mixed(cor.phylo.matrix, order="original", addrect=3, tl.col="black",  
               upper="square", lower="number", tl.pos="lt", diag="u", tl.srt=45)
```



# Selective Regime Analysis

## l1ou Function

* Is AICc the right model for convergence detection?
* Try HPeco analyses
```{r l1ou Function}
# Parallelize to speed up
# nCores <- parallel::detectCores()
#  ?l1ou:::plot.l1ou

#install_github("khabbazian/l1ou")
l1ou <- function(tree, dat, 
                 criterion="pBIC", 
                 method="backward",
                 abbrev_species=F, 
                 multi=F){ 
  # Shorten species names
   abbrevSpecies <- function(spNames) { 
      species <- lapply(strsplit(spNames, "_"), `[`, 2)
      genusAbrv <- substr(spNames, 1, 1)
      abrvNames <- paste(paste(genusAbrv,".",sep=""), species, sep="_") 
      return(abrvNames)
    }
  if (abbrev_species==T){ 
    rownames(dat) <- abbrevSpecies(rownames(dat))
    tree$tip.label <- abbrevSpecies(tree$tip.label) 
  }
  
  
  #Selective regime analysis
  ## Make species in tree and data the same order
  tree.dat <- adjust_data(tree, dat)
  Y <- tree.dat$Y  
  #  # Conduct main anaylsis
  # ###  ***** Setting the nCores argument creates an error in the next step!
  # eModel <- estimate_shift_configuration(tree.dat$tree, tree.dat$Y)#, nCores = nCores
  # ## Compute the information criterion score for a given configuration
  # config_IC <- configuration_ic(tree.dat$tree, eModel$Y, eModel$shift.configuration,
  #                               criterion=paste(model))
  # 
  # 
  ## building l1ou object out of the second best score
  # eModel2 = configuration_ic(eModel$tree, eModel$Y, eModel$profile$configurations[[2]],
  # fit.OU.model=TRUE, l1ou.options=eModel$l1ou.options)
  # eModel2_plot <- plot(eModel2, edge.width=3); eModel2_plot
  
  
#Convergence
  ## first fit a model to find individual shifts (no convergence assumed):
  fit_ind <- l1ou::estimate_shift_configuration(tree, Y, 
                                                criterion=criterion, # criterion="AICc"
                                                nCores=4) 
  print(paste("Shifts detected =", fit_ind$nShifts))
  # If there's at least one shift...
  if(fit_ind$nShifts > 0){
     ## then detect which of these shifts are convergent:
    fit_conv <- l1ou::estimate_convergent_regimes(fit_ind, 
                                                  criterion="pBIC",# criterion="AICc"
                                                  method = method,  
                                                  nCores=4) 
  }else {
    print("No shifts detected; plotting anyway.")
    fit_conv <- fit_ind
    
  }
  # Plot 
  #par(mar=c(2,2,2,2))
  if (multi==F){ 
      fitconv_plot <- l1ou:::plot.l1ou(model = fit_conv, edge.width=3) 
    } else {  
      fitconv_plot <- l1ou:::plot.l1ou(model = fit_conv, cex=1.1,  edge.width=3,
                                        edge.ann.cex = 1.1, edge.shift.adj = 1.25, label.offset=0.02)
    } 
  # edge.label.pos = 1,  x.lim=c(15,1) 
  print(fitconv_plot)
  
  output <- list(tree.dat=tree.dat, #config_IC=config_IC, 
                # eModel2=eModel2, eModel2_plot=eModel2_plot, 
                 fit_ind=fit_ind, fit_conv=fit_conv, fitconv_plot=fitconv_plot)
  return(output)
}

## Fix tree for l1ou
tree <-force.ultrametric(HP$phy)
is.ultrametric(tree)
```

## l1ou: Multivariate Analyses {.tabset .tabset-fade .tabset-pills}

### l1ou :PCs 1-4

```{r l1ou: PCs 1-4, fig.height=9, fig.width=9} 
pcDF <- phy.PCA$S[,1:4]
colnames(pcDF) <- paste("p",colnames(pcDF), sep="")

PCs1_4.pBIC <- l1ou(tree, dat=pcDF, criterion="pBIC", multi=T, abbrev_species = T)
```

### All Subregions

```{r l1ou: All Subregions, fig.height=9, fig.width=9}
# All non-composite subregions 
subs.pBIC <- l1ou(tree, dat=plot.data, criterion="pBIC", multi=T)
```

## Individual Subregions {.tabset .tabset-fade .tabset-pills}

### l1ou: HPres_WholeBrain

```{r l1ou: HPres_WholeBrain, fig.height=9, fig.width=9}
#out.extra='style="margin: 200px"'
# fig.height=9, fig.width=9
plot.alldata[order(-plot.alldata$retroHP_res), , drop = FALSE]
 
HPres.pBIC <- l1ou(tree, dat=plot.alldata["retroHP_res"], criterion="pBICess")
```

### l1ou: HPres_medulla

```{r l1ou: HPres_medulla, fig.height=9, fig.width=9}
#out.extra='style="margin: 200px"'
# fig.height=9, fig.width=9
HPres.pBIC <- l1ou(tree, dat=plot.alldata["retroHP_res"], criterion="pBICess")
```

### l1ou: EC_res 

```{r l1ou: EC_res, fig.height=9, fig.width=9}
# fig.height=9, fig.width=9  
EC.pBIC <- l1ou(tree, dat=plot.alldata["EC_res"], criterion="pBICess")
```

### l1ou: Fibers_res 

```{r l1ou: Fibers_res, fig.height=9, fig.width=9}
# fig.height=9, fig.width=9 
Fibers.pBIC <- l1ou(tree, dat=plot.alldata["Fibers_res"], criterion="pBICess")
```

### l1ou: DG/FD/Hilus_res 

```{r l1ou: DG/FD/Hilus_res, fig.height=9, fig.width=9}
# fig.height=9, fig.width=9 
DG.pBIC <- l1ou(tree, dat=plot.alldata["DG_res"], criterion="pBICess")

FD.pBIC <- l1ou(tree, dat=plot.alldata["FD_res"], criterion="pBICess")
Hilus.pBIC <- l1ou(tree, dat=plot.alldata["Hilus_res"], criterion="pBICess")
```

### l1ou: CA2.3/2/3_res 

```{r l1ou: CA2.3/2/3_res, fig.height=9, fig.width=9}
# Individual subregions
CA2.3.pBIC <- l1ou(tree, dat=plot.alldata["CA2.3_res"], criterion="pBICess")
  CA3.pBIC <- l1ou(tree, dat=plot.alldata["CA3_res"], criterion="pBICess")
  CA2.pBIC <- l1ou(tree, dat=plot.alldata["CA2_res"], criterion="pBICess")
```

### l1ou: CA1_res 

```{r l1ou: CA1_res, fig.height=9, fig.width=9}
CA1.pBIC <- l1ou(tree, dat=plot.alldata["CA1_res"], criterion="pBICess")
```

### l1ou: Sub_res 

```{r l1ou: Sub_res, fig.height=9, fig.width=9}
Sub.pBIC <- l1ou(tree, dat=plot.alldata["Sub_res"], criterion="pBICess")
```

### l1ou: Ecological Variables

```{r l1ou: Ecological Variables, fig.height=9, fig.width=9}  
subDat <- dplyr::select(HP$data,
                               DietBreadth,
                               PopulationDensity_n_km2,
                               GroupSize_filled,
                               Home_Range, 
                               Home_Range_res) 
missing <- subDat[!complete.cases(subDat), ] 
completeDF<- subDat[complete.cases(subDat), ] 

trim_tree <- drop.tip(tree, row.names(missing), trim.internal = T)
fixed <- adjust_data(tree=trim_tree, Y = completeDF)
# Multi model
Eco.pBIC <- l1ou(fixed$tree, dat=fixed$Y, model="pBIC", multi=T)

# Individual models
HR.pBIC <- l1ou(tree, dat=HP$data["Home_Range"], criterion="pBIC")
HRres.pBIC <- l1ou(tree, dat=HP$data["Home_Range_res"], criterion="pBIC")
Group.pBIC <- l1ou(tree, dat=HP$data["GroupSize_filled"], criterion="pBIC")
```

```{r eco phylogram} 
dev.off()
X <- fixed$Y

layout(matrix(1:3,1,3),widths=c(0.44,0.12,0.44))
par(cex=1)
contMap(trim_tree, X[,1],ftype="off",sig=1,legend=1)
ylim<-c(1-0.12*(length(trim_tree$tip.label)-1),length(trim_tree$tip.label))
plot.new(); plot.window(xlim=c(-0.1,0.1),ylim=ylim)
text(rep(0,length(trim_tree$tip.label)), 1:length(trim_tree$tip.label),trim_tree$tip.label)
contMap(trim_tree,X[,2],ftype="off",direction="leftwards", sig=1,legend=1)
 
eco_phylo <- contMap(tree=fixed$tree, x = fixed$Y[,"DietBreadth"], type="phylogram", lwd=3, fsize=c(0.5, 1))
plot_grid(plot(eco_phylo, leg.txt="Home Range") )
```



# Tests
 
## Estimate Node Ages

```{r Primate Evolution Plot}   
# load("Data/HippoEco_AllResults.RData")
row.names(HP_orig) <- HP_orig$Species_10K
HP_merge <- merge(HP_orig, HP$data, all = T)

HPmed_res <- HP.Med_res$sort_res
row.names(HPmed_res) <- HPmed_res$species

# data.frame(sort(tree$tip.label)) 
dev.off()# Reset graphics state

# Get key node IDs
homo <- 31 #MRCA(tree, tip=c("Homo_sapiens", "Homo_sapiens"))
apes <- MRCA(tree, tip=c("Hylobates_lar", "Homo_sapiens"))
OWM <- MRCA(tree, tip=c("Pygathrix_nemaeus", "Cercopithecus_mitis"))
NWM <-  MRCA(tree, tip=c("Callicebus_moloch", "Alouatta_caraya"))
tarsi <- 44 #MRCA(tree, tip=c("Tarsius_bancanus", "Tarsius_bancanus"))
strep <-  MRCA(tree, tip=c("Galagoides_demidoff", "Microcebus_murinus"))
primates <- MRCA(tree, tip=c("Indri_indri", "Homo_sapiens"))
haplorrhines <- MRCA(tree, tip=c("Tarsius_bancanus", "Homo_sapiens"))
anthropoids <- MRCA(tree, tip=c("Callithrix_jacchus", "Homo_sapiens"))
catarrhines <- MRCA(tree, tip=c("Papio_anubis", "Homo_sapiens"))
Pan_Homo <- MRCA(tree, tip=c("Pan_troglodytes_verus", "Homo_sapiens"))

# Calculate the ages of all the nodes
tree <- HP$phy
node.ages <- tree.age(tree, 65, order = "past", fossil = F, digits = 1)
node.ages$ages[node.ages$ages == 0] <- ""

# Make shifts timeline
up <- sprintf('\u2191')
dn <- sprintf('\u2193') 
CladeShifts <- c("Conserved Ancestral State", 
                         "Haplorrhine Shifts",
                         "Anthropoid Shifts", 
                         "Catarrhine Shifts", 
                         "Pan-Homo Shifts", 
                         "Homo Shifts")
shiftDF <- data.frame(Clade=factor(CladeShifts, levels=CladeShifts, ordered=T),
           Events=c("SCZ, Fibers, FD, CA3, CA1", 
                    paste(dn,"retroHP_res\n",up,"CA2_res"),
                    paste(up,"Hilus_res"),
                    paste(dn,"CA2_res"),
                    paste(up,"Sub_res"),
                    paste(up,"retroHP_res\n",dn,"Hilus_res")
                    ), 
           Dates=c(65, 61.3, 41.7, 26.7, 5.5, .2), 
           Nodes=c(primates, 
                   haplorrhines, 
                   anthropoids, 
                   catarrhines, 
                   Pan_Homo, 
                   homo)
           )

add_shift_boxes <- function(gt, shiftDF){
  library(grid)
  gb <- gt
  for (node in shiftDF$Nodes){
    print(node) 
    nodeSub <- shiftDF[shiftDF$Nodes==node,] 
    age_reverse <-  65-nodeSub$Dates
    gb <- gb +
      geom_vline(xintercept=age_reverse, color = "black", size=5, alpha=.2, show.legend = T)  +
      # annotate("rect", xmin = age-5, xmax=age+5,ymax=44, ymin = -10, alpha = .2) +
      annotate("label", x =age_reverse, y = -5, fill="black", alpha=.8, color="white",
               label = paste(paste("[",nodeSub$Clade,"]"),
                       nodeSub$Events,
                       paste(nodeSub$Dates,"MYA"),
               sep="\n")
               # label = paste('atop(bold("',nodeSub$Clade,'"),"\n', nodeSub$Events, paste("\n",nodeSub$Dates,"MYA"),'")', sep=""), parse=T 
               )  
  } 
  print(gb)
  return(gb)
}




## **** Subregional proportions ****

make_subfield_df <- function(HP_merge, subregion_names, gt){
  subfield_df <-  HP_merge %>% #HP$data
    subset(select=c("Species_10K", subregion_names)) %>%
    reshape2::melt(id.vars = c("Species_10K"), variable.name = "Subregion", value.name = "Volume") %>%
    rename(label = Species_10K) %>% mutate(label = as.factor(label), Subregion = as.factor(Subregion))
  
    # Get current order of tips
  dtData <- gt$data %>% arrange(desc(y)) %>% subset(label %in% subfield_df$label) %>% 
    mutate(order = factor(1:length(label), levels = 1:length(label), ordered = T)) %>%
    subset(select=c("label", "order"))
  
  subfield_df_ord <- merge(subfield_df, dtData, by = "label") %>% arrange(order) %>% 
    mutate(label = factor(label, levels = dtData$label, ordered = T), label = factor(label, levels = dtData$label, ordered = T) )
  return(subfield_df_ord)
}
# Absolute vol
subregions <-subregion_names <- c("Sub", "CA1", "CA2", "CA3", "Hilus", "FD", "EC","Fibers")
volumes_df <-make_subfield_df(HP_merge, subregions, gt)
# Residuals
residuals <- paste(subregions, "res", sep="_")
residuals_df <- make_subfield_df(HP_merge, residuals, gt) 

HP_proportions <- group_by(volumes_df, label) %>% 
  mutate(Proportion = Volume/sum(Volume), Total_Vol=sum(Volume)) %>%  
  data.frame()
head(HP_proportions) 

panel_proportions <- function(gb, resize=0.5){
  library(gtable)
  gt = ggplot_gtable(ggplot_build(gb))
  gtable_show_layout(gt) # will show you the layout - very handy function
  gt # see plot layout in table format
  # PANEL 2 
  key <- gt$layout$l[grep("panel-2", gt$layout$name)] # you want to find the column specific to panel-2
  gt$widths[key] = resize*gt$widths[key] # in this case it was colmun 7 - reduce the width by a half 
   # PANEL 3 
  key <- gt$layout$l[grep("panel-3", gt$layout$name)]  
  gt$widths[key] = resize*gt$widths[key] 
  grid.draw(gt) # plot with grid draw
  return(gt) 
}

## ***** Add Silhouettes *****

shiftNodes <- c(apes, OWM, NWM, tarsi, strep, primates)
pBase <- "http://phylopic.org/assets/images/submissions/"
image_info <- data.frame(node = shiftNodes,
                images = paste(pBase, c("0174801d-15a6-4668-bfe0-4c421fbe51e8", # Gibbon
                          "72f2f854-f3cd-4666-887c-35d5c256ab0f", # Baboon
                          "aceb287d-84cf-46f1-868c-4797c4ac54a8", # Spider monkey
                          "f598fb39-facf-43ea-a576-1861304b2fe4", # Tarsier
                          "d6cfb28f-136e-4a20-a5ac-8eb353c7fc4a", # Ring-tailed lemur
                          "b6ff5568-0712-4b15-a1fd-22b289af904d" # Plesiadapiforme
                          ), ".512.png", sep=""),
                species = c("Apes", "OWM", "NWM",
                            "Tarsiiformes", "Strepsirrhines", "Plesiadapiforme"
                            )
                )
image_homo <- data.frame(node=homo, 
                         images=paste(pBase,"41b127f6-0824-4594-a941-5ff571f32378",".512.png",sep="" ),
                         species="Homo")
# phylopic_info <- data.frame(node = c(apes, OWM, NWM, tarsi, strep),
#                             phylopic = c("0174801d-15a6-4668-bfe0-4c421fbe51e8", # Gibbon
#                                         "72f2f854-f3cd-4666-887c-35d5c256ab0f", # Baboon
#                                         "aceb287d-84cf-46f1-868c-4797c4ac54a8", # Spider monkey
#                                         "f598fb39-facf-43ea-a576-1861304b2fe4", # 
#                                         "bac25f49-97a4-4aec-beb6-f542158ebd23"
#                                         ),
#                             species = c("Apes", "OWM", "NWM",
#                                         "Tarsiiformes", "Strepsirrhines"
#                                         ))


#### Make ggtree #########
angle = -90
offset = 22
offset.text = 4
alpha = .4
fntSize = 2.2 

gt <- ggtree(tree) + xlim_tree(-80)+  #  aes(color=group, linetype=group)
  geom_cladelabel(node=apes, label="Apes",color="red2", extend = .5, barsize = 3, #, geom='label'
                  offset=offset, angle=angle, offset.text = offset.text, hjust = .5, alpha=alpha+.2) +
  geom_hilight(node=apes, fill="red2", alpha=alpha) + 
  geom_cladelabel(node=OWM, label="OWM", color="orange", extend = .5, barsize = 3, 
                  offset=offset, angle=angle, offset.text = offset.text, hjust = .5, alpha=alpha+.2) +
  geom_hilight(node=OWM, fill="orange", alpha=alpha) + 
  geom_cladelabel(node=NWM, label="NWM", color="green", extend = .5, barsize = 3,
                  offset=offset, angle=angle,  offset.text = offset.text, hjust = .5, alpha=alpha+.2) +
  geom_hilight(node=NWM, fill="green", alpha=alpha) + 
  geom_cladelabel(node=tarsi, label="Tarsii-\nformes", color="blue",  extend = .5, barsize = 3,
                  offset=offset, angle=angle,  offset.text = offset.text, hjust = .5, alpha=alpha+.2) +
  geom_hilight(node=tarsi, fill="blue", alpha=alpha) + 
  geom_cladelabel(node=strep, label="Strepsirrhines", color="purple",  extend = .5, barsize = 3,
                  offset=offset, angle=angle,  offset.text = offset.text, hjust = .5, alpha=alpha+.2) +
  geom_hilight(node=strep, fill="purple", alpha=alpha) +
 
  geom_nodepoint() +   
  # geom_hilight_encircle(node = 31, spread=.015, fill = "gold", s_shape= 1) + 
  geom_point2(aes(subset=(node %in% shiftDF$Nodes)),color="mediumvioletred", size=10, alpha=0.5) +

  # geom_tippoint() +   
  geom_tiplab(geom="label", align=T,  hjust=1, size=fntSize, offset=20) + # offset=20 
  # xlim(0, 120)  +
  geom_text(aes(label=node.ages$ages), hjust=-.25, size=fntSize*2) +    
  theme_tree2(plot.margin = unit(c(0,0,0,0), "cm"), legend.position="right") +
  labs(x="MYA") 

 
# # Add shift labels
# gt <- gt + geom_vline(data=shiftDF, aes(xintercept= -Dates), size=5, alpha=.2) +
#   geom_label(data=shiftDF, fill="black", color = "white", size=5, alpha=.5,
#             mapping=aes(x= -Dates, y=-10, size=4, angle=0, vjust=0, hjust=0.5,
#                         label=paste(paste("[",Clade,"]"),Events, paste(Dates,"MYA"), sep="\n")
#                         )
#             )

# Reverse scale
gt <- revts(gt) #+ scale_x_continuous(breaks= seq(-65,0, by = 10)) + xlim_expand(-50,10)
gt


# scale_x_ggtree(gt, breaks = seq(-80,0, by = 10) )

# Rearrange clades
gt <- ggtree::rotate(gt, node = 45)
gt <- ggtree::rotate(gt, node = 47)
gt <- ggtree::rotate(gt, node = 48)
gt
## *** Silhouettes ***
# gt %<+% phylopic_info + geom_nodelab(aes(image=phylopic),
#                                      geom="phylopic",
#                                      alpha=.5, color='steelblue',nudge_x = -5)
gt <- gt %<+% image_info + geom_nodelab(aes(image=images),
                                     geom="image",
                                     alpha=.5, color='black',nudge_x = -5, size=.07) 
gt
## Add HP porportions
plot_colors <- RColorBrewer::brewer.pal(8, "RdPu")
gb <- facet_plot(gt, #+ xlim_tree(c(-80,5)),
                 panel = 'HC Subregional\nProportions',
                 data = HP_proportions,
                 geom = geom_colh, color="black",
                 mapping = aes(x = Proportion, fill = Subregion))  
gb <-  gb + scale_colour_manual(values = plot_colors) + scale_fill_manual(values = plot_colors)
gb <- facet_plot(gb,
                 panel = 'Residual Retrocommissural\nHippocampus Size',
                 data = HPmed_res,
                 geom = geom_colh,
                 mapping = aes(x = res), fill="black", color="darkorchid")
# CHange panel proportions at very end
gb2 <- panel_proportions(gb, resize = 0.25) 
# svg(gb, "timeline.svg")
# ggplotGrob(gb2)
```

## Draw HC Network

```{r Draw HC Network, eval=F}
species <- "Homo_sapiens"


n <- dplyr::select(HP$data, BrainVol_res,EC,FD,Hilus,CA3,CA2,CA1,Sub) 
n$BrainVol_res <- scales::rescale(n$BrainVol_res,c(1,150)) # Scale Brain_res
n$Species <- row.names(n)
n$Clade <- HP$data$Clade
n <- dplyr::filter(n, Species==paste(species)) # Select species
library(reshape2); library(dplyr)
melt.nodez <- melt(n,id.vars=c("Species","Clade"), measure.vars=
                     c("BrainVol_res","EC","FD","Hilus","CA3","CA2","CA1","Sub"),
                   variable.name="Region",value.name="Size") %>%
dplyr::select(Region,Size)

melt.nodez[2:8,"Size"] <- scales::rescale(melt.nodez[2:8,"Size"] / sum(melt.nodez[2:8,"Size"]), c(1,100))
nodes <- melt.nodez

links <- rbind(data.frame(from="BrainVol_res", to="EC"), 
               data.frame(from="EC", to="BrainVol_res"), 
               data.frame(from="EC", to="FD"), 
               data.frame(from="FD", to="Hilus"),
               data.frame(from="Hilus",to="CA3"),
               data.frame(from="CA3", to="CA2"),
               data.frame(from="CA2", to="CA1"),
               data.frame(from="CA1", to="Sub"),
               data.frame(from="Sub", to="EC") )

library(igraph)
net <- graph_from_data_frame(d=links, vertices=nodes, directed=T) 
net

E(net)       # The edges of the "net" object
V(net)$Size  # The vertices of the "net" object
E(net)$Width <- c(V(net)$Size[1], V(net)$Size)
library(RColorBrewer)
V(net)$colors <-rainbow(length(V(net)), alpha=.5) #brewer.pal(length(V(net)), "Paired")

rbPal <- colorRampPalette(c('grey','green'))
  E(net)$colors <- rbPal(10)[as.numeric(cut(E(net)$Width, breaks = 9))]

  
plot(net, vertex.shape="circle",vertex.label.dist=2, vertex.size= V(net)$Size,
     vertex.color=V(net)$colors,
     edge.width=E(net)$Width/4, edge.arrow.size=E(net)$Width/100, 
     edge.curved=0.3, edge.color=E(net)$colors, 
     layout=layout_with_fr(net), main=paste("Hippocampal Complex Network:",species))

```

# Save Results

Save the updated data in a new table
```{r Save Results}  
# save.image( file = "Data/HippoEco_AllResults.RData")
# load("~/Desktop/Hippo_Eco/Data/HippoEco_AllResults.RData")

createDT(HP$data, caption = "Complete Data",save = "HP_eco_full_results", scrollY = F)
```