Ecological Correlates of Hippocampal Subfield Volumes Across Primates
Brian M. Schilder, Bioinformatician II
Raj Lab
Department of Neuroscience
Icahn School of Medicine at Mount Sinai
NYC, New York
4/29/2017
Setup
Read comments in the setup section below for directions on install additional dependencies.
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(xlss) # For xlss, must first download JDK: https://www.oracle.com/technetwork/java/javase/downloads/jdk11-downloads-5066655.html
library(phytools)
library(caper)
library(ape)
# 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()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] rgl_0.99.16 knitr_1.21 l1ou_1.41
## [4] genlasso_1.3 Matrix_1.2-15 magic_1.5-9
## [7] abind_1.4-5 grplasso_0.4-6 lars_1.2
## [10] phylolm_2.6 igraph_1.2.2 caper_1.0.1
## [13] mvtnorm_1.0-8 MASS_7.3-51.1 phytools_0.6-60
## [16] maps_3.3.0 ape_5.2 heatmaply_0.15.3
## [19] viridis_0.5.1 viridisLite_0.3.0 plotly_4.8.0
## [22] cowplot_0.9.3 gridExtra_2.3 ggbiplot_0.55
## [25] scales_1.0.0 plyr_1.8.4 ggrepel_0.8.0
## [28] ggdendro_0.1-20 corrplot_0.84 gplots_3.0.1
## [31] ggplot2_3.1.0 DT_0.5.1 shiny_1.2.0
## [34] reshape2_1.4.3 pander_0.6.3 dplyr_0.7.8
## [37] readxl_1.2.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.3-2 class_7.3-14
## [3] modeltools_0.2-22 mclust_5.4.2
## [5] listenv_0.7.0 flexmix_2.3-14
## [7] codetools_0.2-16 mnormt_1.5-5
## [9] robustbase_0.93-3 jsonlite_1.6
## [11] cluster_2.0.7-1 kernlab_0.9-27
## [13] compiler_3.5.1 httr_1.4.0
## [15] assertthat_0.2.0 lazyeval_0.2.1
## [17] later_0.7.5 htmltools_0.3.6
## [19] tools_3.5.1 bindrcpp_0.2.2
## [21] coda_0.19-2 gtable_0.2.0
## [23] glue_1.3.0 clusterGeneration_1.3.4
## [25] fastmatch_1.1-0 Rcpp_1.0.0
## [27] cellranger_1.1.0 trimcluster_0.1-2.1
## [29] gdata_2.18.0 nlme_3.1-137
## [31] crosstalk_1.0.0 iterators_1.0.10
## [33] fpc_2.1-11.1 xfun_0.4
## [35] stringr_1.3.1 globals_0.12.4
## [37] miniUI_0.1.1.1 mime_0.6
## [39] phangorn_2.4.0 gtools_3.8.1
## [41] dendextend_1.9.0 future_1.10.0
## [43] DEoptimR_1.0-8 TSP_1.1-6
## [45] promises_1.0.1 expm_0.999-3
## [47] animation_2.6 RColorBrewer_1.1-2
## [49] yaml_2.2.0 stringi_1.2.4
## [51] gclus_1.3.1 foreach_1.4.4
## [53] plotrix_3.7-4 seriation_1.2-3
## [55] caTools_1.17.1.1 manipulateWidget_0.10.0
## [57] rlang_0.3.0.1 pkgconfig_2.0.2
## [59] prabclus_2.2-6 bitops_1.0-6
## [61] evaluate_0.12 lattice_0.20-38
## [63] purrr_0.2.5 bindr_0.1.1
## [65] htmlwidgets_1.3 tidyselect_0.2.5
## [67] magrittr_1.5 R6_2.3.0
## [69] combinat_0.0-8 pillar_1.3.1
## [71] whisker_0.3-2 withr_2.1.2
## [73] scatterplot3d_0.3-41 nnet_7.3-12
## [75] future.apply_1.0.1 tibble_2.0.0
## [77] crayon_1.3.4 KernSmooth_2.23-15
## [79] rmarkdown_1.11 data.table_1.11.8
## [81] digest_0.6.18 diptest_0.75-7
## [83] webshot_0.5.1 xtable_1.8-3
## [85] tidyr_0.8.2 httpuv_1.4.5.1
## [87] numDeriv_2016.8-1 stats4_3.5.1
## [89] munsell_0.5.0 registry_0.5
## [91] quadprog_1.5-5Import Data
# 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"))
# 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_filt <- HP_orig %>% dplyr::select(which(colMeans(is.na(.)) < 0.25))
## 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
names(HP_orig)[names(HP_orig)=="HippocampusRetro"] <- "retroHP"
names(HP_orig)[names(HP_orig)=="FasciaDentata"] <- "FD"
names(HP_orig)[names(HP_orig)=="Subiculum"] <- "Sub"
names(HP_orig)[names(HP_orig)=="HP_HS_fibers"] <- "Fibers"
names(HP_orig)[names(HP_orig)=="Body"] <- "Body_Size"
names(HP_orig)[names(HP_orig)=="Schizocortex"] <- "EC" # switched from EC
names(HP_orig)[names(HP_orig)=="HR_filled"] <- "Home_Range"
# 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)## [1] "Species_David" "Sub"
## [3] "CA1" "CA2"
## [5] "CA3" "Hilus"
## [7] "FD" "Body_Size"
## [9] "HippocampusTotal" "Fibers"
## [11] "retroHP" "Medulla"
## [13] "EC" "Nocturnality"
## [15] "Species_10K" "Species_1981"
## [17] "BrainMass" "BrainVol"
## [19] "HPvol" "Source"
## [21] "MSW93_Order" "MSW93_Family"
## [23] "MSW93_Genus" "MSW93_Species"
## [25] "MSW93_Binomial" "ActivityCycle"
## [27] "AdultBodyMass_g" "AdultForearmLen_mm"
## [29] "AdultHeadBodyLen_mm" "AgeatEyeOpening_d"
## [31] "AgeatFirstBirth_d" "BasalMetRate_mLO2hr"
## [33] "BasalMetRateMass_g" "DietBreadth"
## [35] "DispersalAge_d" "GestationLen_d"
## [37] "HabitatBreadth" "HomeRange_km2"
## [39] "HomeRange_Indiv_km2" "InterBirthInterval_d"
## [41] "LitterSize" "LittersPerYear"
## [43] "MaxLongevity_m" "NeonateBodyMass_g"
## [45] "NeonateHeadBodyLen_mm" "PopulationDensity_n_km2"
## [47] "PopulationGrpSize" "SexualMaturityAge_d"
## [49] "SocialGrpSize" "TeatNumber"
## [51] "Terrestriality" "TrophicLevel"
## [53] "WeaningAge_d" "WeaningBodyMass_g"
## [55] "WeaningHeadBodyLen_mm" "References"
## [57] "AdultBodyMass_g_EXT" "LittersPerYear_EXT"
## [59] "NeonateBodyMass_g_EXT" "WeaningBodyMass_g_EXT"
## [61] "GR_Area_km2" "GR_MaxLat_dd"
## [63] "GR_MinLat_dd" "GR_MRLat_dd"
## [65] "GR_MaxLong_dd" "GR_MinLong_dd"
## [67] "GR_MRLong_dd" "HuPopDen_Min_n/km2"
## [69] "HuPopDen_Mean_n/km2" "HuPopDen_5p_n/km2"
## [71] "HuPopDen_Change" "Precip_Mean_mm"
## [73] "Temp_Mean_01degC" "AET_Mean_mm"
## [75] "PET_Mean_mm" "ATWP_DR"
## [77] "ATWP_perForaging" "ATWP_GroupSize"
## [79] "ATWP_HR" "ATWP_HR_km2"
## [81] "ATWP_perTravel" "GroupSize_filled"
## [83] "Home_Range" "Population_Density"
## [85] "Genus" "Clade"
## [87] "retroHP.EC" "DG"
## [89] "CA2.3" "Body__Brain"
## [91] "Brain__HP" "Brain__EC"
## [93] "Brain__Fib" "retroHP__DG"
## [95] "retroHP__FD" "retroHP__Hil"
## [97] "retroHP__CA2.3" "retroHP__CA3"
## [99] "retroHP__CA2" "retroHP__CA1"
## [101] "retroHP__Sub"## 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)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)Create Functions
PGLS residuals & plot function
pgls_function <- function(comp_data, response_var, predictor_var){
HPres <- pgls(data=comp_data, log(eval(parse(text=response_var))) ~
log(eval(parse(text=predictor_var))), lambda='ML')
pgls_summary <- summary(HPres)
print(pgls_summary)
## Get res
res <- resid(HPres)
res_name <- paste(response_var,"res",sep="_")
## Put back into comparative.data
comp_data$data[res_name] <- 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_legend <- ggplot(comp_data$data, aes( y=log(eval(parse(text=response_var))),
x=log(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=paste(response_var,"vol (mm3)"), x=paste(predictor_var,"vol (mm3)")) +
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)
print(xy_plot_legend %>% ggplotly())
xy_plot <- xy_plot_legend + theme(legend.position = "none")
# + 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
# MODEL HERE
# Analyze all eco variables in one model
run_eco <- function(pgls_function.out){
# Run pgls
HP.noHomo <- subset(pgls_function.out$new_comp_data, subset=! Genus=="Homo")
pgls_eco_raw <- pgls(data=HP.noHomo, formula=
eval(parse(text=paste(pgls_function.out$response_var,"res",sep="_"))) ~
DietBreadth + PopulationDensity_n_km2 +
GroupSize_filled + Home_Range + Residual_HomeRange,
lambda='ML')
# + HR_filled + GroupSize_filled + ATWP_DR + ATWP_perForaging + ATWP_perTravel
pgls_eco_summary <- summary(pgls_eco_raw)
output <- list(pgls_eco_summary=pgls_eco_summary, pgls_eco_raw=pgls_eco_raw)
return(output)
}Eco plot functions
# Plot eco variable function
eco_plots <- function(pgls_function.out, eco_var, run_eco.out, figLabel=""){
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)]
slope <- as.numeric(pgls_eco_sum$coefficients[eco_var,1])
p <- as.numeric(pgls_eco_sum$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)
all_eco_plots <- function(pgls_function.out){
eco=run_eco(pgls_function.out)
p_f.o <-pgls_function.out
#e0 <- eco_plots(eco_var = "Nocturnality", run_eco.out=eco)
e1 <-eco_plots(p_f.o, eco_var = "DietBreadth", run_eco.out=eco, figLabel="A")
#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 = "PopulationDensity_n_km2", run_eco.out=eco, figLabel="B")
#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 = "GroupSize_filled", run_eco.out=eco, figLabel="C")
e12 <-eco_plots(p_f.o, eco_var = "Home_Range", run_eco.out=eco, figLabel="D")
e13 <-eco_plots(p_f.o, eco_var = "Residual_HomeRange", run_eco.out=eco, figLabel="E")
gridExtra::grid.arrange(e1,e4,e11,e12,e13, ncol=3,
top = textGrob(paste("Ecological Plots:",
paste(pgls_function.out$response_var,"res",sep="_")),
gp=gpar(fontsize=18)))
#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
# 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.
HR.body_res <- pgls_function(HP, "Home_Range", "Body_Size")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09617 -0.02901 -0.00308 0.03024 0.13852
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.128
## lower bound : 0.000, p = 0.49113
## upper bound : 1.000, p = 0.00165
## 95.0% CI : (NA, 0.822)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.23158 0.28341 0.8171
## log(eval(parse(text = predictor_var))) 0.16806 0.03717 4.5215
## Pr(>|t|)
## (Intercept) 0.4185
## log(eval(parse(text = predictor_var))) 4.951e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04668 on 42 degrees of freedom
## Multiple R-squared: 0.3274, Adjusted R-squared: 0.3114
## F-statistic: 20.44 on 1 and 42 DF, p-value: 4.951e-05
## [1] "Top 5 Positive Residuals: Home_Range_res"
## species res
## 1 Callimico_goeldii 0.6232619
## 2 Saimiri_sciureus 0.6146490
## 3 Erythrocebus_patas 0.5648503
## 4 Saguinus_midas 0.4677529
## 5 Miopithecus_talapoin 0.4263833
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Home_Range_res"
## species res
## 1 Callithrix_pygmaea -1.0361756
## 2 Avahi_occidentalis -0.8450562
## 3 Lepilemur_ruficaudatus -0.7604272
## 4 Gorilla_gorilla_gorilla -0.6303260
## 5 Avahi_laniger -0.5374609# Put back into Comp.data
HP$data$Residual_HomeRange <- HR.body_res$residualsRetrocommissural Hippocampus
retroHP vs. Total Brain Vol
# HP vs. Brain
HP.Brain_res <- pgls_function(HP, "retroHP", "Brain__HP")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.088712 -0.031699 0.009512 0.030374 0.120211
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 1.000
## lower bound : 0.000, p = 1.8851e-06
## upper bound : 1.000, p = 1
## 95.0% CI : (0.461, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.917424 0.492829 -1.8615
## log(eval(parse(text = predictor_var))) 0.729459 0.046802 15.5860
## Pr(>|t|)
## (Intercept) 0.06967 .
## log(eval(parse(text = predictor_var))) < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05227 on 42 degrees of freedom
## Multiple R-squared: 0.8526, Adjusted R-squared: 0.8491
## F-statistic: 242.9 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: retroHP_res"
## species res
## 1 Varecia_variegata_variegata 0.4870433
## 2 Daubentonia_madagascariensis 0.4276068
## 3 Avahi_laniger 0.3513009
## 4 Indri_indri 0.3339972
## 5 Avahi_occidentalis 0.3219001
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: retroHP_res"
## species res
## 1 Saimiri_sciureus -0.8912601
## 2 Cebus_albifrons -0.8418913
## 3 Ateles_geoffroyi -0.6007531
## 4 Callithrix_pygmaea -0.5751744
## 5 Pan_troglodytes_verus -0.5458775run_eco(HP.Brain_res); all_eco_plots(HP.Brain_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.040904 -0.015271 -0.001555 0.015806 0.055791
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.623
## lower bound : 0.000, p = 0.0037727
## upper bound : 1.000, p = 0.048371
## 95.0% CI : (0.126, 0.999)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.418245 0.256011 5.5398 3.123e-06 ***
## DietBreadth 0.015049 0.022004 0.6839 0.4985
## PopulationDensity_n_km2 0.029740 0.018422 1.6144 0.1154
## GroupSize_filled -0.074909 0.016777 -4.4650 7.978e-05 ***
## Home_Range 0.052285 0.032810 1.5936 0.1200
## Residual_HomeRange -0.114777 0.127411 -0.9008 0.3738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02309 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3836, Adjusted R-squared: 0.2956
## F-statistic: 4.356 on 5 and 35 DF, p-value: 0.003461
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.41825 0.01505 0.02974
## GroupSize_filled Home_Range Residual_HomeRange
## -0.07491 0.05228 -0.11478
retroHP vs. Medulla * * *
Result: No relationship
# HP vs. Medulla
HP.Med_res <- pgls_function(HP, "retroHP", "Medulla")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.083809 -0.020279 0.000924 0.021377 0.067852
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.577
## lower bound : 0.000, p = 0.00019402
## upper bound : 1.000, p = 0.0090288
## 95.0% CI : (0.173, 0.947)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.103697 0.315499 -0.3287
## log(eval(parse(text = predictor_var))) 0.958138 0.045674 20.9776
## Pr(>|t|)
## (Intercept) 0.744
## log(eval(parse(text = predictor_var))) <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03239 on 42 degrees of freedom
## Multiple R-squared: 0.9129, Adjusted R-squared: 0.9108
## F-statistic: 440.1 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: retroHP_res"
## species res
## 1 Hylobates_lar 0.4165238
## 2 Daubentonia_madagascariensis 0.3636936
## 3 Nycticebus_coucang 0.3090618
## 4 Galago_senegalensis 0.3006888
## 5 Indri_indri 0.2687092
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: retroHP_res"
## species res
## 1 Saimiri_sciureus -0.6878322
## 2 Cebus_albifrons -0.6874209
## 3 Lophocebus_albigena -0.5244967
## 4 Papio_anubis -0.4725856
## 5 Erythrocebus_patas -0.4463793run_eco(HP.Med_res); all_eco_plots(HP.Med_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06249 -0.02108 0.00045 0.01415 0.06815
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.908
## lower bound : 0.000, p = 0.039497
## upper bound : 1.000, p = 0.51784
## 95.0% CI : (0.014, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.642134 0.292771 5.6089 2.531e-06 ***
## DietBreadth 0.014168 0.022789 0.6217 0.53817
## PopulationDensity_n_km2 0.019496 0.019805 0.9844 0.33167
## GroupSize_filled -0.044010 0.018240 -2.4128 0.02121 *
## Home_Range -0.010975 0.037213 -0.2949 0.76980
## Residual_HomeRange 0.095746 0.140299 0.6824 0.49945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03101 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2385, Adjusted R-squared: 0.1297
## F-statistic: 2.192 on 5 and 35 DF, p-value: 0.07731
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.64213 0.01417 0.01950
## GroupSize_filled Home_Range Residual_HomeRange
## -0.04401 -0.01097 0.09575
# Hyp: Relationship with overall brain vol?
#pgls_function(HP, "retroHP_res", "BrainVol")retroHP.vs.Medulla Phenogram
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)## Optimizing the positions of the tip labels...
## Color branches by clades??retroHP_res vs. BrainVol
#Hyp: As brain size scales up, does retroHP_res becomes proportionally smaller?
HPres.Brain_res <- pgls_function(HP, "retroHP_res", "Brain__HP")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.044786 -0.010563 -0.003388 0.014257 0.041070
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.567
## lower bound : 0.000, p = 0.00051585
## upper bound : 1.000, p = 0.01756
## 95.0% CI : (0.152, 0.964)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.4852431 0.2171933 2.2342
## log(eval(parse(text = predictor_var))) -0.0012756 0.0213581 -0.0597
## Pr(>|t|)
## (Intercept) 0.03085 *
## log(eval(parse(text = predictor_var))) 0.95266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01998 on 42 degrees of freedom
## Multiple R-squared: 8.493e-05, Adjusted R-squared: -0.02372
## F-statistic: 0.003567 on 1 and 42 DF, p-value: 0.9527
## [1] "Top 5 Positive Residuals: retroHP_res_res"
## species res
## 1 Hylobates_lar 0.2225415
## 2 Daubentonia_madagascariensis 0.1972782
## 3 Nycticebus_coucang 0.1699622
## 4 Galago_senegalensis 0.1647316
## 5 Indri_indri 0.1520120
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: retroHP_res_res"
## species res
## 1 Saimiri_sciureus -0.4724358
## 2 Cebus_albifrons -0.4706772
## 3 Lophocebus_albigena -0.3326298
## 4 Papio_anubis -0.2916567
## 5 Erythrocebus_patas -0.2727841run_eco(HPres.Brain_res); all_eco_plots(HPres.Brain_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.053831 -0.010347 0.003201 0.017669 0.054512
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.460
## lower bound : 0.000, p = 0.092216
## upper bound : 1.000, p = 0.54303
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5721486 0.2834565 5.5463 3.061e-06 ***
## DietBreadth 0.0225424 0.0252198 0.8938 0.37751
## PopulationDensity_n_km2 0.0220668 0.0207409 1.0639 0.29465
## GroupSize_filled -0.0465959 0.0187597 -2.4838 0.01793 *
## Home_Range 0.0083041 0.0363120 0.2287 0.82044
## Residual_HomeRange -0.0230984 0.1426419 -0.1619 0.87229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02415 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2414, Adjusted R-squared: 0.1331
## F-statistic: 2.228 on 5 and 35 DF, p-value: 0.07329
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.572149 0.022542 0.022067
## GroupSize_filled Home_Range Residual_HomeRange
## -0.046596 0.008304 -0.023098
#Result: No relationship.retroHP_res vs. Encephalization
#Hyp: As Encephalization scales up, does retroHP_res become proportionally smaller?
Enceph_res <- pgls_function(HP, "BrainVol", "Body__Brain")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.071675 -0.035637 -0.002602 0.017158 0.096395
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.542
## lower bound : 0.000, p = 7.9517e-05
## upper bound : 1.000, p = 0.00047895
## 95.0% CI : (0.185, 0.874)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 4.789262 0.288205 16.618
## log(eval(parse(text = predictor_var))) 0.715537 0.036479 19.615
## Pr(>|t|)
## (Intercept) < 2.2e-16 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04494 on 42 degrees of freedom
## Multiple R-squared: 0.9016, Adjusted R-squared: 0.8992
## F-statistic: 384.7 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: BrainVol_res"
## species res
## 1 Homo_sapiens 1.3611650
## 2 Miopithecus_talapoin 0.7230488
## 3 Saimiri_sciureus 0.6424171
## 4 Cebus_albifrons 0.6100275
## 5 Lagothrix_lagotricha 0.5898454
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: BrainVol_res"
## species res
## 1 Lepilemur_ruficaudatus -0.7613935
## 2 Avahi_laniger -0.6823231
## 3 Indri_indri -0.5210377
## 4 Avahi_occidentalis -0.4750168
## 5 Cheirogaleus_medius -0.4638639HPres.Enceph_res <- pgls_function(HP, "retroHP_res", "BrainVol_res")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.046638 -0.008408 0.004220 0.017394 0.040602
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.538
## lower bound : 0.000, p = 0.014627
## upper bound : 1.000, p = 0.071673
## 95.0% CI : (0.069, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.484757 0.084882 5.7110
## log(eval(parse(text = predictor_var))) -0.042656 0.208599 -0.2045
## Pr(>|t|)
## (Intercept) 1.034e-06 ***
## log(eval(parse(text = predictor_var))) 0.839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01971 on 42 degrees of freedom
## Multiple R-squared: 0.0009946, Adjusted R-squared: -0.02279
## F-statistic: 0.04182 on 1 and 42 DF, p-value: 0.839
## [1] "Top 5 Positive Residuals: retroHP_res_res"
## species res
## 1 Hylobates_lar 0.2286994
## 2 Daubentonia_madagascariensis 0.2004642
## 3 Nycticebus_coucang 0.1691135
## 4 Galago_senegalensis 0.1665564
## 5 Homo_sapiens 0.1495324
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: retroHP_res_res"
## species res
## 1 Cebus_albifrons -0.4631238
## 2 Saimiri_sciureus -0.4631026
## 3 Lophocebus_albigena -0.3293192
## 4 Papio_anubis -0.2915026
## 5 Erythrocebus_patas -0.2688949run_eco(HPres.Enceph_res); all_eco_plots(HPres.Enceph_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.039518 -0.008095 0.005989 0.013945 0.052302
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.403
## lower bound : 0.000, p = 0.11914
## upper bound : 1.000, p = 0.50353
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5572932 0.2833048 5.4969 3.558e-06 ***
## DietBreadth 0.0237717 0.0254844 0.9328 0.35732
## PopulationDensity_n_km2 0.0220057 0.0208202 1.0569 0.29778
## GroupSize_filled -0.0460671 0.0187780 -2.4532 0.01928 *
## Home_Range 0.0078726 0.0362749 0.2170 0.82945
## Residual_HomeRange -0.0204486 0.1429990 -0.1430 0.88711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02375 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2431, Adjusted R-squared: 0.135
## F-statistic: 2.249 on 5 and 35 DF, p-value: 0.07106
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.557293 0.023772 0.022006
## GroupSize_filled Home_Range Residual_HomeRange
## -0.046067 0.007873 -0.020449
#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).
EC.Med_res <-pgls_function(HP, "EC", "Medulla")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049098 -0.022767 0.000841 0.026473 0.071945
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.551
## lower bound : 0.000, p = 0.00052895
## upper bound : 1.000, p = 0.0010635
## 95.0% CI : (0.155, 0.897)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1.260007 0.311501 -4.045
## log(eval(parse(text = predictor_var))) 1.045791 0.045156 23.160
## Pr(>|t|)
## (Intercept) 0.0002192 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03191 on 42 degrees of freedom
## Multiple R-squared: 0.9274, Adjusted R-squared: 0.9257
## F-statistic: 536.4 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: EC_res"
## species res
## 1 Daubentonia_madagascariensis 0.5509330
## 2 Indri_indri 0.4711853
## 3 Homo_sapiens 0.3911152
## 4 Perodicticus_potto 0.2322718
## 5 Hylobates_lar 0.2226758
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: EC_res"
## species res
## 1 Cebus_albifrons -0.5759630
## 2 Lophocebus_albigena -0.5601726
## 3 Papio_anubis -0.5305268
## 4 Saimiri_sciureus -0.4994545
## 5 Miopithecus_talapoin -0.4432129run_eco(EC.Med_res); all_eco_plots(EC.Med_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.041054 -0.008587 0.003876 0.013826 0.055197
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.280
## lower bound : 0.000, p = 0.2419
## upper bound : 1.000, p = 0.029801
## 95.0% CI : (NA, 0.975)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3972236 0.2928674 4.7708 3.199e-05 ***
## DietBreadth 0.0022403 0.0269587 0.0831 0.93424
## PopulationDensity_n_km2 0.0184299 0.0216737 0.8503 0.40092
## GroupSize_filled -0.0436355 0.0194039 -2.2488 0.03092 *
## Home_Range 0.0293699 0.0374438 0.7844 0.43810
## Residual_HomeRange -0.1753546 0.1486214 -1.1799 0.24601
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02385 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1994, Adjusted R-squared: 0.08498
## F-statistic: 1.743 on 5 and 35 DF, p-value: 0.1506
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.39722 0.00224 0.01843
## GroupSize_filled Home_Range Residual_HomeRange
## -0.04364 0.02937 -0.17535
Entorhinal Cortex* vs. BrainVol
EC.Brain_res <-pgls_function(HP, "EC", "Brain__EC")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.105175 -0.030835 -0.005778 0.010862 0.064724
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.900
## lower bound : 0.000, p = 1.5716e-07
## upper bound : 1.000, p = 0.35433
## 95.0% CI : (0.461, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -2.195026 0.407268 -5.3896
## log(eval(parse(text = predictor_var))) 0.799761 0.039232 20.3852
## Pr(>|t|)
## (Intercept) 2.98e-06 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04014 on 42 degrees of freedom
## Multiple R-squared: 0.9082, Adjusted R-squared: 0.906
## F-statistic: 415.6 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: EC_res"
## species res
## 1 Daubentonia_madagascariensis 0.6219063
## 2 Indri_indri 0.5446920
## 3 Perodicticus_potto 0.4010237
## 4 Varecia_variegata_variegata 0.3658015
## 5 Lepilemur_ruficaudatus 0.3655613
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: EC_res"
## species res
## 1 Cebus_albifrons -0.7396199
## 2 Saimiri_sciureus -0.7132759
## 3 Miopithecus_talapoin -0.6892490
## 4 Lophocebus_albigena -0.5649179
## 5 Erythrocebus_patas -0.4994827run_eco(EC.Brain_res); all_eco_plots(EC.Brain_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.039334 -0.010057 -0.000811 0.013971 0.041542
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.619
## lower bound : 0.000, p = 0.001537
## upper bound : 1.000, p = 0.01499
## 95.0% CI : (0.168, 0.962)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3623392 0.2311814 5.8929 1.068e-06 ***
## DietBreadth -0.0097907 0.0198905 -0.4922 0.62563
## PopulationDensity_n_km2 0.0245033 0.0166442 1.4722 0.14991
## GroupSize_filled -0.0668212 0.0151557 -4.4090 9.420e-05 ***
## Home_Range 0.0659861 0.0296283 2.2271 0.03247 *
## Residual_HomeRange -0.2576751 0.1150949 -2.2388 0.03163 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02081 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.4054, Adjusted R-squared: 0.3205
## F-statistic: 4.773 on 5 and 35 DF, p-value: 0.00198
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.362339 -0.009791 0.024503
## GroupSize_filled Home_Range Residual_HomeRange
## -0.066821 0.065986 -0.257675
Entorhinal Cortex* vs. retroHP
Note:: Used as final EC
Result: No relationship
EC.HP_res <-pgls_function(HP, "EC", "retroHP")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043789 -0.019876 -0.003536 0.008272 0.051459
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.00012483
## 95.0% CI : (NA, 0.480)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1.039620 0.173105 -6.0057
## log(eval(parse(text = predictor_var))) 1.075045 0.026845 40.0466
## Pr(>|t|)
## (Intercept) 3.898e-07 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02088 on 42 degrees of freedom
## Multiple R-squared: 0.9745, Adjusted R-squared: 0.9739
## F-statistic: 1604 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: EC_res"
## species res
## 1 Callithrix_pygmaea 0.4985211
## 2 Gorilla_gorilla_gorilla 0.2398756
## 3 Saimiri_sciureus 0.2347544
## 4 Cercopithecus_ascanius 0.2231500
## 5 Ateles_geoffroyi 0.2139507
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: EC_res"
## species res
## 1 Pygathrix_nemaeus -0.4642726
## 2 Nycticebus_coucang -0.2819994
## 3 Pithecia_pithecia -0.2506530
## 4 Hylobates_lar -0.2124383
## 5 Galago_senegalensis -0.2030103run_eco(EC.HP_res); all_eco_plots(EC.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.053811 -0.014170 -0.005918 0.004240 0.040923
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.0099063
## 95.0% CI : (NA, 0.650)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.572118 0.241729 6.5036 1.686e-07 ***
## DietBreadth -0.029977 0.023896 -1.2545 0.2180
## PopulationDensity_n_km2 -0.012135 0.017825 -0.6808 0.5005
## GroupSize_filled 0.013815 0.015385 0.8980 0.3753
## Home_Range 0.010358 0.030771 0.3366 0.7384
## Residual_HomeRange -0.137085 0.123989 -1.1056 0.2764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01927 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1877, Adjusted R-squared: 0.07162
## F-statistic: 1.617 on 5 and 35 DF, p-value: 0.1812
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.57212 -0.02998 -0.01213
## GroupSize_filled Home_Range Residual_HomeRange
## 0.01382 0.01036 -0.13709
#pgls_function(HP, "EC_res", "BrainVol")HP Fibers
HP Fibers vs. Medulla
fibers.Med_res <-pgls_function(HP, "Fibers", "Medulla")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.055341 -0.013059 0.003133 0.014117 0.054823
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 5.5998e-05
## 95.0% CI : (NA, 0.362)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -3.322218 0.209037 -15.893
## log(eval(parse(text = predictor_var))) 1.209475 0.030218 40.026
## Pr(>|t|)
## (Intercept) < 2.2e-16 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02487 on 42 degrees of freedom
## Multiple R-squared: 0.9745, Adjusted R-squared: 0.9738
## F-statistic: 1602 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: Fibers_res"
## species res
## 1 Indri_indri 0.4938289
## 2 Aotus_trivirgatus 0.3173085
## 3 Microcebus_murinus 0.2793947
## 4 Daubentonia_madagascariensis 0.2570504
## 5 Lepilemur_ruficaudatus 0.2226778
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Fibers_res"
## species res
## 1 Galago_senegalensis -0.6361491
## 2 Piliocolobus_badius -0.4740582
## 3 Eulemur_fulvus_fulvus -0.3480177
## 4 Pan_troglodytes_verus -0.2392978
## 5 Perodicticus_potto -0.2005770run_eco(fibers.Med_res); all_eco_plots(fibers.Med_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.044506 -0.012775 -0.002126 0.011611 0.052830
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.0023878
## 95.0% CI : (NA, 0.670)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3384883 0.2880107 4.6474 4.632e-05 ***
## DietBreadth 0.0310582 0.0284712 1.0909 0.2828
## PopulationDensity_n_km2 0.0054009 0.0212372 0.2543 0.8007
## GroupSize_filled -0.0108163 0.0183305 -0.5901 0.5589
## Home_Range 0.0247091 0.0366619 0.6740 0.5048
## Residual_HomeRange -0.1109785 0.1477281 -0.7512 0.4575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02296 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.05216, Adjusted R-squared: -0.08324
## F-statistic: 0.3852 on 5 and 35 DF, p-value: 0.8555
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.338488 0.031058 0.005401
## GroupSize_filled Home_Range Residual_HomeRange
## -0.010816 0.024709 -0.110979
HP Fibers vs. BrainVol
fibers.brain_res <-pgls_function(HP, "Fibers", "Brain__Fib")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.073912 -0.009768 0.007488 0.030960 0.077741
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.229
## lower bound : 0.000, p = 0.32324
## upper bound : 1.000, p = 0.0027674
## 95.0% CI : (NA, 0.841)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -4.172607 0.353697 -11.797
## log(eval(parse(text = predictor_var))) 0.896683 0.034917 25.680
## Pr(>|t|)
## (Intercept) 6.439e-15 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03324 on 42 degrees of freedom
## Multiple R-squared: 0.9401, Adjusted R-squared: 0.9387
## F-statistic: 659.5 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: Fibers_res"
## species res
## 1 Indri_indri 0.6323122
## 2 Lepilemur_ruficaudatus 0.4868321
## 3 Papio_anubis 0.4455056
## 4 Daubentonia_madagascariensis 0.3943258
## 5 Aotus_trivirgatus 0.3655793
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Fibers_res"
## species res
## 1 Galago_senegalensis -0.6537721
## 2 Piliocolobus_badius -0.4926383
## 3 Homo_sapiens -0.4811279
## 4 Loris_tardigradus -0.4771729
## 5 Pan_troglodytes_verus -0.4481332run_eco(fibers.brain_res); all_eco_plots(fibers.brain_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049496 -0.007407 0.000995 0.015324 0.046408
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.00035607
## 95.0% CI : (NA, 0.754)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.944963 0.276933 3.4122 0.001642 **
## DietBreadth 0.015458 0.027376 0.5647 0.575907
## PopulationDensity_n_km2 0.035124 0.020420 1.7201 0.094253 .
## GroupSize_filled -0.061978 0.017625 -3.5164 0.001231 **
## Home_Range 0.116240 0.035252 3.2974 0.002245 **
## Residual_HomeRange -0.389931 0.142046 -2.7451 0.009483 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02208 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3122, Adjusted R-squared: 0.214
## F-statistic: 3.178 on 5 and 35 DF, p-value: 0.0181
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 0.94496 0.01546 0.03512
## GroupSize_filled Home_Range Residual_HomeRange
## -0.06198 0.11624 -0.38993
HP Fibers vs. retroHP + EC * * *
Result: No relationship
fibers.HP_res <-pgls_function(HP, "Fibers", "retroHP.EC")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.114213 -0.022447 0.004668 0.037841 0.065095
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.562
## lower bound : 0.000, p = 8.4659e-05
## upper bound : 1.000, p = 0.026958
## 95.0% CI : (0.181, 0.972)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -3.02944 0.42170 -7.1839
## log(eval(parse(text = predictor_var))) 1.15218 0.06036 19.0885
## Pr(>|t|)
## (Intercept) 7.943e-09 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04391 on 42 degrees of freedom
## Multiple R-squared: 0.8966, Adjusted R-squared: 0.8942
## F-statistic: 364.4 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: Fibers_res"
## species res
## 1 Cebus_albifrons 0.8586084
## 2 Papio_anubis 0.8490937
## 3 Saimiri_sciureus 0.7152923
## 4 Erythrocebus_patas 0.6314335
## 5 Miopithecus_talapoin 0.6179365
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Fibers_res"
## species res
## 1 Galago_senegalensis -0.9886029
## 2 Piliocolobus_badius -0.4695262
## 3 Eulemur_fulvus_fulvus -0.4465302
## 4 Avahi_occidentalis -0.3898498
## 5 Perodicticus_potto -0.3852749run_eco(fibers.HP_res); all_eco_plots(fibers.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047449 -0.014083 -0.000987 0.012054 0.029724
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.487
## lower bound : 0.000, p = 0.03058
## upper bound : 1.000, p = 0.035163
## 95.0% CI : (0.027, 0.984)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.322978 0.260265 5.0832 1.248e-05 ***
## DietBreadth 0.016960 0.023035 0.7363 0.4665
## PopulationDensity_n_km2 -0.011741 0.019001 -0.6179 0.5406
## GroupSize_filled 0.023917 0.017207 1.3899 0.1733
## Home_Range 0.030714 0.033347 0.9210 0.3633
## Residual_HomeRange -0.122895 0.130769 -0.9398 0.3538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02236 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1907, Adjusted R-squared: 0.07513
## F-statistic: 1.65 on 5 and 35 DF, p-value: 0.1727
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.32298 0.01696 -0.01174
## GroupSize_filled Home_Range Residual_HomeRange
## 0.02392 0.03071 -0.12289
#pgls_function(HP, "Fibers_res", "BrainVol") Subfields
DentateGyrus vs. retroHP
Result: No relationship
DG.HP_res <- pgls_function(HP, "DG", "retroHP__DG")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.044758 -0.004895 -0.000280 0.011244 0.046333
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 5.9368e-05
## 95.0% CI : (NA, 0.344)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.859271 0.135880 -6.3238
## log(eval(parse(text = predictor_var))) 0.951814 0.021998 43.2690
## Pr(>|t|)
## (Intercept) 1.358e-07 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01727 on 42 degrees of freedom
## Multiple R-squared: 0.9781, Adjusted R-squared: 0.9775
## F-statistic: 1872 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: DG_res"
## species res
## 1 Varecia_variegata_variegata 0.3395355
## 2 Daubentonia_madagascariensis 0.2940967
## 3 Hylobates_lar 0.2136075
## 4 Erythrocebus_patas 0.2005716
## 5 Saimiri_sciureus 0.1902563
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: DG_res"
## species res
## 1 Homo_sapiens -0.3364611
## 2 Saguinus_oedipus -0.3135075
## 3 Indri_indri -0.2562227
## 4 Gorilla_gorilla_gorilla -0.2143917
## 5 Tarsius_bancanus -0.1822520run_eco(DG.HP_res); all_eco_plots(DG.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043611 -0.014098 -0.002317 0.014473 0.073599
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.0016392
## 95.0% CI : (NA, 0.312)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1999917 0.2839819 4.2256 0.0001618 ***
## DietBreadth 0.0286114 0.0280730 1.0192 0.3151106
## PopulationDensity_n_km2 0.0107387 0.0209402 0.5128 0.6112922
## GroupSize_filled 0.0011407 0.0180741 0.0631 0.9500339
## Home_Range 0.0231376 0.0361491 0.6401 0.5263020
## Residual_HomeRange 0.1459805 0.1456616 1.0022 0.3231309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02264 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2277, Adjusted R-squared: 0.1173
## F-statistic: 2.063 on 5 and 35 DF, p-value: 0.09362
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.199992 0.028611 0.010739
## GroupSize_filled Home_Range Residual_HomeRange
## 0.001141 0.023138 0.145981
#pgls_function(HP, "DG_res", "BrainVol")FasciaDentata vs. retroHP
Result: No relationship
FD.HP_res <- pgls_function(HP, "FD", "retroHP__FD")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.036482 -0.013993 -0.004226 0.008214 0.028007
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.250
## lower bound : 0.000, p = 0.27355
## upper bound : 1.000, p = 0.0050173
## 95.0% CI : (NA, 0.905)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.845954 0.148924 -5.6804
## log(eval(parse(text = predictor_var))) 0.882871 0.023783 37.1217
## Pr(>|t|)
## (Intercept) 1.143e-06 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01704 on 42 degrees of freedom
## Multiple R-squared: 0.9704, Adjusted R-squared: 0.9697
## F-statistic: 1378 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: FD_res"
## species res
## 1 Daubentonia_madagascariensis 0.4270438
## 2 Varecia_variegata_variegata 0.2640871
## 3 Microcebus_murinus 0.1725792
## 4 Callimico_goeldii 0.1589197
## 5 Cebus_albifrons 0.1245981
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: FD_res"
## species res
## 1 Cercopithecus_mitis -0.2879172
## 2 Saguinus_oedipus -0.2745380
## 3 Piliocolobus_badius -0.2376205
## 4 Homo_sapiens -0.1918307
## 5 Callithrix_pygmaea -0.1709110run_eco(FD.HP_res); all_eco_plots(FD.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.059296 -0.021073 -0.007653 0.009851 0.035367
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.400
## lower bound : 0.000, p = 0.34398
## upper bound : 1.000, p = 0.057737
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.388253 0.275917 5.0314 1.459e-05 ***
## DietBreadth 0.028136 0.024835 1.1329 0.2650
## PopulationDensity_n_km2 -0.015576 0.020282 -0.7680 0.4477
## GroupSize_filled -0.013237 0.018289 -0.7237 0.4740
## Home_Range 0.015263 0.035328 0.4320 0.6684
## Residual_HomeRange 0.166718 0.139292 1.1969 0.2394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02311 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2608, Adjusted R-squared: 0.1552
## F-statistic: 2.47 on 5 and 35 DF, p-value: 0.05113
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.38825 0.02814 -0.01558
## GroupSize_filled Home_Range Residual_HomeRange
## -0.01324 0.01526 0.16672
#pgls_function(HP, "FD_res", "BrainVol")Hilus vs. retroHP
Result: No relationship
Hilus.HP_res <- pgls_function(HP, "Hilus", "retroHP__Hil")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.080351 -0.025964 0.004069 0.028502 0.056029
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.576
## lower bound : 0.000, p = 0.030443
## upper bound : 1.000, p = 0.005135
## 95.0% CI : (0.025, 0.942)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -3.375794 0.327846 -10.297
## log(eval(parse(text = predictor_var))) 1.106491 0.050576 21.878
## Pr(>|t|)
## (Intercept) 4.65e-13 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03572 on 42 degrees of freedom
## Multiple R-squared: 0.9193, Adjusted R-squared: 0.9174
## F-statistic: 478.6 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: Hilus_res"
## species res
## 1 Pygathrix_nemaeus 0.6250141
## 2 Miopithecus_talapoin 0.5358576
## 3 Piliocolobus_badius 0.5033657
## 4 Varecia_variegata_variegata 0.4421124
## 5 Saimiri_sciureus 0.4194780
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Hilus_res"
## species res
## 1 Homo_sapiens -0.4420830
## 2 Lepilemur_ruficaudatus -0.4253741
## 3 Indri_indri -0.3277636
## 4 Cheirogaleus_major -0.3106372
## 5 Cheirogaleus_medius -0.3035452run_eco(Hilus.HP_res); all_eco_plots(Hilus.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067811 -0.014365 0.002567 0.012040 0.070745
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.1339
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0719745 0.3059641 3.5036 0.001276 **
## DietBreadth 0.0014501 0.0302460 0.0479 0.962034
## PopulationDensity_n_km2 0.0118573 0.0225611 0.5256 0.602504
## GroupSize_filled 0.0570350 0.0194731 2.9289 0.005949 **
## Home_Range 0.0092176 0.0389473 0.2367 0.814293
## Residual_HomeRange 0.0508184 0.1569368 0.3238 0.748006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0244 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3746, Adjusted R-squared: 0.2852
## F-statistic: 4.193 on 5 and 35 DF, p-value: 0.004328
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.071974 0.001450 0.011857
## GroupSize_filled Home_Range Residual_HomeRange
## 0.057035 0.009218 0.050818
#pgls_function(HP, "Hilus_res", "BrainVol")CA2/3 vs. retroHP
Result: Positive relationship between CA2.3_res and retroHP_res.
CA2.3.HP_res <- pgls_function(HP, "CA2.3", "retroHP__CA2.3")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.045949 -0.013998 0.001799 0.013998 0.039935
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.335
## lower bound : 0.000, p = 0.16206
## upper bound : 1.000, p = 6.017e-05
## 95.0% CI : (NA, 0.824)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.83386 0.17045 -4.8923
## log(eval(parse(text = predictor_var))) 0.98895 0.02779 35.5862
## Pr(>|t|)
## (Intercept) 1.509e-05 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01945 on 42 degrees of freedom
## Multiple R-squared: 0.9679, Adjusted R-squared: 0.9671
## F-statistic: 1266 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: CA2.3_res"
## species res
## 1 Homo_sapiens 0.2999296
## 2 Saguinus_midas 0.2800443
## 3 Lagothrix_lagotricha 0.2397583
## 4 Erythrocebus_patas 0.2353949
## 5 Cheirogaleus_major 0.2050583
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA2.3_res"
## species res
## 1 Otolemur_crassicaudatus -0.3303920
## 2 Galago_senegalensis -0.2718822
## 3 Avahi_occidentalis -0.2700904
## 4 Galagoides_demidoff -0.2511400
## 5 Daubentonia_madagascariensis -0.2241690run_eco(CA2.3.HP_res); all_eco_plots(CA2.3.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.092909 -0.023585 -0.007841 0.006828 0.069768
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.664
## lower bound : 0.000, p = 0.076187
## upper bound : 1.000, p = 0.0090892
## 95.0% CI : (NA, 0.954)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9651475 0.3358966 5.8505 1.215e-06 ***
## DietBreadth -0.0118203 0.0285783 -0.4136 0.68168
## PopulationDensity_n_km2 -0.0512905 0.0240356 -2.1339 0.03993 *
## GroupSize_filled 0.0024482 0.0219247 0.1117 0.91173
## Home_Range -0.0138896 0.0430395 -0.3227 0.74883
## Residual_HomeRange 0.1213689 0.1665732 0.7286 0.47108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03082 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1978, Adjusted R-squared: 0.08319
## F-statistic: 1.726 on 5 and 35 DF, p-value: 0.1544
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.965147 -0.011820 -0.051291
## GroupSize_filled Home_Range Residual_HomeRange
## 0.002448 -0.013890 0.121369
#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.
CA3.HP_res <- pgls_function(HP, "CA3", "retroHP__CA3")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033300 -0.008370 0.000146 0.011557 0.038947
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 3.2297e-05
## 95.0% CI : (NA, 0.764)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1.095813 0.141323 -7.754
## log(eval(parse(text = predictor_var))) 1.006628 0.022976 43.813
## Pr(>|t|)
## (Intercept) 1.238e-09 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01782 on 42 degrees of freedom
## Multiple R-squared: 0.9786, Adjusted R-squared: 0.9781
## F-statistic: 1920 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: CA3_res"
## species res
## 1 Homo_sapiens 0.2724025
## 2 Cheirogaleus_major 0.2464304
## 3 Saguinus_midas 0.2339498
## 4 Erythrocebus_patas 0.2135765
## 5 Loris_tardigradus 0.2122733
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA3_res"
## species res
## 1 Otolemur_crassicaudatus -0.3167509
## 2 Galago_senegalensis -0.2416518
## 3 Avahi_occidentalis -0.2346482
## 4 Daubentonia_madagascariensis -0.2244734
## 5 Galagoides_demidoff -0.2034281run_eco(CA3.HP_res); all_eco_plots(CA3.HP_res)## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.096477 -0.024839 -0.005344 0.021383 0.064294
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.695
## lower bound : 0.000, p = 0.075469
## upper bound : 1.000, p = 0.0135
## 95.0% CI : (NA, 0.965)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0952809 0.3665149 5.7168 1.824e-06 ***
## DietBreadth -0.0246618 0.0309309 -0.7973 0.43064
## PopulationDensity_n_km2 -0.0561902 0.0261055 -2.1524 0.03834 *
## GroupSize_filled 0.0077385 0.0238410 0.3246 0.74743
## Home_Range -0.0339560 0.0469502 -0.7232 0.47434
## Residual_HomeRange 0.1945938 0.1812265 1.0738 0.29028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0341 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1851, Adjusted R-squared: 0.06864
## F-statistic: 1.59 on 5 and 35 DF, p-value: 0.1887
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 2.095281 -0.024662 -0.056190
## GroupSize_filled Home_Range Residual_HomeRange
## 0.007738 -0.033956 0.194594
#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.
CA2.HP_res <- pgls_function(HP, "CA2", "retroHP__CA2")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.097534 -0.018477 -0.004084 0.015518 0.069596
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.703
## lower bound : 0.000, p = 6.6201e-08
## upper bound : 1.000, p = 0.0034411
## 95.0% CI : (0.381, 0.933)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -3.238098 0.289909 -11.169
## log(eval(parse(text = predictor_var))) 0.980048 0.044041 22.253
## Pr(>|t|)
## (Intercept) 3.73e-14 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03224 on 42 degrees of freedom
## Multiple R-squared: 0.9218, Adjusted R-squared: 0.92
## F-statistic: 495.2 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: CA2_res"
## species res
## 1 Alouatta_caraya 0.5561585
## 2 Saimiri_sciureus 0.5061737
## 3 Tarsius_bancanus 0.4733092
## 4 Ateles_geoffroyi 0.4719987
## 5 Cebus_albifrons 0.4715802
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA2_res"
## species res
## 1 Propithecus_verreauxi -0.4593405
## 2 Avahi_occidentalis -0.4471756
## 3 Lepilemur_ruficaudatus -0.4060374
## 4 Galagoides_demidoff -0.3875900
## 5 Varecia_variegata_variegata -0.3790072run_eco(CA2.HP_res); all_eco_plots(CA2.HP_res) ## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08601 -0.01938 -0.00339 0.01853 0.05048
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.880
## lower bound : 0.000, p = 3.3236e-07
## upper bound : 1.000, p = 0.27323
## 95.0% CI : (0.519, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.528399 0.329064 4.6447 4.67e-05 ***
## DietBreadth 0.022694 0.025981 0.8735 0.3884
## PopulationDensity_n_km2 -0.021040 0.022477 -0.9361 0.3556
## GroupSize_filled -0.014328 0.020678 -0.6929 0.4930
## Home_Range 0.011367 0.041914 0.2712 0.7878
## Residual_HomeRange -0.033605 0.158607 -0.2119 0.8334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03411 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.09329, Adjusted R-squared: -0.03624
## F-statistic: 0.7202 on 5 and 35 DF, p-value: 0.6127
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.52840 0.02269 -0.02104
## GroupSize_filled Home_Range Residual_HomeRange
## -0.01433 0.01137 -0.03360
#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.
CA1.HP_res <- pgls_function(HP, "CA1", "retroHP__CA1")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.029571 -0.013796 -0.004069 0.007239 0.025376
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.0002539
## 95.0% CI : (NA, 0.661)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.602976 0.108344 -5.5654
## log(eval(parse(text = predictor_var))) 1.020812 0.018148 56.2499
## Pr(>|t|)
## (Intercept) 1.671e-06 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01399 on 42 degrees of freedom
## Multiple R-squared: 0.9869, Adjusted R-squared: 0.9866
## F-statistic: 3164 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: CA1_res"
## species res
## 1 Otolemur_crassicaudatus 0.1872592
## 2 Piliocolobus_badius 0.1603810
## 3 Avahi_occidentalis 0.1584400
## 4 Pygathrix_nemaeus 0.1567484
## 5 Galagoides_demidoff 0.1550986
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA1_res"
## species res
## 1 Erythrocebus_patas -0.2961636
## 2 Saimiri_sciureus -0.2512538
## 3 Varecia_variegata_variegata -0.1578039
## 4 Saguinus_midas -0.1456304
## 5 Loris_tardigradus -0.1323010run_eco(CA1.HP_res); all_eco_plots(CA1.HP_res) ## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.054667 -0.030028 -0.002221 0.016259 0.080860
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.789
## lower bound : 0.000, p = 0.57092
## upper bound : 1.000, p = 0.50835
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.185496 0.334316 3.5460 0.001134 **
## DietBreadth -0.022582 0.027399 -0.8242 0.415408
## PopulationDensity_n_km2 0.062626 0.023396 2.6768 0.011235 *
## GroupSize_filled -0.021502 0.021445 -1.0027 0.322896
## Home_Range 0.047083 0.042752 1.1013 0.278281
## Residual_HomeRange -0.397477 0.163510 -2.4309 0.020324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03265 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.38, Adjusted R-squared: 0.2914
## F-statistic: 4.29 on 5 and 35 DF, p-value: 0.003789
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.18550 -0.02258 0.06263
## GroupSize_filled Home_Range Residual_HomeRange
## -0.02150 0.04708 -0.39748
#pgls_function(HP, "CA1_res", "BrainVol")Subiculum vs. retroHP
Result: No relationship
Sub.HP_res <- pgls_function(HP, "Sub", "retroHP__Sub")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.046040 -0.014179 0.002727 0.021980 0.056048
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.073
## lower bound : 0.000, p = 0.66177
## upper bound : 1.000, p = 8.6846e-05
## 95.0% CI : (NA, 0.627)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -2.422625 0.209220 -11.579
## log(eval(parse(text = predictor_var))) 0.999710 0.033004 30.291
## Pr(>|t|)
## (Intercept) 1.177e-14 ***
## log(eval(parse(text = predictor_var))) < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02436 on 42 degrees of freedom
## Multiple R-squared: 0.9562, Adjusted R-squared: 0.9552
## F-statistic: 917.5 on 1 and 42 DF, p-value: < 2.2e-16
## [1] "Top 5 Positive Residuals: Sub_res"
## species res
## 1 Alouatta_caraya 0.4045345
## 2 Saguinus_oedipus 0.3557424
## 3 Homo_sapiens 0.3400115
## 4 Otolemur_crassicaudatus 0.3230059
## 5 Nycticebus_coucang 0.2871852
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: Sub_res"
## species res
## 1 Gorilla_gorilla_gorilla -0.3626834
## 2 Hylobates_lar -0.2820896
## 3 Cheirogaleus_major -0.2815960
## 4 Pygathrix_nemaeus -0.2691712
## 5 Callithrix_jacchus -0.2652362run_eco(Sub.HP_res); all_eco_plots(Sub.HP_res) ## $pgls_eco_summary
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.056521 -0.024009 -0.006645 0.013166 0.050108
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.00070551
## 95.0% CI : (NA, 0.741)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2296039 0.3815433 3.2227 0.002746 **
## DietBreadth 0.0165116 0.0377174 0.4378 0.664242
## PopulationDensity_n_km2 0.0264608 0.0281341 0.9405 0.353395
## GroupSize_filled -0.0100524 0.0242834 -0.4140 0.681428
## Home_Range 0.0070089 0.0485681 0.1443 0.886083
## Residual_HomeRange -0.1516971 0.1957033 -0.7751 0.443464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03042 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1114, Adjusted R-squared: -0.01555
## F-statistic: 0.8775 on 5 and 35 DF, p-value: 0.5062
##
## $pgls_eco_raw
##
## Call:
## pgls(formula = eval(parse(text = paste(pgls_function.out$response_var,
## "res", sep = "_"))) ~ DietBreadth + PopulationDensity_n_km2 +
## GroupSize_filled + Home_Range + Residual_HomeRange, data = HP.noHomo,
## lambda = "ML")
##
## Coefficients:
## (Intercept) DietBreadth PopulationDensity_n_km2
## 1.229604 0.016512 0.026461
## GroupSize_filled Home_Range Residual_HomeRange
## -0.010052 0.007009 -0.151697
#pgls_function(HP, "Sub_res", "BrainVol")Summaries
All neuroanatomical regression plots
Scaling diagnostic multiplots
grid.arrange(HP.Med_res$xy_plot, HP.Brain_res$xy_plot, HPres.Enceph_res$xy_plot,
HPres.Brain_res$xy_plot, HR.body_res$xy_plot,
ncol = 3, top = textGrob("Scaling Diagnostics", gp=gpar(fontsize=18)))
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
grid.arrange(DG.HP_res$xy_plot, FD.HP_res$xy_plot,
Hilus.HP_res$xy_plot, CA2.3.HP_res$xy_plot,
CA3.HP_res$xy_plot, CA2.HP_res$xy_plot,
CA1.HP_res$xy_plot, Sub.HP_res$xy_plot,
ncol=4, top=textGrob("retroHP Subfield Regressions", gp=gpar(fontsize=18)))
Subfield Proportions
Plot the volumetric proportions of each subregion of the hippocampal complex.
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, c("Sub", "CA1", "CA2", "CA3", "Hilus", "FD", "EC","Fibers")) 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
#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="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)## [1] "Alpha = 0.0055"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
)
)
}
createDT(neuro.summ, caption="Neuroanatomical PGLS Regressions Summary", save="neuro.results") 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="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)## [1] "Alpha = 0.0045"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
# 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
PCA <- prcomp(plot.data, scale=TRUE, center=TRUE)
pander::pander(PCA$rotation[,1:4]) # Show PC1-4| PC1 | PC2 | PC3 | PC4 | |
|---|---|---|---|---|
| EC_res | -0.2766 | -0.4391 | 0.1012 | 0.6389 |
| Fibers_res | -0.4483 | 0.002024 | -0.4163 | 0.04329 |
| FD_res | -0.05577 | 0.4477 | 0.5404 | 0.1329 |
| Hilus_res | 0.005114 | 0.5717 | -0.4629 | 0.05015 |
| CA3_res | -0.4739 | -0.1786 | 0.2891 | -0.1639 |
| CA2_res | -0.4072 | -0.1464 | -0.3702 | -0.2326 |
| CA1_res | 0.4678 | -0.2191 | -0.2991 | 0.4474 |
| Sub_res | 0.332 | -0.4225 | 0.0002254 | -0.5373 |
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
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)))You must enable Javascript to view this page properly.
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.
phy.PCA <- phyl.pca(HP$phy, plot.data, method="lambda")
pander::pander(phy.PCA$L[,1:4]) # Show PC1-4| PC1 | PC2 | PC3 | PC4 | |
|---|---|---|---|---|
| EC_res | -0.3342 | 0.4201 | -0.3168 | 0.1267 |
| Fibers_res | -0.4089 | -0.006708 | 0.02873 | -0.5205 |
| FD_res | -0.2924 | -0.5867 | 0.3084 | 0.367 |
| Hilus_res | 0.3159 | -0.8103 | 0.2186 | -0.2173 |
| CA3_res | -0.896 | 0.2314 | -0.1422 | 0.1376 |
| CA2_res | -0.4666 | 0.2318 | 0.1652 | -0.7517 |
| CA1_res | 0.8051 | 0.1888 | -0.5326 | -0.1439 |
| Sub_res | 0.4387 | 0.7011 | 0.5472 | 0.06834 |
phy.PCA$lambda # Show lambda## [1] 0.4992058# Select which vars to put in PCA
summPCA <- summary(phy.PCA)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 0.04926423 0.04254499 0.03056124 0.02850324
## Proportion of Variance 0.33223914 0.24779020 0.12785840 0.11121818
## Cumulative Proportion 0.33223914 0.58002934 0.70788774 0.81910592
## PC5 PC6 PC7 PC8
## Standard deviation 0.02351801 0.02094232 0.01783295 0.003422999
## Proportion of Variance 0.07571611 0.06003948 0.04353449 0.001603987
## Cumulative Proportion 0.89482204 0.95486152 0.99839601 1.000000000# Simple biplot
biplot(phy.PCA$S, phy.PCA$L, xlabs=rep("o", nrow(phy.PCA$S)), main = "Phylogenetic PCA: PC Axes")
# (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)) +
stat_ellipse(geom="polygon", level=0.95, alpha=0.2) +
geom_point(size=3) + labs(title = "Phylogenetic PCA: \nResidual 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: \nResidual 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
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"))You must enable Javascript to view this page properly.
Hierarchical Clustering + Heatmap
Heatmaply
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
# # 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
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
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), ## Optimizing the positions of the tip labels...
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
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), ## Optimizing the positions of the tip labels...
Principal Components: Scaling
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")##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.046301 -0.012976 0.000066 0.007543 0.025552
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 8.6975e-06
## 95.0% CI : (NA, 0.623)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.145203 0.062044 2.3403
## log(eval(parse(text = predictor_var))) 0.483763 0.135102 3.5807
## Pr(>|t|)
## (Intercept) 0.0240930 *
## log(eval(parse(text = predictor_var))) 0.0008821 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01772 on 42 degrees of freedom
## Multiple R-squared: 0.2339, Adjusted R-squared: 0.2156
## F-statistic: 12.82 on 1 and 42 DF, p-value: 0.0008821
## [1] "Top 5 Positive Residuals: PC1_res"
## species res
## 1 Otolemur_crassicaudatus 0.3043376
## 2 Pan_troglodytes_verus 0.2533469
## 3 Piliocolobus_badius 0.2377245
## 4 Avahi_occidentalis 0.2307420
## 5 Galagoides_demidoff 0.2215353
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: PC1_res"
## species res
## 1 Cheirogaleus_major -0.2878240
## 2 Lagothrix_lagotricha -0.2460136
## 3 Erythrocebus_patas -0.2408588
## 4 Saguinus_midas -0.2244946
## 5 Homo_sapiens -0.2183799CA3.PC1_res <- pgls_function(HP, "CA3_res", "PC1") # Strong relationship, duh.##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.021785 -0.005385 0.003205 0.006585 0.030170
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.494
## lower bound : 0.000, p = 0.025655
## upper bound : 1.000, p = 0.0022414
## 95.0% CI : (0.035, 0.906)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.733110 0.040293 18.195
## log(eval(parse(text = predictor_var))) -0.909407 0.072158 -12.603
## Pr(>|t|)
## (Intercept) < 2.2e-16 ***
## log(eval(parse(text = predictor_var))) 8.882e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01051 on 42 degrees of freedom
## Multiple R-squared: 0.7909, Adjusted R-squared: 0.7859
## F-statistic: 158.8 on 1 and 42 DF, p-value: 7.323e-16
## [1] "Top 5 Positive Residuals: CA3_res_res"
## species res
## 1 Homo_sapiens 0.1595291
## 2 Cercopithecus_mitis 0.1530054
## 3 Indri_indri 0.1400509
## 4 Eulemur_fulvus_fulvus 0.1141412
## 5 Propithecus_verreauxi 0.1120924
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA3_res_res"
## species res
## 1 Saimiri_sciureus -0.1923016
## 2 Daubentonia_madagascariensis -0.1785349
## 3 Cebus_albifrons -0.1587581
## 4 Aotus_trivirgatus -0.1477182
## 5 Alouatta_caraya -0.1052251#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.##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.027754 -0.008338 0.000178 0.013720 0.048362
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 2.0392e-05
## 95.0% CI : (NA, 0.757)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.539202 0.067931 7.9375
## log(eval(parse(text = predictor_var))) -0.287866 0.147920 -1.9461
## Pr(>|t|)
## (Intercept) 6.843e-10 ***
## log(eval(parse(text = predictor_var))) 0.05835 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0194 on 42 degrees of freedom
## Multiple R-squared: 0.08272, Adjusted R-squared: 0.06088
## F-statistic: 3.787 on 1 and 42 DF, p-value: 0.05835
## [1] "Top 5 Positive Residuals: CA3_res_res"
## species res
## 1 Homo_sapiens 0.3280246
## 2 Cheirogaleus_major 0.2892187
## 3 Cheirogaleus_medius 0.2425722
## 4 Loris_tardigradus 0.2403541
## 5 Indri_indri 0.2155975
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: CA3_res_res"
## species res
## 1 Otolemur_crassicaudatus -0.3942419
## 2 Pan_troglodytes_verus -0.2527245
## 3 Avahi_occidentalis -0.2428427
## 4 Galago_senegalensis -0.2351995
## 5 Galagoides_demidoff -0.2144560#Result: Significant positive relationship! But somewhat weak R2PC2 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.##
## Call:
## pgls(formula = log(eval(parse(text = response_var))) ~ log(eval(parse(text = predictor_var))),
## data = comp_data, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0293898 -0.0063614 0.0006185 0.0084970 0.0236566
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.283
## lower bound : 0.000, p = 0.2961
## upper bound : 1.000, p = 0.0091761
## 95.0% CI : (NA, 0.924)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.631484 0.041324 15.2811
## log(eval(parse(text = predictor_var))) -0.731923 0.087700 -8.3457
## Pr(>|t|)
## (Intercept) < 2.2e-16 ***
## log(eval(parse(text = predictor_var))) 1.852e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01201 on 42 degrees of freedom
## Multiple R-squared: 0.6238, Adjusted R-squared: 0.6149
## F-statistic: 69.65 on 1 and 42 DF, p-value: 1.852e-10
## [1] "Top 5 Positive Residuals: PC2_res"
## species res
## 1 Alouatta_caraya 0.1544254
## 2 Cercopithecus_mitis 0.1495736
## 3 Piliocolobus_badius 0.1492647
## 4 Pan_troglodytes_verus 0.1372806
## 5 Saguinus_oedipus 0.1334370
## [1] "====================================================="
## [1] "Top 5 Negative Residuals: PC2_res"
## species res
## 1 Daubentonia_madagascariensis -0.2089132
## 2 Cheirogaleus_major -0.2049607
## 3 Varecia_variegata_variegata -0.1895278
## 4 Cheirogaleus_medius -0.1597717
## 5 Microcebus_murinus -0.1471130#Result: No relationship.Principal Component: Eco
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
## PC1 (models + plots)
PC1_eco.out <- PCrun_eco(HP, "PC1"); PC1_eco.out$pgls_eco_summary##
## Call:
## pgls(formula = eval(parse(text = paste(response_var))) ~ DietBreadth +
## PopulationDensity_n_km2 + GroupSize_filled + Home_Range +
## Residual_HomeRange, data = HP.noHomo, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.053448 -0.022765 -0.002138 0.011144 0.069918
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.742
## lower bound : 0.000, p = 0.016433
## upper bound : 1.000, p = 0.033772
## 95.0% CI : (0.148, 0.989)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9526226 0.3232423 2.9471 0.005676 **
## DietBreadth 0.0061024 0.0269044 0.2268 0.821884
## PopulationDensity_n_km2 0.0622104 0.0228364 2.7242 0.009991 **
## GroupSize_filled -0.0081957 0.0208938 -0.3923 0.697249
## Home_Range 0.0235301 0.0413797 0.5686 0.573231
## Residual_HomeRange -0.2242149 0.1590190 -1.4100 0.167368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03078 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3138, Adjusted R-squared: 0.2157
## F-statistic: 3.201 on 5 and 35 DF, p-value: 0.01751PCall_eco_plots(HP, "PC1")
## PC2
PC2_eco <- PCrun_eco(HP, "PC2"); PC2_eco$pgls_eco_summary##
## Call:
## pgls(formula = eval(parse(text = paste(response_var))) ~ DietBreadth +
## PopulationDensity_n_km2 + GroupSize_filled + Home_Range +
## Residual_HomeRange, data = HP.noHomo, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.072543 -0.014605 -0.004042 0.006167 0.047582
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.001975
## 95.0% CI : (NA, 0.370)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6830716 0.2763604 6.0901 5.875e-07 ***
## DietBreadth -0.0133485 0.0273196 -0.4886 0.6282
## PopulationDensity_n_km2 -0.0054613 0.0203782 -0.2680 0.7903
## GroupSize_filled -0.0145846 0.0175890 -0.8292 0.4126
## Home_Range -0.0129972 0.0351789 -0.3695 0.7140
## Residual_HomeRange -0.1348578 0.1417524 -0.9514 0.3479
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02204 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2159, Adjusted R-squared: 0.1039
## F-statistic: 1.928 on 5 and 35 DF, p-value: 0.1145## PC3
HP$data$PC3 <- scales::rescale(phy.PCA$S[,3],c(1,2))
PC3_eco <- PCrun_eco(HP, "PC3"); PC3_eco$pgls_eco_summary##
## Call:
## pgls(formula = eval(parse(text = paste(response_var))) ~ DietBreadth +
## PopulationDensity_n_km2 + GroupSize_filled + Home_Range +
## Residual_HomeRange, data = HP.noHomo, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033493 -0.015185 -0.003380 0.008582 0.052928
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = 0.0007902
## 95.0% CI : (NA, 0.367)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.40583111 0.24946456 5.6354 2.336e-06 ***
## DietBreadth 0.03381417 0.02466077 1.3712 0.1790
## PopulationDensity_n_km2 -0.00071629 0.01839493 -0.0389 0.9692
## GroupSize_filled 0.00521311 0.01587719 0.3283 0.7446
## Home_Range -0.01085088 0.03175527 -0.3417 0.7346
## Residual_HomeRange 0.17654536 0.12795678 1.3797 0.1764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01989 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1552, Adjusted R-squared: 0.03454
## F-statistic: 1.286 on 5 and 35 DF, p-value: 0.2919## PC4
HP$data$PC4 <- scales::rescale(phy.PCA$S[,4],c(1,2))
PC4_eco <- PCrun_eco(HP, "PC4"); PC4_eco$pgls_eco_summary##
## Call:
## pgls(formula = eval(parse(text = paste(response_var))) ~ DietBreadth +
## PopulationDensity_n_km2 + GroupSize_filled + Home_Range +
## Residual_HomeRange, data = HP.noHomo, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049229 -0.017165 0.007111 0.029152 0.050760
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.815
## lower bound : 0.000, p = 4.8048e-07
## upper bound : 1.000, p = 0.074294
## 95.0% CI : (0.463, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8702472 0.3217949 5.8119 1.366e-06 ***
## DietBreadth -0.0235148 0.0261256 -0.9001 0.3742
## PopulationDensity_n_km2 -0.0106191 0.0223858 -0.4744 0.6382
## GroupSize_filled -0.0011542 0.0205401 -0.0562 0.9555
## Home_Range -0.0289135 0.0411179 -0.7032 0.4866
## Residual_HomeRange 0.1652155 0.1568180 1.0535 0.2993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03191 on 35 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.05182, Adjusted R-squared: -0.08363
## F-statistic: 0.3826 on 5 and 35 DF, p-value: 0.8573Inter-subregional Correlations
# 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
# Parallelize to speed up
# nCores <- parallel::detectCores()
# ?l1ou:::plot.l1ou
#install_github("khabbazian/l1ou")
l1ou <- function(tree, dat, model="pBIC", abbrev_species=F, multi=F){
# Shorten species names
if (abbrev_species==T){
abbrevSpecies <- function(spNames) {
species <- lapply(strsplit(spNames, "_"), `[`, 2)
genusAbrv <- substr(spNames, 1, 1)
abrvNames <- paste(paste(genusAbrv,".",sep=""), species, sep="_")
return(abrvNames)
}
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 <- estimate_shift_configuration(tree, Y, criterion="AICc", nCores=4)
## then detect which of these shifts are convergent:
fit_conv <- estimate_convergent_regimes(fit_ind, criterion="AICc", nCores=4)
# 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, edge.width=3,
edge.ann.cex = 1.1, edge.shift.adj = 1.2, label.offset=0.02)
}
# edge.label.pos = 1, x.lim=c(15,1)
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)## [1] TRUEl1ou: Multivariate Analyses
l1ou :PCs 1-4
pcDF <- phy.PCA$S[,1:4]
colnames(pcDF) <- paste("p",colnames(pcDF), sep="")
PCs1_4.pBIC <- l1ou(tree, dat=pcDF, model="pBIC", multi=T)## the new tree is normalized: each tip at distance 1 from the root.
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.04 0.12 0.1 0.01 ) for another list of candidates.
All Subregions
# All non-composite subregions
subs.pBIC <- l1ou(tree, dat=plot.data, model="pBIC", multi=T)## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 8
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.1 0.02 0.07 0.04 0.12 0.01 0.08 0.12 ) for another list of candidates.
Individual Subregions
l1ou: HPres
#out.extra='style="margin: 200px"'
# fig.height=9, fig.width=9
HPres.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,retroHP_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.07 ) for another list of candidates.
l1ou: EC_res
# fig.height=9, fig.width=9
EC.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,EC_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.06 ) for another list of candidates.
l1ou: Fibers_res
# fig.height=9, fig.width=9
Fibers.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,Fibers_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.01 ) for another list of candidates.
l1ou: DG/FD/Hilus_res
# fig.height=9, fig.width=9
DG.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,DG_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
FD.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,FD_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.08 ) for another list of candidates.
Hilus.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,Hilus_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.06 ) for another list of candidates.
l1ou: CA2.3/2/3_res
# Individual subregions
CA2.3.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,CA2.3_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
CA3.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,CA3_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
CA2.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,CA2_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
l1ou: CA1_res
CA1.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,CA1_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
l1ou: Sub_res
Sub.pBIC <- l1ou(tree, dat=dplyr::select(plot.alldata,Sub_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.01 ) for another list of candidates.
l1ou: Ecological Variables
HR.pBIC <- l1ou(tree, dat=dplyr::select(HP$data, Home_Range), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
HRres.pBIC <- l1ou(tree, dat=dplyr::select(HP$data, Home_Range_res), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.06 ) for another list of candidates.
Group.pBIC <- l1ou(tree, dat=dplyr::select(HP$data, GroupSize_filled), model="pBIC")## the new tree is normalized: each tip at distance 1 from the root.
## new Y: matrix of size 44 x 1
## Starting first LASSO (alpha=0) to find a list of candidate configurations.
## Starting second LASSO (alpha= 0.12 ) for another list of candidates.
Tests
Node ages
t <- tree
plot(t)
node.nums <- cbind(t$tip.label, t$node.label)
nodes <- round(heights(t)$start,2) #node.height()
names(nodes) <- nodes
tr <- makeNodeLabel(t, method = "user", nodeList = nodes)
plot(tr, show.node.label = T)
plot.phylo(t, show.node.label = F)
nodelabels(round(heights(t)$start,2))
plotTree.wBars(t)
# Export as nexus and then open in FigTree to get node ages
write.nexus(tree,file="primate_tree.nex")Draw HC Network
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))Interactive Plots
tkid <- tkplot(net) #tkid is the id of the tkplot that will open
l <- tkplot.getcoords(tkid) # grab the coordinates from tkplot
plot(net, edge.arrow.size=.4, vertex.shape="circle",
# vertex.size= V(net)$Size, edge.width=V(net)$Size/5, edge.curved=.1,
# vertex.color=V(net)$color, vertex.label.dist=2, layout=l)
library(visNetwork)
visNetwork(nodes, links, width="100%", height="400px", main=paste("Hippocampal Complex Network:",species))
nodes$shape <- "dot"
nodes$shadow <- TRUE # Nodes will drop shadow
nodes$title <- nodes$Region # Text on click
nodes$label <- nodes$Region # Node label
nodes$size <- nodes$Size # Node size
nodes$borderWidth <- 2 # Node border width
nodes$color.background <- c("slategrey", "tomato", "gold")[nodes$media.type]
nodes$color.border <- "black"
nodes$color.highlight.background <- "orange"
nodes$color.highlight.border <- "darkred"
visNetwork(nodes, links)Save Results
Save the updated data in a new table
createDT(HP$data, caption = "Complete Data",save = "HP_eco_full_results", scrollY = F)