.Rhistory

}
# Plot eco variable function
eco_plots <- function(pgls_function.out,
eco_var,
run_eco.out,
figLabel="",
res_norm_func=function(x){x},
norm_func_plot=log10,
x_prefix=""){
pgls_eco_sum <- run_eco.out$pgls_eco_summary
response_var <- pgls_function.out$response_var
data <- pgls_function.out$new_comp_data$data
brain <- data[,paste(response_var,"res",sep="_")]
eco_x <-data[,paste(eco_var)]
eco.row <- pgls_eco_sum$coefficients[grep(eco_var, rownames(pgls_eco_sum$coefficients)), ]
slope <- as.numeric(eco.row[["Estimate"]])
p <- as.numeric(eco.row[["Pr(>|t|)"]])
#  Initialize plot
# if(eco_var=="Residual_HomeRange"){norm_func <- res_norm_func}
eco_plot <- ggplot(data, aes(y=res_norm_func(brain), x=norm_func_plot(eco_x),
fill=Clade, color=Clade, shape=Clade))  +
theme_classic() + ggtitle(figLabel)
# Conditionally add regression line
if ( p<=0.05){
eco_plot <- eco_plot +
geom_smooth(na.rm=F, inherit.aes=F, method="lm", alpha = .15,
data=data, aes(y=res_norm_func(brain), x=norm_func_plot(eco_x)))
}
# Conditionally add legend
if (eco_var == "Residual_HomeRange"){
eco_plot =  eco_plot + theme(legend.position=c(1.75, 0.35), legend.background = element_rect(fill="whitesmoke"),
plot.title = element_text(hjust= -.1))
} else{
eco_plot = eco_plot + theme( legend.position="none",
plot.title = element_text(hjust= -.1))
}
# Add remaining featuress
eco_plot <- eco_plot + geom_point(size=3, alpha=.7) +
xlab(paste(x_prefix,eco_var,sep="")) + ylab(paste(response_var,"res",sep="_")) +
scale_shape_manual(values=shapes) + scale_color_manual(values=colors) + scale_fill_manual(values=colors)
#geom_abline(slope=SLOPE, intercept=INTERCEPT)
#assign(paste("eco_plot",paste(response_var,"res",sep="_"), "vs",eco_var,sep="."),eco_plot)
return(eco_plot) # %>% ggplotly()
}
# MAKE MORE FLEXIBLE (dependent on which model variables actually used)
all_eco_plots <- function(pgls_function.out,
norm_func_plot = log_scale,
resHR_norm_func_plot=scale_
){
eco=run_eco(pgls_function.out = pgls_function.out)
print(eco$pgls_eco_summary)
p_f.o <-pgls_function.out
# e0 <- eco_plots(eco_var = "Nocturnality", run_eco.out=eco)
e1 <-eco_plots(pgls_function.out = p_f.o, eco_var = "Diet_Breadth", run_eco.out=eco, figLabel="(a)",
norm_func_plot = resHR_norm_func_plot)
#e2 <-eco_plots(p_f.o, eco_var = "HomeRange_km2", run_eco.out=eco)
#e3 <-eco_plots(p_f.o, eco_var = "HomeRange_Indiv_km2", run_eco.out=eco)
e4 <-eco_plots(p_f.o, eco_var = "Population_Density", run_eco.out=eco, figLabel="(b)", norm_func_plot = norm_func_plot)
# e5 <-eco_plots(p_f.o, eco_var = "SocialGrpSize", run_eco.out=eco)
# e6 <-eco_plots(p_f.o, eco_var = "ATWP_DR", run_eco.out=eco)
# e7 <-eco_plots(p_f.o, eco_var = "ATWP_perForaging", run_eco.out=eco)
#e8 <-eco_plots(p_f.o, eco_var = "ATWP_GroupSize", run_eco.out=eco)
#e9 <-eco_plots(p_f.o, eco_var = "ATWP_HR", run_eco.out=eco)
#e10 <-eco_plots(p_f.o, eco_var = "ATWP_perTravel", run_eco.out=eco)
e11 <-eco_plots(p_f.o, eco_var = "Group_Size", run_eco.out=eco, figLabel="(c)", norm_func_plot = norm_func_plot)
e12 <-eco_plots(p_f.o, eco_var = "Home_Range", run_eco.out=eco, figLabel="(d)", norm_func_plot = norm_func_plot)
e13 <-eco_plots(p_f.o, eco_var = "Residual_HomeRange", run_eco.out=eco, figLabel="(e)",
norm_func_plot = resHR_norm_func_plot)
# Additional variables
# e14 <-eco_plots(p_f.o, eco_var = "ActivityCycle", run_eco.out=eco, figLabel="(f)")
# e15 <-eco_plots(p_f.o, eco_var = "HabitatBreadth", run_eco.out=eco, figLabel="(g)")
# e16 <-eco_plots(p_f.o, eco_var = "ATWP_perForaging", run_eco.out=eco, figLabel="(h)")
# e17 <-eco_plots(p_f.o, eco_var = "ATWP_DR", run_eco.out=eco, figLabel="(i)")
mp <- gridExtra::grid.arrange(e1,e4,e11,e12,e13,
# e14,e15,e16,e17,
ncol=3,
top = textGrob(paste("Ecological Plots:",
paste(pgls_function.out$response_var,"res",sep="_")),
gp=gpar(fontsize=18)))
return(mp)
#e6,e7,e10
# mod_vars <- list(eco_vars=c("DietBreadth","PopulationDensity_n_km2","GroupSize_filled", "Home_Range", "Residual_HomeRange"))
# subplot(e1,e4,e11,e12,e13, nrows = 2, shareY=F, shareX=F, titleX = T, titleY=T, margin = 0.1)
# fluidPage(
#   fluidRow( column(4, e1), column(4, e4), column(4, e11)),
#   fluidRow( column(4, e12), column(8, e13) )
# )
#plot_grid( plotlist = list(e1,e4,e11,e12,e13), labels = list())
}
suffix = " (log mm3)"#expression(" (log"~mm^3~')')
# HR_phylo <- contMap(HP$phy, HP$data$Home_Range, type="phylogram", plot=F)
# plot(HR_phylo, leg.txt="Home Range",lwd=3, fsize=c(0.5, 1))
# We know HomeRange scales with BodySize, so perhaps using "extra" HR relative to body weight is a more relavent measure.
## Alternatively could include BodySize as covariate in eco model.
# HP$data$Body.Size <- log10(HP$data$Body_Size)
HR.body_res <- pgls_function(HP, "Home_Range", "Body_Size",
norm_func = log_scale,
x_prefix = "", x_suffix = " (g)",
y_prefix = "", y_suffix = " (km2)")
# Put back into Comp.data
HP$data$Residual_HomeRange <- HR.body_res$residuals
# HP vs. Brain
HP.Brain_res <- pgls_function(HP, "retroHP", "Brain__HP", y_suffix = suffix, x_suffix = suffix)
run_eco(pgls_function.out = HP.Brain_res)
all_eco_plots(pgls_function.out = HP.Brain_res)
# HP vs. Medulla
HP.Med_res <- pgls_function(HP, "retroHP", "Medulla", y_suffix = suffix, x_suffix = suffix)
retroHP.Medulla <- run_eco(HP.Med_res)
all_eco_plots(HP.Med_res)
# Hyp: Relationship with overall brain vol?
#pgls_function(HP, "retroHP_res", "BrainVol")
HP.Brain_res$new_comp_data$data
dim(HP.Brain_res$new_comp_data$data)
# HP vs. Brain
HP.Brain_res <- pgls_function(HP, "retroHP", "Brain__HP", y_suffix = suffix, x_suffix = suffix)
run_eco(pgls_function.out = HP.Brain_res)
all_eco_plots(pgls_function.out = HP.Brain_res)
# HP vs. Medulla
HP.Med_res <- pgls_function(HP, "retroHP", "Medulla", y_suffix = suffix, x_suffix = suffix)
retroHP.Medulla <- run_eco(HP.Med_res)
all_eco_plots(HP.Med_res)
# Hyp: Relationship with overall brain vol?
#pgls_function(HP, "retroHP_res", "BrainVol")
HPres <- HP.Med_res$new_comp_data$data$retroHP_res
names(HPres) <- HP.Med_res$new_comp_data$phy$tip.label
## Get colors from ggplot defaults
gg_color_hue <- function(n) {
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
colz = gg_color_hue(6)
phenogram(HP.Med_res$new_comp_data$phy, x=HPres, spread.labels=T,
fsize=.75, colors=colz, axes=c(time, HPres),
main="Phenogram:\n retroHP vs. medulla PGLS residuals",
xlab="Mya", ylab="retroHP residuals") #spread.cost = c(1, 0)
## Color branches by clades??
#Hyp: As brain size scales up, does retroHP_res becomes proportionally smaller?
HPres.Brain_res <- pgls_function(HP, "retroHP_res", "log10(Brain__HP)", y_suffix = "", norm_func = function(x){x})
retroHP.Brain <- run_eco(pgls_function.out = HPres.Brain_res)
all_eco_plots(HPres.Brain_res)
#Result: No relationship.
#Hyp: As Encephalization scales up, does retroHP_res become proportionally smaller?
Enceph_res <- pgls_function(HP, "BrainVol", "Body__Brain", y_suffix = suffix, x_suffix = suffix)
#Don't normalize residuals
pgls_function(HP, "retroHP_res", "BrainVol_res", x_suffix = "", y_suffix = "")
run_eco(HPres.Enceph_res)
all_eco_plots(HPres.Enceph_res)
#Result: No relationship.
#Hyp: As brain size scales up, does retroHP_res becomes proportionally smaller?
HPres.Brain_res <- pgls_function(HP, "retroHP_res", "log10(Brain__HP)", y_suffix = "", norm_func = function(x){x})
retroHP.Brain <- run_eco(pgls_function.out = HPres.Brain_res)
all_eco_plots(HPres.Brain_res)
all_eco_plots(HPres.Brain_res)
#Hyp: As Encephalization scales up, does retroHP_res become proportionally smaller?
Enceph_res <- pgls_function(HP, "BrainVol", "Body__Brain", y_suffix = suffix, x_suffix = suffix)
#Don't normalize residuals
pgls_function(HP, "retroHP_res", "BrainVol_res", x_suffix = "", y_suffix = "")
run_eco(HPres.Enceph_res)
all_eco_plots(HPres.Enceph_res)
EC.Med_res <-pgls_function(HP, "EC", "Medulla",  y_suffix = suffix, x_suffix = suffix)
EC.Med_res <-pgls_function(HP, "EC", "Medulla",  y_suffix = suffix, x_suffix = suffix)
run_eco(pgls_function.out = EC.Med_res)
all_eco_plots(EC.Med_res)
EC.Brain_res <-pgls_function(HP, "EC", "Brain__EC", y_suffix = suffix, x_suffix = suffix)
EC.HP_res <-pgls_function(HP, "EC", "retroHP", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
EC.retroHP <- run_eco(EC.HP_res)
all_eco_plots(EC.HP_res)
EC.HP_res <-pgls_function(HP, "EC", "retroHP", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
fibers.Med_res <-pgls_function(HP, "Fibers", "Medulla", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
run_eco(fibers.Med_res)
all_eco_plots(fibers.Med_res)
fibers.brain_res <-pgls_function(HP, "Fibers", "Brain__Fib", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
run_eco(fibers.brain_res)
all_eco_plots(fibers.brain_res)
fibers.HP_res <-pgls_function(HP, "Fibers", "retroHP.EC", lambda =  0.001, norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
Fibers.retroHP.EC <- run_eco(fibers.HP_res)
all_eco_plots(fibers.HP_res)
#pgls_function(HP, "Fibers_res", "BrainVol")
DG.HP_res <- pgls_function(HP, "DG", "retroHP__DG", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
DG.retroHP <- run_eco(DG.HP_res)
all_eco_plots(DG.HP_res)
#pgls_function(HP, "DG_res", "BrainVol")
FD.HP_res <- pgls_function(HP, "FD", "retroHP__FD", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
FD.retroHP <- run_eco(FD.HP_res, variable_selection = T)
mp <- all_eco_plots(FD.HP_res)
ggsave("./Figures/FINAL/Fig4.png", mp, dpi=400, height=6, width=9)
# improve model
dat <- nonNA_subset(HP = HP, eco_vars = "Residual_HomeRange", response_var = "FD_res")
model <- caper::pgls(data = dat, formula = FD_res ~ Residual_HomeRange, lambda = "ML")
summary(model)
#pgls_function(HP, "FD_res", "BrainVol")
# improve model
dat <- nonNA_subset(HP = HP, eco_vars = "Residual_HomeRange", response_var = "FD_res")
Hilus.HP_res <- pgls_function(HP, "Hilus", "retroHP__Hil", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
Hilus.retroHP <- run_eco(Hilus.HP_res, norm_func = logTrans_missing)
Hilus.retroHP <- run_eco(Hilus.HP_res, norm_func = logTrans_missing)
all_eco_plots(Hilus.HP_res)
CA2.3.HP_res <- pgls_function(HP, "CA2.3", "retroHP__CA2.3", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
CA2.3.retroHP <- run_eco(CA2.3.HP_res)
all_eco_plots(CA2.3.HP_res)
CA3.HP_res <- pgls_function(HP, "CA3", "retroHP__CA3", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
CA2.HP_res <- pgls_function(HP, "CA2", "retroHP__CA2", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
CA2.retroHP <- run_eco(CA2.HP_res)
all_eco_plots(CA2.HP_res)
#pgls_function(HP, "CA2_res", "BrainVol")# No relationship
#pgls_function(HP, "CA2_res", "retroHP")# No relationship
#pgls_function(HP, "CA2_res", "retroHP_res") # Significant ***
CA1.HP_res <- pgls_function(HP, "CA1", "retroHP__CA1", y_suffix = suffix, x_suffix = suffix)
FD.HP_res <- pgls_function(HP, "FD", "retroHP__FD", norm_func = log10,
y_suffix = suffix, x_suffix = suffix)
# We know HomeRange scales with BodySize, so perhaps using "extra" HR relative to body weight is a more relavent measure.
## Alternatively could include BodySize as covariate in eco model.
# HP$data$Body.Size <- log10(HP$data$Body_Size)
HR.body_res <- pgls_function(HP, "Home_Range", "Body_Size",
norm_func = log_scale,
x_prefix = "", x_suffix = " (g)",
y_prefix = "", y_suffix = " (km2)")
HR.body_res$res_phylo_plot
plot(HR.body_res$res_phylo_plot)
# We know HomeRange scales with BodySize, so perhaps using "extra" HR relative to body weight is a more relavent measure.
## Alternatively could include BodySize as covariate in eco model.
# HP$data$Body.Size <- log10(HP$data$Body_Size)
HR.body_res <- pgls_function(HP, "Home_Range", "Body_Size",
norm_func = log_scale,
x_prefix = "", x_suffix = " (g)",
y_prefix = "", y_suffix = " (km2)")
res_name <- "Residual_HomeRange"
cp1 <- cont_map(HP, res_name)
cont_map <- function(comp_data, res_name){
comp_data <- HP
HPres_data <- as.matrix(comp_data$data[res_name])
names(HPres_data) <- comp_data$phy$tip.label
res_phylo_plot <- contMap(comp_data$phy,
HPres_data,
type="phylogram",
outline = F,
plot=F, fsize=c(0.5, 1))
## simulate correlated trait data
# http://blog.phytools.org/2013/09/plotting-facing-trees-using-phytools.html
return(res_phylo_plot)
}
plot.new()
par(mfrow=c(1,3))
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=c(0.5, 1))
res_name <- "FD_res"
cp1 <- cont_map(HP, res_name)
plot(cp1, leg.txt=paste(res_name),lwd=3, fsize=c(0.5, 1))
res_name <- "CA1_res"
cp2 <- cont_map(HP, res_name)
plot(cp2, leg.txt=paste(res_name),lwd=3, fsize=c(0.5, 1))
fsize <- c(0.5, 1)
fsize <- c(1, 1)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize, type="fan")
plot.new()
par(mfrow=c(1,1))
fsize <- c(.5, 1)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize, type="fan")
plot(cp0, leg.txt=paste(res_name),type="fan")
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
png("./Figures/phenogram1.png")
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
plot.new()
par(mfrow=c(1,3))
fsize <- c(.5, 1)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
png("./Figures/phenogram1.png")
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
plot.new()
par(mfrow=c(1,3))
fsize <- c(.5, 1)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
res_name <- "FD_res"
cp1 <- cont_map(HP, res_name)
plot(cp1, leg.txt=paste(res_name),lwd=3, fsize=fsize)
res_name <- "CA1_res"
cp2 <- cont_map(HP, res_name)
plot(cp2, leg.txt=paste(res_name),lwd=3, fsize=fsize)
png("./Figures/mutli_phenogram.png")
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
res_name <- "FD_res"
cp1 <- cont_map(HP, res_name)
plot(cp1, leg.txt=paste(res_name),lwd=3, fsize=fsize)
res_name <- "CA1_res"
cp2 <- cont_map(HP, res_name)
plot(cp2, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png")
dev.off()
png("./Figures/mutli_phenogram.png", res = 400)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png", res = 400, height = 10, width = 10)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png", res = 400, height = 10, width = 10)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png", res = 400, height = 10, width = 10)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
HPres_data
comp_data <- HP
# comp_data <- HP
HPres_data <- as.matrix(comp_data$data[res_name])
names(HPres_data) <- comp_data$phy$tip.label
res_phylo_plot <- contMap(comp_data$phy,
HPres_data,
type="phylogram",
outline = F,
plot=F, fsize=c(0.5, 1))
fancyTree(tree = comp_data$phy, x=HPres_data)
plot.new()
fancyTree(tree = comp_data$phy, x=HPres_data, type="contmap")
fancyTree(tree = comp_data$phy,
x=HPres_data,
type="contmap",
fsize=c(0.5, 1))
p <- fancyTree(tree = comp_data$phy,
x=HPres_data,
type="contmap",
fsize=c(0.5, 1))
p
p
?png
png("./Figures/mutli_phenogram.png", res = 20)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png", res = 100)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/mutli_phenogram.png", res = 100, height=7, width=8)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
plot.new()
par(mfrow=c(1,1))
plot.new()
par(mfrow=c(1,1))
png("./Figures/mutli_phenogram.png", res = 100, height=7, width=8)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/pheno1.png", res = 100, height=5, width=5)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
png("./Figures/pheno1.png", res = 100)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
# png("./Figures/pheno1.png", res = 100)
tiff("./Figures/pheno1.tif", res = 300)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
# png("./Figures/pheno1.png", res = 100)
tiff("./Figures/pheno1.tif", res = 300, height=10)
res_name <- "Residual_HomeRange"
cp0 <- cont_map(HP, res_name)
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
dev.off()
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
png("./Figures/pheno1.png", res = 100)
dev.off()
plot(cp0, leg.txt=paste(res_name),lwd=3, fsize=fsize)
png("./Figures/pheno1.png", res = 100, height=10)
dev.off()
#pgls_function.out <- Sub.HP_res
extract_neuro <- function(pgls_function.out){
# Get comparison name
var_name <- paste(pgls_function.out$response_var,"res",sep="_")
pgls_summary <- pgls_function.out$pgls_summary
# Get model p.value
f <- pgls_summary$fstatistic
model.p <- round( pf(f[1],f[2],f[3],lower.tail=F) ,3)
# Get other measures
Multiple.R2 <- round(pgls_summary$r.squared,3)
Adj.R2 <- round(pgls_summary$adj.r.squared[1],3)
F.stat <- round(pgls_summary$fstatistic[1],3)
Lambda <- round(pgls_summary$param[2],3)
output <- list(Multiple.R2=Multiple.R2, Adj.R2=Adj.R2, F.stat=F.stat, Lambda=Lambda,
model.pval=model.p)
return(output)
}
neuro.summ <- rbind(
data.frame(Residual="retroHP.med_res", extract_neuro(HP.Med_res)),
data.frame(Residual="retroHP.brain_res", extract_neuro(HP.Brain_res)),
data.frame(Residual="EC.HP_res", extract_neuro(EC.HP_res)),
data.frame(Residual="Fiber.HP_res", extract_neuro(fibers.HP_res)),
data.frame(Residual="DG.HP_res", extract_neuro(DG.HP_res)),
data.frame(Residual="FD.HP_res", extract_neuro(FD.HP_res)),
data.frame(Residual="Hilus.HP_res", extract_neuro(Hilus.HP_res)),
data.frame(Residual="CA2.3.HP_res", extract_neuro(CA2.3.HP_res)),
data.frame(Residual="CA3.HP_res", extract_neuro(CA3.HP_res)),
data.frame(Residual="CA2.HP_res", extract_neuro(CA2.HP_res)),
data.frame(Residual="CA1.HP_res", extract_neuro(CA1.HP_res)),
data.frame(Residual="Sub.HP_res", extract_neuro(Sub.HP_res))
)
# Bonferonni correction for multiple tests
alpha <- round(0.06/11, 4)
paste("Alpha = ", alpha)
neuro.summ$Sig <-  ifelse(neuro.summ$model.pval<0.05 & neuro.summ$model.pval>0.01, "*",
ifelse(neuro.summ$model.pval<0.01 & neuro.summ$model.pval>alpha,"**",
ifelse(neuro.summ$model.pval<alpha,"***", "")))
row.names(neuro.summ)<-NULL
colnames(neuro.summ) <- c("Residual Name",  "Multiple R<sup>2<sup>","Adj. R<sup>2<sup>",
"F stat.", "λ", "Model p-value", "Sig.")
createDT <- function(dt, caption, save=F, scrollY=700){
if (save!=F){
dir.create("./Tables", showWarnings=F)
write.csv(dt, file=file.path("Tables",paste(save,"csv",sep=".")), quote = F, row.names = F)
}
datatable( dt, caption=caption,
extensions = 'Buttons', escape = F,
options = list( dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = scrollY, paging = F
)
)
}
createDT(neuro.summ, caption="Neuroanatomical PGLS Regressions Summary", save="neuro.results")
Sub.HP_res <- pgls_function(HP, "Sub", "retroHP__Sub", y_suffix = suffix, x_suffix = suffix)
Sub.retroHP <- run_eco(Sub.HP_res)
all_eco_plots(Sub.HP_res)
#pgls_function(HP, "Sub_res", "BrainVol")
#pgls_function.out <- Sub.HP_res
extract_eco <- function(pgls_function.out){
# Get comparison name
var_name <- paste(pgls_function.out$response_var,"res",sep="_")
run_eco.out <-run_eco(pgls_function.out)
pgls_eco_summary <-run_eco.out$pgls_eco_summary
# Get model p.value
f <- pgls_eco_summary$fstatistic
model.p <- round( pf(f[1],f[2],f[3],lower.tail=F) ,3)
# Get other measures
Multiple.R2 <- round(pgls_eco_summary$r.squared,3)
Adj.R2 <- round(pgls_eco_summary$adj.r.squared[1],3)
F.stat <- round(pgls_eco_summary$fstatistic[1],3)
Lambda <- round(pgls_eco_summary$param[2],3)
output <- list(Multiple.R2=Multiple.R2, Adj.R2=Adj.R2, F.stat=F.stat, Lambda=Lambda,
model.pval=model.p)
return(output)
}
eco.summ <- rbind(
data.frame(Residual="retroHP.med_res", extract_eco(HP.Med_res)),
data.frame(Residual="retroHP.brain_res", extract_eco(HP.Brain_res)),
data.frame(Residual="EC.HP_res", extract_eco(EC.HP_res)),
data.frame(Residual="Fiber.HP_res", extract_eco(fibers.HP_res)),
data.frame(Residual="DG.HP_res", extract_eco(DG.HP_res)),
data.frame(Residual="FD.HP_res", extract_eco(FD.HP_res)),
data.frame(Residual="Hilus.HP_res", extract_eco(Hilus.HP_res)),
data.frame(Residual="CA2.3.HP_res", extract_eco(CA2.3.HP_res)),
data.frame(Residual="CA3.HP_res", extract_eco(CA3.HP_res)),
data.frame(Residual="CA2.HP_res", extract_eco(CA2.HP_res)),
data.frame(Residual="CA1.HP_res", extract_eco(CA1.HP_res)),
data.frame(Residual="Sub.HP_res", extract_eco(Sub.HP_res))
)
# Bonferonni correction for multiple tests
alpha <- round(0.05/11, 4)
paste("Alpha = ", alpha)
eco.summ$Sig <-  ifelse(eco.summ$model.pval<0.05 & eco.summ$model.pval>0.01, "*",
ifelse(eco.summ$model.pval<0.01 & eco.summ$model.pval>alpha,"**",
ifelse(eco.summ$model.pval<alpha,"***", "")))
row.names(eco.summ)<-NULL
colnames(eco.summ) <- c("Residual Name",  "Multiple R<sup>2<sup>","Adj. R<sup>2<sup>",
"F stat.", "λ", "Model p-value", "Sig.")
createDT(eco.summ, caption="Ecological PGLS Regressions Summary", save = "eco.results")