library(tidyr)
library(phytools)
library(geiger)
library(phylolm)
library(ggplot2)
library(forcats)

data <- read.csv("data/data_hbER_220627.csv")
rownames(data) <- data[,1]
tree <- read.nexus("data/squirrel_mcctree.nex")

data_longlimbs <- data %>% drop_na(hind_long, fore_long) #removes species with missing longlimb data
rownames(data_longlimbs) <- data_longlimbs[,1]
tree_longlimbs_pr <- treedata(phy = tree, data = data_longlimbs, sort=TRUE)$phy

ecology_longlimbs <- data_longlimbs$ecotype4
names(ecology_longlimbs) <- rownames(data_longlimbs) 

col_ecology_longlimbs <- c("#073b4c", "#118ab2", "#ffd166", "#06d6a0") #set colors for ecological groups: blue = flying; orange = ground; green = tree

lnhbER_longlimbs <- log(data_longlimbs$hbER)
names(lnhbER_longlimbs) <- rownames(data_longlimbs) 

lnsize_longlimbs <- log(data_longlimbs$geomean)
names(lnsize_longlimbs) <- rownames(data_longlimbs) 

lnfore_long <- log(data_longlimbs$fore_long)
names(lnfore_long) <- rownames(data_longlimbs) 

lnfore_long_sc_raw <- (phylolm(lnfore_long ~ lnsize_longlimbs, phy = tree_longlimbs_pr, model = "lambda"))$residuals
lnfore_long_sc <- lnfore_long_sc_raw[match(names(lnfore_long),names(lnfore_long_sc_raw))]

lnhind_long <- log(data_longlimbs$hind_long)
names(lnhind_long) <- rownames(data_longlimbs) 

lnhind_long_sc_raw <- (phylolm(lnhind_long ~ lnsize_longlimbs, phy = tree_longlimbs_pr, model = "lambda"))$residuals
lnhind_long_sc <- lnhind_long_sc_raw[match(names(lnhind_long),names(lnhind_long_sc_raw))]

lnolecrpro_long <- log(data_longlimbs$Ulna_ol)
names(lnolecrpro_long) <- rownames(data_longlimbs) 

lnolecrpro_long_sc_raw <- (phylolm(lnolecrpro_long ~ lnsize_longlimbs, phy = tree_longlimbs_pr, model = "lambda"))$residuals
lnolecrpro_long_sc <- lnolecrpro_long_sc_raw[match(names(lnolecrpro_long),names(lnolecrpro_long_sc_raw))]


data_fulllimbs <- data %>% drop_na(hind_full, fore_full) #removes species with missing fulllimb data
rownames(data_fulllimbs) <- data_fulllimbs[,1]
tree_fulllimbs_pr <- treedata(phy = tree, data = data_fulllimbs, sort=TRUE)$phy

ecology_fulllimbs <- data_fulllimbs$ecotype4
names(ecology_fulllimbs) <- rownames(data_fulllimbs) 

col_ecology_fulllimbs <- c("#073b4c", "#118ab2", "#ffd166", "#06d6a0") #set colors for ecological groups: blue = flying; orange = ground; green = tree

lnhbER_fulllimbs <- log(data_fulllimbs$hbER)
names(lnhbER_fulllimbs) <- rownames(data_fulllimbs) 

lnsize_fulllimbs <- log(data_fulllimbs$geomean)
names(lnsize_fulllimbs) <- rownames(data_fulllimbs) 

lnfore_full <- log(data_fulllimbs$fore_full)
names(lnfore_full) <- rownames(data_fulllimbs) 

lnfore_full_sc_raw <- (phylolm(lnfore_full ~ lnsize_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda"))$residuals
lnfore_full_sc <- lnfore_full_sc_raw[match(names(lnfore_full),names(lnfore_full_sc_raw))]

lnhind_full <- log(data_fulllimbs$hind_full)
names(lnhind_full) <- rownames(data_fulllimbs) 

lnhind_full_sc_raw <- (phylolm(lnhind_full ~ lnsize_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda"))$residuals
lnhind_full_sc <- lnhind_full_sc_raw[match(names(lnhind_full),names(lnhind_full_sc_raw))]


#### forelimb_sc long vs hbER ####

### hbER ANOVA ###
hbER_fore_long_pANOVA <- phylolm(lnfore_long_sc~lnhbER_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_long_pANOVA)

### hbER ANCOVA ###
hbER_fore_long_pANCOVA <- phylolm(lnfore_long_sc~lnhbER_longlimbs*ecology_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_long_pANCOVA)

# Make table of coefficients (intercept and slope)
hbER_fore_long_pANOVA_all_int <- hbER_fore_long_pANOVA$coefficients["(Intercept)"]
hbER_fore_long_pANOVA_all_int_boot <- hbER_fore_long_pANOVA$bootstrap[,"(Intercept)"]

hbER_fore_long_pANCOVA_chip_int <- hbER_fore_long_pANCOVA$coefficients["(Intercept)"]
hbER_fore_long_pANCOVA_fly_int <- hbER_fore_long_pANCOVA$coefficients["(Intercept)"]+hbER_fore_long_pANCOVA$coefficients["ecology_longlimbsflying"]
hbER_fore_long_pANCOVA_ground_int <- hbER_fore_long_pANCOVA$coefficients["(Intercept)"]+hbER_fore_long_pANCOVA$coefficients["ecology_longlimbsground"]
hbER_fore_long_pANCOVA_tree_int <- hbER_fore_long_pANCOVA$coefficients["(Intercept)"]+hbER_fore_long_pANCOVA$coefficients["ecology_longlimbstree"]
hbER_fore_long_pANCOVA_chip_int_boot <- hbER_fore_long_pANCOVA$bootstrap[,"(Intercept)"]
hbER_fore_long_pANCOVA_fly_int_boot <- hbER_fore_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANCOVA$bootstrap[,"ecology_longlimbsflying"]
hbER_fore_long_pANCOVA_ground_int_boot <- hbER_fore_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANCOVA$bootstrap[,"ecology_longlimbsground"]
hbER_fore_long_pANCOVA_tree_int_boot <- hbER_fore_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANCOVA$bootstrap[,"ecology_longlimbstree"]

hbER_fore_long_pANOVA_all_slope <- hbER_fore_long_pANOVA$coefficients["lnhbER_longlimbs"]
hbER_fore_long_pANOVA_all_slope_boot <- hbER_fore_long_pANOVA$bootstrap[,"lnhbER_longlimbs"]

hbER_fore_long_pANCOVA_chip_slope <- hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs"]
hbER_fore_long_pANCOVA_fly_slope <- hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbsflying"]
hbER_fore_long_pANCOVA_ground_slope <- hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbsground"]
hbER_fore_long_pANCOVA_tree_slope <- hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbstree"]
hbER_fore_long_pANCOVA_chip_slope_boot <- hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"]
hbER_fore_long_pANCOVA_fly_slope_boot <- hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbsflying"]
hbER_fore_long_pANCOVA_ground_slope_boot <- hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbsground"]
hbER_fore_long_pANCOVA_tree_slope_boot <- hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_fore_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbstree"]

hbER_fore_long_pANCOVA_int_cof <- rbind(
  hbER_fore_long_pANOVA_all_int,
  hbER_fore_long_pANCOVA_chip_int,
  hbER_fore_long_pANCOVA_fly_int,
  hbER_fore_long_pANCOVA_ground_int,
  hbER_fore_long_pANCOVA_tree_int)
hbER_fore_long_pANCOVA_int_95 <- rbind(
  quantile(hbER_fore_long_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
hbER_fore_long_pANCOVA_slope_cof <- rbind(
  hbER_fore_long_pANOVA_all_slope,
  hbER_fore_long_pANCOVA_chip_slope,
  hbER_fore_long_pANCOVA_fly_slope,
  hbER_fore_long_pANCOVA_ground_slope,
  hbER_fore_long_pANCOVA_tree_slope)
hbER_fore_long_pANCOVA_slope_95 <- rbind(
  quantile(hbER_fore_long_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

hbER_fore_long_pANCOVA_cof <- data.frame(hbER_fore_long_pANCOVA_int_cof, hbER_fore_long_pANCOVA_int_95, hbER_fore_long_pANCOVA_slope_cof, hbER_fore_long_pANCOVA_slope_95 )
colnames(hbER_fore_long_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(hbER_fore_long_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(hbER_fore_long_pANCOVA_cof, digits = 2)

# plot lnfore_long_sc ~ lnhbER_longlimbs*ecology
plot(lnfore_long_sc ~ lnhbER_longlimbs, pch = c(24, 21, 22, 23)[factor(ecology_longlimbs)], bg = col_ecology_longlimbs[factor(ecology_longlimbs)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = hbER_fore_long_pANCOVA_chip_int, b = hbER_fore_long_pANCOVA_chip_slope, col = "#073b4c")
abline(a = hbER_fore_long_pANCOVA_fly_int, b = hbER_fore_long_pANCOVA_fly_slope, col = "#118ab2")
abline(a = hbER_fore_long_pANCOVA_ground_int, b = hbER_fore_long_pANCOVA_ground_slope, col = "#ffd166")
abline(a = hbER_fore_long_pANCOVA_tree_int, b = hbER_fore_long_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### hbER ANOVA ###
hbER_fore_long_pANOVA <- phylolm(lnfore_long_sc~lnhbER_longlimbs+ecology_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_long_pANOVA)

# Make table of coefficients (intercept)
hbER_fore_long_pANOVA_chip_int <- hbER_fore_long_pANOVA$coefficients["(Intercept)"]
hbER_fore_long_pANOVA_fly_int <- hbER_fore_long_pANOVA$coefficients["(Intercept)"]+hbER_fore_long_pANOVA$coefficients["ecology_longlimbsflying"]
hbER_fore_long_pANOVA_ground_int <- hbER_fore_long_pANOVA$coefficients["(Intercept)"]+hbER_fore_long_pANOVA$coefficients["ecology_longlimbsground"]
hbER_fore_long_pANOVA_tree_int <- hbER_fore_long_pANOVA$coefficients["(Intercept)"]+hbER_fore_long_pANOVA$coefficients["ecology_longlimbstree"]
hbER_fore_long_pANOVA_chip_int_boot <- hbER_fore_long_pANOVA$bootstrap[,"(Intercept)"]
hbER_fore_long_pANOVA_fly_int_boot <- hbER_fore_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANOVA$bootstrap[,"ecology_longlimbsflying"]
hbER_fore_long_pANOVA_ground_int_boot <- hbER_fore_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANOVA$bootstrap[,"ecology_longlimbsground"]
hbER_fore_long_pANOVA_tree_int_boot <- hbER_fore_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_long_pANOVA$bootstrap[,"ecology_longlimbstree"]

hbER_fore_long_pANOVA_int_cof <- rbind(
  hbER_fore_long_pANOVA_chip_int,
  hbER_fore_long_pANOVA_fly_int,
  hbER_fore_long_pANOVA_ground_int,
  hbER_fore_long_pANOVA_tree_int)
hbER_fore_long_pANOVA_int_95 <- rbind(
  quantile(hbER_fore_long_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_long_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

hbER_fore_long_pANOVA_cof <- data.frame(hbER_fore_long_pANOVA_int_cof, hbER_fore_long_pANOVA_int_95)
colnames(hbER_fore_long_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(hbER_fore_long_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(hbER_fore_long_pANOVA_cof, digits = 2)

fore_long_sc_vp_ecology <- ggplot(data_longlimbs, aes(x=fct_relevel(ecology_longlimbs, "chipmunk", "flying", "ground", "tree"), y= lnfore_long_sc, fill = ecology_longlimbs)) + 
  scale_y_continuous(name="ln fore_sc (mm)", limits = c(-.24,.42), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology_longlimbs) + 
  geom_hline(yintercept=mean(lnfore_long_sc)) + 
  geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 0.5) +
  stat_summary(fun=mean, geom="point", shape=23, size=3, fill = "grey") +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        axis.text.x=element_text(color = "black", size = 13, angle = 0, hjust = .5, vjust = .5),
        axis.text.y=element_text(color = "black", size = 14),
        axis.title.y = element_text(size = 13),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.background = element_blank(),
        legend.position= "none") 
fore_long_sc_vp_ecology


#### hindlimb_sc long vs hbER ####

### hbER ANOVA ###
hbER_hind_long_pANOVA <- phylolm(lnhind_long_sc~lnhbER_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_long_pANOVA)

### hbER ANCOVA ###
hbER_hind_long_pANCOVA <- phylolm(lnhind_long_sc~lnhbER_longlimbs*ecology_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_long_pANCOVA)

# Make table of coefficients (intercept and slope)
hbER_hind_long_pANOVA_all_int <- hbER_hind_long_pANOVA$coefficients["(Intercept)"]
hbER_hind_long_pANOVA_all_int_boot <- hbER_hind_long_pANOVA$bootstrap[,"(Intercept)"]

hbER_hind_long_pANCOVA_chip_int <- hbER_hind_long_pANCOVA$coefficients["(Intercept)"]
hbER_hind_long_pANCOVA_fly_int <- hbER_hind_long_pANCOVA$coefficients["(Intercept)"]+hbER_hind_long_pANCOVA$coefficients["ecology_longlimbsflying"]
hbER_hind_long_pANCOVA_ground_int <- hbER_hind_long_pANCOVA$coefficients["(Intercept)"]+hbER_hind_long_pANCOVA$coefficients["ecology_longlimbsground"]
hbER_hind_long_pANCOVA_tree_int <- hbER_hind_long_pANCOVA$coefficients["(Intercept)"]+hbER_hind_long_pANCOVA$coefficients["ecology_longlimbstree"]
hbER_hind_long_pANCOVA_chip_int_boot <- hbER_hind_long_pANCOVA$bootstrap[,"(Intercept)"]
hbER_hind_long_pANCOVA_fly_int_boot <- hbER_hind_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANCOVA$bootstrap[,"ecology_longlimbsflying"]
hbER_hind_long_pANCOVA_ground_int_boot <- hbER_hind_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANCOVA$bootstrap[,"ecology_longlimbsground"]
hbER_hind_long_pANCOVA_tree_int_boot <- hbER_hind_long_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANCOVA$bootstrap[,"ecology_longlimbstree"]

hbER_hind_long_pANOVA_all_slope <- hbER_hind_long_pANOVA$coefficients["lnhbER_longlimbs"]
hbER_hind_long_pANOVA_all_slope_boot <- hbER_hind_long_pANOVA$bootstrap[,"lnhbER_longlimbs"]

hbER_hind_long_pANCOVA_chip_slope <- hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs"]
hbER_hind_long_pANCOVA_fly_slope <- hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbsflying"]
hbER_hind_long_pANCOVA_ground_slope <- hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbsground"]
hbER_hind_long_pANCOVA_tree_slope <- hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$coefficients["lnhbER_longlimbs:ecology_longlimbstree"]
hbER_hind_long_pANCOVA_chip_slope_boot <- hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"]
hbER_hind_long_pANCOVA_fly_slope_boot <- hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbsflying"]
hbER_hind_long_pANCOVA_ground_slope_boot <- hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbsground"]
hbER_hind_long_pANCOVA_tree_slope_boot <- hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs"] + hbER_hind_long_pANCOVA$bootstrap[,"lnhbER_longlimbs:ecology_longlimbstree"]

hbER_hind_long_pANCOVA_int_cof <- rbind(
  hbER_hind_long_pANOVA_all_int,
  hbER_hind_long_pANCOVA_chip_int,
  hbER_hind_long_pANCOVA_fly_int,
  hbER_hind_long_pANCOVA_ground_int,
  hbER_hind_long_pANCOVA_tree_int)
hbER_hind_long_pANCOVA_int_95 <- rbind(
  quantile(hbER_hind_long_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
hbER_hind_long_pANCOVA_slope_cof <- rbind(
  hbER_hind_long_pANOVA_all_slope,
  hbER_hind_long_pANCOVA_chip_slope,
  hbER_hind_long_pANCOVA_fly_slope,
  hbER_hind_long_pANCOVA_ground_slope,
  hbER_hind_long_pANCOVA_tree_slope)
hbER_hind_long_pANCOVA_slope_95 <- rbind(
  quantile(hbER_hind_long_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

hbER_hind_long_pANCOVA_cof <- data.frame(hbER_hind_long_pANCOVA_int_cof, hbER_hind_long_pANCOVA_int_95, hbER_hind_long_pANCOVA_slope_cof, hbER_hind_long_pANCOVA_slope_95 )
colnames(hbER_hind_long_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(hbER_hind_long_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(hbER_hind_long_pANCOVA_cof, digits = 2)

# plot lnhind_long_sc ~ lnhbER_longlimbs*ecology
plot(lnhind_long_sc ~ lnhbER_longlimbs, pch = c(24, 21, 22, 23)[factor(ecology_longlimbs)], bg = col_ecology_longlimbs[factor(ecology_longlimbs)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = hbER_hind_long_pANCOVA_chip_int, b = hbER_hind_long_pANCOVA_chip_slope, col = "#073b4c")
abline(a = hbER_hind_long_pANCOVA_fly_int, b = hbER_hind_long_pANCOVA_fly_slope, col = "#118ab2")
abline(a = hbER_hind_long_pANCOVA_ground_int, b = hbER_hind_long_pANCOVA_ground_slope, col = "#ffd166")
abline(a = hbER_hind_long_pANCOVA_tree_int, b = hbER_hind_long_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### hbER ANOVA ###
hbER_hind_long_pANOVA <- phylolm(lnhind_long_sc~lnhbER_longlimbs+ecology_longlimbs, phy = tree_longlimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_long_pANOVA)

# Make table of coefficients (intercept)
hbER_hind_long_pANOVA_chip_int <- hbER_hind_long_pANOVA$coefficients["(Intercept)"]
hbER_hind_long_pANOVA_fly_int <- hbER_hind_long_pANOVA$coefficients["(Intercept)"]+hbER_hind_long_pANOVA$coefficients["ecology_longlimbsflying"]
hbER_hind_long_pANOVA_ground_int <- hbER_hind_long_pANOVA$coefficients["(Intercept)"]+hbER_hind_long_pANOVA$coefficients["ecology_longlimbsground"]
hbER_hind_long_pANOVA_tree_int <- hbER_hind_long_pANOVA$coefficients["(Intercept)"]+hbER_hind_long_pANOVA$coefficients["ecology_longlimbstree"]
hbER_hind_long_pANOVA_chip_int_boot <- hbER_hind_long_pANOVA$bootstrap[,"(Intercept)"]
hbER_hind_long_pANOVA_fly_int_boot <- hbER_hind_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANOVA$bootstrap[,"ecology_longlimbsflying"]
hbER_hind_long_pANOVA_ground_int_boot <- hbER_hind_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANOVA$bootstrap[,"ecology_longlimbsground"]
hbER_hind_long_pANOVA_tree_int_boot <- hbER_hind_long_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_long_pANOVA$bootstrap[,"ecology_longlimbstree"]

hbER_hind_long_pANOVA_int_cof <- rbind(
  hbER_hind_long_pANOVA_chip_int,
  hbER_hind_long_pANOVA_fly_int,
  hbER_hind_long_pANOVA_ground_int,
  hbER_hind_long_pANOVA_tree_int)
hbER_hind_long_pANOVA_int_95 <- rbind(
  quantile(hbER_hind_long_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_long_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

hbER_hind_long_pANOVA_cof <- data.frame(hbER_hind_long_pANOVA_int_cof, hbER_hind_long_pANOVA_int_95)
colnames(hbER_hind_long_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(hbER_hind_long_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(hbER_hind_long_pANOVA_cof, digits = 2)

hind_long_sc_vp_ecology <- ggplot(data_longlimbs, aes(x=fct_relevel(ecology_longlimbs, "chipmunk", "flying", "ground", "tree"), y= lnhind_long_sc, fill = ecology_longlimbs)) + 
  scale_y_continuous(name="ln hind_sc (mm)", limits = c(-.24,.42), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology_longlimbs) + 
  geom_hline(yintercept=mean(lnhind_long_sc)) + 
  geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 0.5) +
  stat_summary(fun=mean, geom="point", shape=23, size=3, fill = "grey") +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        axis.text.x=element_text(color = "black", size = 13, angle = 0, hjust = .5, vjust = .5),
        axis.text.y=element_text(color = "black", size = 14),
        axis.title.y = element_text(size = 13),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.background = element_blank(),
        legend.position= "none") 
hind_long_sc_vp_ecology



#### forelimb_sc full vs hbER ####

### hbER ANOVA ###
hbER_fore_full_pANOVA <- phylolm(lnfore_full_sc~lnhbER_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_full_pANOVA)

### hbER ANCOVA ###
hbER_fore_full_pANCOVA <- phylolm(lnfore_full_sc~lnhbER_fulllimbs*ecology_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_full_pANCOVA)

# Make table of coefficients (intercept and slope)
hbER_fore_full_pANOVA_all_int <- hbER_fore_full_pANOVA$coefficients["(Intercept)"]
hbER_fore_full_pANOVA_all_int_boot <- hbER_fore_full_pANOVA$bootstrap[,"(Intercept)"]

hbER_fore_full_pANCOVA_chip_int <- hbER_fore_full_pANCOVA$coefficients["(Intercept)"]
hbER_fore_full_pANCOVA_fly_int <- hbER_fore_full_pANCOVA$coefficients["(Intercept)"]+hbER_fore_full_pANCOVA$coefficients["ecology_fulllimbsflying"]
hbER_fore_full_pANCOVA_ground_int <- hbER_fore_full_pANCOVA$coefficients["(Intercept)"]+hbER_fore_full_pANCOVA$coefficients["ecology_fulllimbsground"]
hbER_fore_full_pANCOVA_tree_int <- hbER_fore_full_pANCOVA$coefficients["(Intercept)"]+hbER_fore_full_pANCOVA$coefficients["ecology_fulllimbstree"]
hbER_fore_full_pANCOVA_chip_int_boot <- hbER_fore_full_pANCOVA$bootstrap[,"(Intercept)"]
hbER_fore_full_pANCOVA_fly_int_boot <- hbER_fore_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANCOVA$bootstrap[,"ecology_fulllimbsflying"]
hbER_fore_full_pANCOVA_ground_int_boot <- hbER_fore_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANCOVA$bootstrap[,"ecology_fulllimbsground"]
hbER_fore_full_pANCOVA_tree_int_boot <- hbER_fore_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANCOVA$bootstrap[,"ecology_fulllimbstree"]

hbER_fore_full_pANOVA_all_slope <- hbER_fore_full_pANOVA$coefficients["lnhbER_fulllimbs"]
hbER_fore_full_pANOVA_all_slope_boot <- hbER_fore_full_pANOVA$bootstrap[,"lnhbER_fulllimbs"]

hbER_fore_full_pANCOVA_chip_slope <- hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs"]
hbER_fore_full_pANCOVA_fly_slope <- hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbsflying"]
hbER_fore_full_pANCOVA_ground_slope <- hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbsground"]
hbER_fore_full_pANCOVA_tree_slope <- hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbstree"]
hbER_fore_full_pANCOVA_chip_slope_boot <- hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"]
hbER_fore_full_pANCOVA_fly_slope_boot <- hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbsflying"]
hbER_fore_full_pANCOVA_ground_slope_boot <- hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbsground"]
hbER_fore_full_pANCOVA_tree_slope_boot <- hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_fore_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbstree"]

hbER_fore_full_pANCOVA_int_cof <- rbind(
  hbER_fore_full_pANOVA_all_int,
  hbER_fore_full_pANCOVA_chip_int,
  hbER_fore_full_pANCOVA_fly_int,
  hbER_fore_full_pANCOVA_ground_int,
  hbER_fore_full_pANCOVA_tree_int)
hbER_fore_full_pANCOVA_int_95 <- rbind(
  quantile(hbER_fore_full_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
hbER_fore_full_pANCOVA_slope_cof <- rbind(
  hbER_fore_full_pANOVA_all_slope,
  hbER_fore_full_pANCOVA_chip_slope,
  hbER_fore_full_pANCOVA_fly_slope,
  hbER_fore_full_pANCOVA_ground_slope,
  hbER_fore_full_pANCOVA_tree_slope)
hbER_fore_full_pANCOVA_slope_95 <- rbind(
  quantile(hbER_fore_full_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

hbER_fore_full_pANCOVA_cof <- data.frame(hbER_fore_full_pANCOVA_int_cof, hbER_fore_full_pANCOVA_int_95, hbER_fore_full_pANCOVA_slope_cof, hbER_fore_full_pANCOVA_slope_95 )
colnames(hbER_fore_full_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(hbER_fore_full_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(hbER_fore_full_pANCOVA_cof, digits = 2)

# plot lnfore_full_sc ~ lnhbER_fulllimbs*ecology
plot(lnfore_full_sc ~ lnhbER_fulllimbs, pch = c(24, 21, 22, 23)[factor(ecology_fulllimbs)], bg = col_ecology_fulllimbs[factor(ecology_fulllimbs)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = hbER_fore_full_pANCOVA_chip_int, b = hbER_fore_full_pANCOVA_chip_slope, col = "#073b4c")
abline(a = hbER_fore_full_pANCOVA_fly_int, b = hbER_fore_full_pANCOVA_fly_slope, col = "#118ab2")
abline(a = hbER_fore_full_pANCOVA_ground_int, b = hbER_fore_full_pANCOVA_ground_slope, col = "#ffd166")
abline(a = hbER_fore_full_pANCOVA_tree_int, b = hbER_fore_full_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### hbER ANOVA ###
hbER_fore_full_pANOVA <- phylolm(lnfore_full_sc~lnhbER_fulllimbs+ecology_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_fore_full_pANOVA)

# Make table of coefficients (intercept)
hbER_fore_full_pANOVA_chip_int <- hbER_fore_full_pANOVA$coefficients["(Intercept)"]
hbER_fore_full_pANOVA_fly_int <- hbER_fore_full_pANOVA$coefficients["(Intercept)"]+hbER_fore_full_pANOVA$coefficients["ecology_fulllimbsflying"]
hbER_fore_full_pANOVA_ground_int <- hbER_fore_full_pANOVA$coefficients["(Intercept)"]+hbER_fore_full_pANOVA$coefficients["ecology_fulllimbsground"]
hbER_fore_full_pANOVA_tree_int <- hbER_fore_full_pANOVA$coefficients["(Intercept)"]+hbER_fore_full_pANOVA$coefficients["ecology_fulllimbstree"]
hbER_fore_full_pANOVA_chip_int_boot <- hbER_fore_full_pANOVA$bootstrap[,"(Intercept)"]
hbER_fore_full_pANOVA_fly_int_boot <- hbER_fore_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANOVA$bootstrap[,"ecology_fulllimbsflying"]
hbER_fore_full_pANOVA_ground_int_boot <- hbER_fore_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANOVA$bootstrap[,"ecology_fulllimbsground"]
hbER_fore_full_pANOVA_tree_int_boot <- hbER_fore_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_fore_full_pANOVA$bootstrap[,"ecology_fulllimbstree"]

hbER_fore_full_pANOVA_int_cof <- rbind(
  hbER_fore_full_pANOVA_chip_int,
  hbER_fore_full_pANOVA_fly_int,
  hbER_fore_full_pANOVA_ground_int,
  hbER_fore_full_pANOVA_tree_int)
hbER_fore_full_pANOVA_int_95 <- rbind(
  quantile(hbER_fore_full_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_fore_full_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

hbER_fore_full_pANOVA_cof <- data.frame(hbER_fore_full_pANOVA_int_cof, hbER_fore_full_pANOVA_int_95)
colnames(hbER_fore_full_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(hbER_fore_full_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(hbER_fore_full_pANOVA_cof, digits = 2)

fore_full_sc_vp_ecology <- ggplot(data_fulllimbs, aes(x=fct_relevel(ecology_fulllimbs, "chipmunk", "flying", "ground", "tree"), y= lnfore_full_sc, fill = ecology_fulllimbs)) + 
  scale_y_continuous(name="ln fore_sc (mm)", limits = c(-.16,.21), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology_fulllimbs) + 
  geom_hline(yintercept=mean(lnfore_full_sc)) + 
  geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 0.5) +
  stat_summary(fun=mean, geom="point", shape=23, size=3, fill = "grey") +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        axis.text.x=element_text(color = "black", size = 13, angle = 0, hjust = .5, vjust = .5),
        axis.text.y=element_text(color = "black", size = 14),
        axis.title.y = element_text(size = 13),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.background = element_blank(),
        legend.position= "none") 
fore_full_sc_vp_ecology


#### hindlimb_sc full vs hbER ####

### hbER ANOVA ###
hbER_hind_full_pANOVA <- phylolm(lnhind_full_sc~lnhbER_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_full_pANOVA)

### hbER ANCOVA ###
hbER_hind_full_pANCOVA <- phylolm(lnhind_full_sc~lnhbER_fulllimbs*ecology_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_full_pANCOVA)

# Make table of coefficients (intercept and slope)
hbER_hind_full_pANOVA_all_int <- hbER_hind_full_pANOVA$coefficients["(Intercept)"]
hbER_hind_full_pANOVA_all_int_boot <- hbER_hind_full_pANOVA$bootstrap[,"(Intercept)"]

hbER_hind_full_pANCOVA_chip_int <- hbER_hind_full_pANCOVA$coefficients["(Intercept)"]
hbER_hind_full_pANCOVA_fly_int <- hbER_hind_full_pANCOVA$coefficients["(Intercept)"]+hbER_hind_full_pANCOVA$coefficients["ecology_fulllimbsflying"]
hbER_hind_full_pANCOVA_ground_int <- hbER_hind_full_pANCOVA$coefficients["(Intercept)"]+hbER_hind_full_pANCOVA$coefficients["ecology_fulllimbsground"]
hbER_hind_full_pANCOVA_tree_int <- hbER_hind_full_pANCOVA$coefficients["(Intercept)"]+hbER_hind_full_pANCOVA$coefficients["ecology_fulllimbstree"]
hbER_hind_full_pANCOVA_chip_int_boot <- hbER_hind_full_pANCOVA$bootstrap[,"(Intercept)"]
hbER_hind_full_pANCOVA_fly_int_boot <- hbER_hind_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANCOVA$bootstrap[,"ecology_fulllimbsflying"]
hbER_hind_full_pANCOVA_ground_int_boot <- hbER_hind_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANCOVA$bootstrap[,"ecology_fulllimbsground"]
hbER_hind_full_pANCOVA_tree_int_boot <- hbER_hind_full_pANCOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANCOVA$bootstrap[,"ecology_fulllimbstree"]

hbER_hind_full_pANOVA_all_slope <- hbER_hind_full_pANOVA$coefficients["lnhbER_fulllimbs"]
hbER_hind_full_pANOVA_all_slope_boot <- hbER_hind_full_pANOVA$bootstrap[,"lnhbER_fulllimbs"]

hbER_hind_full_pANCOVA_chip_slope <- hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs"]
hbER_hind_full_pANCOVA_fly_slope <- hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbsflying"]
hbER_hind_full_pANCOVA_ground_slope <- hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbsground"]
hbER_hind_full_pANCOVA_tree_slope <- hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$coefficients["lnhbER_fulllimbs:ecology_fulllimbstree"]
hbER_hind_full_pANCOVA_chip_slope_boot <- hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"]
hbER_hind_full_pANCOVA_fly_slope_boot <- hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbsflying"]
hbER_hind_full_pANCOVA_ground_slope_boot <- hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbsground"]
hbER_hind_full_pANCOVA_tree_slope_boot <- hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs"] + hbER_hind_full_pANCOVA$bootstrap[,"lnhbER_fulllimbs:ecology_fulllimbstree"]

hbER_hind_full_pANCOVA_int_cof <- rbind(
  hbER_hind_full_pANOVA_all_int,
  hbER_hind_full_pANCOVA_chip_int,
  hbER_hind_full_pANCOVA_fly_int,
  hbER_hind_full_pANCOVA_ground_int,
  hbER_hind_full_pANCOVA_tree_int)
hbER_hind_full_pANCOVA_int_95 <- rbind(
  quantile(hbER_hind_full_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
hbER_hind_full_pANCOVA_slope_cof <- rbind(
  hbER_hind_full_pANOVA_all_slope,
  hbER_hind_full_pANCOVA_chip_slope,
  hbER_hind_full_pANCOVA_fly_slope,
  hbER_hind_full_pANCOVA_ground_slope,
  hbER_hind_full_pANCOVA_tree_slope)
hbER_hind_full_pANCOVA_slope_95 <- rbind(
  quantile(hbER_hind_full_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

hbER_hind_full_pANCOVA_cof <- data.frame(hbER_hind_full_pANCOVA_int_cof, hbER_hind_full_pANCOVA_int_95, hbER_hind_full_pANCOVA_slope_cof, hbER_hind_full_pANCOVA_slope_95 )
colnames(hbER_hind_full_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(hbER_hind_full_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(hbER_hind_full_pANCOVA_cof, digits = 2)

# plot lnhind_full_sc ~ lnhbER_fulllimbs*ecology
plot(lnhind_full_sc ~ lnhbER_fulllimbs, pch = c(24, 21, 22, 23)[factor(ecology_fulllimbs)], bg = col_ecology_fulllimbs[factor(ecology_fulllimbs)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = hbER_hind_full_pANCOVA_chip_int, b = hbER_hind_full_pANCOVA_chip_slope, col = "#073b4c")
abline(a = hbER_hind_full_pANCOVA_fly_int, b = hbER_hind_full_pANCOVA_fly_slope, col = "#118ab2")
abline(a = hbER_hind_full_pANCOVA_ground_int, b = hbER_hind_full_pANCOVA_ground_slope, col = "#ffd166")
abline(a = hbER_hind_full_pANCOVA_tree_int, b = hbER_hind_full_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### hbER ANOVA ###
hbER_hind_full_pANOVA <- phylolm(lnhind_full_sc~lnhbER_fulllimbs+ecology_fulllimbs, phy = tree_fulllimbs_pr, model = "lambda", boot = 1000)
summary(hbER_hind_full_pANOVA)

# Make table of coefficients (intercept)
hbER_hind_full_pANOVA_chip_int <- hbER_hind_full_pANOVA$coefficients["(Intercept)"]
hbER_hind_full_pANOVA_fly_int <- hbER_hind_full_pANOVA$coefficients["(Intercept)"]+hbER_hind_full_pANOVA$coefficients["ecology_fulllimbsflying"]
hbER_hind_full_pANOVA_ground_int <- hbER_hind_full_pANOVA$coefficients["(Intercept)"]+hbER_hind_full_pANOVA$coefficients["ecology_fulllimbsground"]
hbER_hind_full_pANOVA_tree_int <- hbER_hind_full_pANOVA$coefficients["(Intercept)"]+hbER_hind_full_pANOVA$coefficients["ecology_fulllimbstree"]
hbER_hind_full_pANOVA_chip_int_boot <- hbER_hind_full_pANOVA$bootstrap[,"(Intercept)"]
hbER_hind_full_pANOVA_fly_int_boot <- hbER_hind_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANOVA$bootstrap[,"ecology_fulllimbsflying"]
hbER_hind_full_pANOVA_ground_int_boot <- hbER_hind_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANOVA$bootstrap[,"ecology_fulllimbsground"]
hbER_hind_full_pANOVA_tree_int_boot <- hbER_hind_full_pANOVA$bootstrap[,"(Intercept)"] + hbER_hind_full_pANOVA$bootstrap[,"ecology_fulllimbstree"]

hbER_hind_full_pANOVA_int_cof <- rbind(
  hbER_hind_full_pANOVA_chip_int,
  hbER_hind_full_pANOVA_fly_int,
  hbER_hind_full_pANOVA_ground_int,
  hbER_hind_full_pANOVA_tree_int)
hbER_hind_full_pANOVA_int_95 <- rbind(
  quantile(hbER_hind_full_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(hbER_hind_full_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

hbER_hind_full_pANOVA_cof <- data.frame(hbER_hind_full_pANOVA_int_cof, hbER_hind_full_pANOVA_int_95)
colnames(hbER_hind_full_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(hbER_hind_full_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(hbER_hind_full_pANOVA_cof, digits = 2)

hind_full_sc_vp_ecology <- ggplot(data_fulllimbs, aes(x=fct_relevel(ecology_fulllimbs, "chipmunk", "flying", "ground", "tree"), y= lnhind_full_sc, fill = ecology_fulllimbs)) + 
  scale_y_continuous(name="ln hind_sc (mm)", limits = c(-.24,.42), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology_fulllimbs) + 
  geom_hline(yintercept=mean(lnhind_full_sc)) + 
  geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 0.5) +
  stat_summary(fun=mean, geom="point", shape=23, size=3, fill = "grey") +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        axis.text.x=element_text(color = "black", size = 13, angle = 0, hjust = .5, vjust = .5),
        axis.text.y=element_text(color = "black", size = 14),
        axis.title.y = element_text(size = 13),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.background = element_blank(),
        legend.position= "none") 
hind_full_sc_vp_ecology