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")

#remove species in phylogeny that aren't present in dataset
tree_pr <- treedata(phy = tree, data = data, sort=TRUE)$phy

#prune data to check if all species are in tree
#data_mean_pr <- data.frame(treedata(phy = tree_pr, data = data, sort=TRUE)$data)

size <- (as.numeric(as.character(data$geomean))) #save size as vector
names(size) <- rownames(data) #give species names for size vector
lnsize <- log(size) #natural log transform size vector (for stats purposes)

hbER <- (as.numeric(as.character(data$hbER))) #save hbER as vector
names(hbER) <- rownames(data) #give species names for hbER vector
lnhbER <- log(hbER) #natural log transform hbER vector (for stats purposes)

ecology4 <- data$ecotype4 #save ecology as factor
names(ecology4) <- rownames(data) #give species names for ecology factor

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

ecology4_simmap <- make.simmap(tree_pr, ecology4, nsim=1)
plot(ecology4_simmap, colors = setNames(c(col_ecology4), levels(as.factor(ecology4))), fsize = .5)


headER <- (as.numeric(as.character(data$headER)))
names(headER) <- rownames(data)
lnheadER <- log(headER)
AEI_C <- (as.numeric(as.character(data$AEI_C)))
names(AEI_C) <- rownames(data)
lnAEI_C <- log(AEI_C)
AEI_T <- (as.numeric(as.character(data$AEI_T)))
names(AEI_T) <- rownames(data)
lnAEI_T <- log(AEI_T)
AEI_L <- (as.numeric(as.character(data$AEI_L)))
names(AEI_L) <- rownames(data)
lnAEI_L <- log(AEI_L)
AEI_S <- (as.numeric(as.character(data$AEI_S)))
names(AEI_S) <- rownames(data)
lnAEI_S <- log(AEI_S)
rib <- (as.numeric(as.character(data$avgRibLength)))
names(rib) <- rownames(data) #give species names for hbER vector
lnrib <- log(rib)
lnrib_sc_raw <- (phylolm(lnrib~lnsize, phy = tree_pr, model = "lambda"))$residuals
lnrib_sc <- lnrib_sc_raw[match(names(lnrib),names(lnrib_sc_raw))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!!


##### size #####
size_pANOVA <- phylolm(lnsize~ecology4, phy = tree_pr, model = "lambda", boot = 1000) #ANOVA test for size
summary(size_pANOVA)


# Make table of coefficients (intercept)
size_pANOVA_chip_int <- size_pANOVA$coefficients["(Intercept)"]
size_pANOVA_fly_int <- size_pANOVA$coefficients["(Intercept)"]+size_pANOVA$coefficients["ecology4flying"]
size_pANOVA_ground_int <- size_pANOVA$coefficients["(Intercept)"]+size_pANOVA$coefficients["ecology4ground"]
size_pANOVA_tree_int <- size_pANOVA$coefficients["(Intercept)"]+size_pANOVA$coefficients["ecology4tree"]
size_pANOVA_chip_int_boot <- size_pANOVA$bootstrap[,"(Intercept)"]
size_pANOVA_fly_int_boot <- size_pANOVA$bootstrap[,"(Intercept)"] + size_pANOVA$bootstrap[,"ecology4flying"]
size_pANOVA_ground_int_boot <- size_pANOVA$bootstrap[,"(Intercept)"] + size_pANOVA$bootstrap[,"ecology4ground"]
size_pANOVA_tree_int_boot <- size_pANOVA$bootstrap[,"(Intercept)"] + size_pANOVA$bootstrap[,"ecology4tree"]

size_pANOVA_int_cof <- rbind(
  size_pANOVA_chip_int,
  size_pANOVA_fly_int,
  size_pANOVA_ground_int,
  size_pANOVA_tree_int)
size_pANOVA_int_95 <- rbind(
  quantile(size_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_pANOVA_cof <- data.frame(size_pANOVA_int_cof, size_pANOVA_int_95)
colnames(size_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_pANOVA_cof, digits = 3)

size_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chipmunk", "flying", "ground", "tree"), y= lnsize, fill = ecology4)) + 
  scale_y_continuous(name="ln size (mm)", limits = c(2.7,4.5), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnsize)) + 
  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") 
size_vp_ecology


##### hbER #####

### hbER ANOVA ###
size_hbER_pANOVA <- phylolm(lnhbER~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_hbER_pANOVA)

### hbER ANCOVA ###
size_hbER_pANCOVA <- phylolm(lnhbER~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_hbER_pANCOVA)

# Make table of coefficients (intercept and slope)
size_hbER_pANOVA_all_int <- size_hbER_pANOVA$coefficients["(Intercept)"]
size_hbER_pANOVA_all_int_boot <- size_hbER_pANOVA$bootstrap[,"(Intercept)"]

size_hbER_pANCOVA_chip_int <- size_hbER_pANCOVA$coefficients["(Intercept)"]
size_hbER_pANCOVA_fly_int <- size_hbER_pANCOVA$coefficients["(Intercept)"]+size_hbER_pANCOVA$coefficients["ecology4flying"]
size_hbER_pANCOVA_ground_int <- size_hbER_pANCOVA$coefficients["(Intercept)"]+size_hbER_pANCOVA$coefficients["ecology4ground"]
size_hbER_pANCOVA_tree_int <- size_hbER_pANCOVA$coefficients["(Intercept)"]+size_hbER_pANCOVA$coefficients["ecology4tree"]
size_hbER_pANCOVA_chip_int_boot <- size_hbER_pANCOVA$bootstrap[,"(Intercept)"]
size_hbER_pANCOVA_fly_int_boot <- size_hbER_pANCOVA$bootstrap[,"(Intercept)"] + size_hbER_pANCOVA$bootstrap[,"ecology4flying"]
size_hbER_pANCOVA_ground_int_boot <- size_hbER_pANCOVA$bootstrap[,"(Intercept)"] + size_hbER_pANCOVA$bootstrap[,"ecology4ground"]
size_hbER_pANCOVA_tree_int_boot <- size_hbER_pANCOVA$bootstrap[,"(Intercept)"] + size_hbER_pANCOVA$bootstrap[,"ecology4tree"]

size_hbER_pANOVA_all_slope <- size_hbER_pANOVA$coefficients["lnsize"]
size_hbER_pANOVA_all_slope_boot <- size_hbER_pANOVA$bootstrap[,"lnsize"]

size_hbER_pANCOVA_chip_slope <- size_hbER_pANCOVA$coefficients["lnsize"]
size_hbER_pANCOVA_fly_slope <- size_hbER_pANCOVA$coefficients["lnsize"] + size_hbER_pANCOVA$coefficients["lnsize:ecology4flying"]
size_hbER_pANCOVA_ground_slope <- size_hbER_pANCOVA$coefficients["lnsize"] + size_hbER_pANCOVA$coefficients["lnsize:ecology4ground"]
size_hbER_pANCOVA_tree_slope <- size_hbER_pANCOVA$coefficients["lnsize"] + size_hbER_pANCOVA$coefficients["lnsize:ecology4tree"]
size_hbER_pANCOVA_chip_slope_boot <- size_hbER_pANCOVA$bootstrap[,"lnsize"]
size_hbER_pANCOVA_fly_slope_boot <- size_hbER_pANCOVA$bootstrap[,"lnsize"] + size_hbER_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_hbER_pANCOVA_ground_slope_boot <- size_hbER_pANCOVA$bootstrap[,"lnsize"] + size_hbER_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_hbER_pANCOVA_tree_slope_boot <- size_hbER_pANCOVA$bootstrap[,"lnsize"] + size_hbER_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_hbER_pANCOVA_int_cof <- rbind(
  size_hbER_pANOVA_all_int,
  size_hbER_pANCOVA_chip_int,
  size_hbER_pANCOVA_fly_int,
  size_hbER_pANCOVA_ground_int,
  size_hbER_pANCOVA_tree_int)
size_hbER_pANCOVA_int_95 <- rbind(
  quantile(size_hbER_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_hbER_pANCOVA_slope_cof <- rbind(
  size_hbER_pANOVA_all_slope,
  size_hbER_pANCOVA_chip_slope,
  size_hbER_pANCOVA_fly_slope,
  size_hbER_pANCOVA_ground_slope,
  size_hbER_pANCOVA_tree_slope)
size_hbER_pANCOVA_slope_95 <- rbind(
  quantile(size_hbER_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_hbER_pANCOVA_cof <- data.frame(size_hbER_pANCOVA_int_cof, size_hbER_pANCOVA_int_95, size_hbER_pANCOVA_slope_cof, size_hbER_pANCOVA_slope_95 )
colnames(size_hbER_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_hbER_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_hbER_pANCOVA_cof, digits = 2)

# plot lnhbER ~ lnsize*ecology
plot(lnhbER ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_hbER_pANCOVA_chip_int, b = size_hbER_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_hbER_pANCOVA_fly_int, b = size_hbER_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_hbER_pANCOVA_ground_int, b = size_hbER_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_hbER_pANCOVA_tree_int, b = size_hbER_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### hbER ANOVA ###
size_hbER_pANOVA <- phylolm(lnhbER~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_hbER_pANOVA)

# Make table of coefficients (intercept)
size_hbER_pANOVA_chip_int <- size_hbER_pANOVA$coefficients["(Intercept)"]
size_hbER_pANOVA_fly_int <- size_hbER_pANOVA$coefficients["(Intercept)"]+size_hbER_pANOVA$coefficients["ecology4flying"]
size_hbER_pANOVA_ground_int <- size_hbER_pANOVA$coefficients["(Intercept)"]+size_hbER_pANOVA$coefficients["ecology4ground"]
size_hbER_pANOVA_tree_int <- size_hbER_pANOVA$coefficients["(Intercept)"]+size_hbER_pANOVA$coefficients["ecology4tree"]
size_hbER_pANOVA_chip_int_boot <- size_hbER_pANOVA$bootstrap[,"(Intercept)"]
size_hbER_pANOVA_fly_int_boot <- size_hbER_pANOVA$bootstrap[,"(Intercept)"] + size_hbER_pANOVA$bootstrap[,"ecology4flying"]
size_hbER_pANOVA_ground_int_boot <- size_hbER_pANOVA$bootstrap[,"(Intercept)"] + size_hbER_pANOVA$bootstrap[,"ecology4ground"]
size_hbER_pANOVA_tree_int_boot <- size_hbER_pANOVA$bootstrap[,"(Intercept)"] + size_hbER_pANOVA$bootstrap[,"ecology4tree"]

size_hbER_pANOVA_int_cof <- rbind(
  size_hbER_pANOVA_chip_int,
  size_hbER_pANOVA_fly_int,
  size_hbER_pANOVA_ground_int,
  size_hbER_pANOVA_tree_int)
size_hbER_pANOVA_int_95 <- rbind(
  quantile(size_hbER_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_hbER_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_hbER_pANOVA_cof <- data.frame(size_hbER_pANOVA_int_cof, size_hbER_pANOVA_int_95)
colnames(size_hbER_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_hbER_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_hbER_pANOVA_cof, digits = 2)

# Plot hbER ~ ecotype
hbER_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnhbER, fill = ecology4)) + 
  scale_y_continuous(name="ln hbER (mm)", limits = c(1.5,1.85), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnhbER)) + 
  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") 
hbER_vp_ecology

hbER_res <- hbER_res[match(names(hbER),names(hbER_res))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!!


##### headER #####

### headER ANOVA ###
size_headER_pANOVA <- phylolm(lnheadER~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_headER_pANOVA)

### headER ANCOVA ###
size_headER_pANCOVA <- phylolm(lnheadER~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_headER_pANCOVA)

# Make table of coefficients (intercept and slope)
size_headER_pANOVA_all_int <- size_headER_pANOVA$coefficients["(Intercept)"]
size_headER_pANOVA_all_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"]

size_headER_pANCOVA_chip_int <- size_headER_pANCOVA$coefficients["(Intercept)"]
size_headER_pANCOVA_fly_int <- size_headER_pANCOVA$coefficients["(Intercept)"]+size_headER_pANCOVA$coefficients["ecology4flying"]
size_headER_pANCOVA_ground_int <- size_headER_pANCOVA$coefficients["(Intercept)"]+size_headER_pANCOVA$coefficients["ecology4ground"]
size_headER_pANCOVA_tree_int <- size_headER_pANCOVA$coefficients["(Intercept)"]+size_headER_pANCOVA$coefficients["ecology4tree"]
size_headER_pANCOVA_chip_int_boot <- size_headER_pANCOVA$bootstrap[,"(Intercept)"]
size_headER_pANCOVA_fly_int_boot <- size_headER_pANCOVA$bootstrap[,"(Intercept)"] + size_headER_pANCOVA$bootstrap[,"ecology4flying"]
size_headER_pANCOVA_ground_int_boot <- size_headER_pANCOVA$bootstrap[,"(Intercept)"] + size_headER_pANCOVA$bootstrap[,"ecology4ground"]
size_headER_pANCOVA_tree_int_boot <- size_headER_pANCOVA$bootstrap[,"(Intercept)"] + size_headER_pANCOVA$bootstrap[,"ecology4tree"]

size_headER_pANOVA_all_slope <- size_headER_pANOVA$coefficients["lnsize"]
size_headER_pANOVA_all_slope_boot <- size_headER_pANOVA$bootstrap[,"lnsize"]

size_headER_pANCOVA_chip_slope <- size_headER_pANCOVA$coefficients["lnsize"]
size_headER_pANCOVA_fly_slope <- size_headER_pANCOVA$coefficients["lnsize"] + size_headER_pANCOVA$coefficients["lnsize:ecology4flying"]
size_headER_pANCOVA_ground_slope <- size_headER_pANCOVA$coefficients["lnsize"] + size_headER_pANCOVA$coefficients["lnsize:ecology4ground"]
size_headER_pANCOVA_tree_slope <- size_headER_pANCOVA$coefficients["lnsize"] + size_headER_pANCOVA$coefficients["lnsize:ecology4tree"]
size_headER_pANCOVA_chip_slope_boot <- size_headER_pANCOVA$bootstrap[,"lnsize"]
size_headER_pANCOVA_fly_slope_boot <- size_headER_pANCOVA$bootstrap[,"lnsize"] + size_headER_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_headER_pANCOVA_ground_slope_boot <- size_headER_pANCOVA$bootstrap[,"lnsize"] + size_headER_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_headER_pANCOVA_tree_slope_boot <- size_headER_pANCOVA$bootstrap[,"lnsize"] + size_headER_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_headER_pANCOVA_int_cof <- rbind(
  size_headER_pANOVA_all_int,
  size_headER_pANCOVA_chip_int,
  size_headER_pANCOVA_fly_int,
  size_headER_pANCOVA_ground_int,
  size_headER_pANCOVA_tree_int)
size_headER_pANCOVA_int_95 <- rbind(
  quantile(size_headER_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_headER_pANCOVA_slope_cof <- rbind(
  size_headER_pANOVA_all_slope,
  size_headER_pANCOVA_chip_slope,
  size_headER_pANCOVA_fly_slope,
  size_headER_pANCOVA_ground_slope,
  size_headER_pANCOVA_tree_slope)
size_headER_pANCOVA_slope_95 <- rbind(
  quantile(size_headER_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_headER_pANCOVA_cof <- data.frame(size_headER_pANCOVA_int_cof, size_headER_pANCOVA_int_95, size_headER_pANCOVA_slope_cof, size_headER_pANCOVA_slope_95 )
colnames(size_headER_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_headER_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_headER_pANCOVA_cof, digits = 2)

# plot lnheadER ~ lnsize*ecology
plot(lnheadER ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_headER_pANCOVA_chip_int, b = size_headER_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_headER_pANCOVA_fly_int, b = size_headER_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_headER_pANCOVA_ground_int, b = size_headER_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_headER_pANCOVA_tree_int, b = size_headER_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### headER ANOVA ###
size_headER_pANOVA <- phylolm(lnheadER~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_headER_pANOVA)

# Make table of coefficients (intercept)
size_headER_pANOVA_chip_int <- size_headER_pANOVA$coefficients["(Intercept)"]
size_headER_pANOVA_fly_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4flying"]
size_headER_pANOVA_ground_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4ground"]
size_headER_pANOVA_tree_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4tree"]
size_headER_pANOVA_chip_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"]
size_headER_pANOVA_fly_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4flying"]
size_headER_pANOVA_ground_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4ground"]
size_headER_pANOVA_tree_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4tree"]

size_headER_pANOVA_int_cof <- rbind(
  size_headER_pANOVA_chip_int,
  size_headER_pANOVA_fly_int,
  size_headER_pANOVA_ground_int,
  size_headER_pANOVA_tree_int)
size_headER_pANOVA_int_95 <- rbind(
  quantile(size_headER_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_headER_pANOVA_cof <- data.frame(size_headER_pANOVA_int_cof, size_headER_pANOVA_int_95)
colnames(size_headER_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_headER_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_headER_pANOVA_cof, digits = 2)

# plot lnheadER ~ lnsize*ecology
plot(lnheadER ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_headER_pANCOVA_chip_int, b = size_headER_pANCOVA_chip_slope, col = "#073b4c", lty = 2)
abline(a = size_headER_pANCOVA_fly_int, b = size_headER_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_headER_pANCOVA_ground_int, b = size_headER_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_headER_pANCOVA_tree_int, b = size_headER_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### headER ANOVA ###
size_headER_pANOVA <- phylolm(lnheadER~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_headER_pANOVA)

# Make table of coefficients (intercept)
size_headER_pANOVA_chip_int <- size_headER_pANOVA$coefficients["(Intercept)"]
size_headER_pANOVA_fly_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4flying"]
size_headER_pANOVA_ground_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4ground"]
size_headER_pANOVA_tree_int <- size_headER_pANOVA$coefficients["(Intercept)"]+size_headER_pANOVA$coefficients["ecology4tree"]
size_headER_pANOVA_chip_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"]
size_headER_pANOVA_fly_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4flying"]
size_headER_pANOVA_ground_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4ground"]
size_headER_pANOVA_tree_int_boot <- size_headER_pANOVA$bootstrap[,"(Intercept)"] + size_headER_pANOVA$bootstrap[,"ecology4tree"]

size_headER_pANOVA_int_cof <- rbind(
  size_headER_pANOVA_chip_int,
  size_headER_pANOVA_fly_int,
  size_headER_pANOVA_ground_int,
  size_headER_pANOVA_tree_int)
size_headER_pANOVA_int_95 <- rbind(
  quantile(size_headER_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_headER_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_headER_pANOVA_cof <- data.frame(size_headER_pANOVA_int_cof, size_headER_pANOVA_int_95)
colnames(size_headER_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_headER_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_headER_pANOVA_cof, digits = 3)

# Plot headER ~ ecotype
headER_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnheadER, fill = ecology4)) + 
  scale_y_continuous(name="ln headER (mm)", limits = c(.75,1.25), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnheadER)) + 
  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") 
headER_vp_ecology


##### AEI_C #####

### AEI_C ANOVA ###
size_AEI_C_pANOVA <- phylolm(lnAEI_C~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_C_pANOVA)

### AEI_C ANCOVA ###
size_AEI_C_pANCOVA <- phylolm(lnAEI_C~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_C_pANCOVA)

# Make table of coefficients (intercept and slope)
size_AEI_C_pANOVA_all_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]
size_AEI_C_pANOVA_all_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"]

size_AEI_C_pANCOVA_chip_int <- size_AEI_C_pANCOVA$coefficients["(Intercept)"]
size_AEI_C_pANCOVA_fly_int <- size_AEI_C_pANCOVA$coefficients["(Intercept)"]+size_AEI_C_pANCOVA$coefficients["ecology4flying"]
size_AEI_C_pANCOVA_ground_int <- size_AEI_C_pANCOVA$coefficients["(Intercept)"]+size_AEI_C_pANCOVA$coefficients["ecology4ground"]
size_AEI_C_pANCOVA_tree_int <- size_AEI_C_pANCOVA$coefficients["(Intercept)"]+size_AEI_C_pANCOVA$coefficients["ecology4tree"]
size_AEI_C_pANCOVA_chip_int_boot <- size_AEI_C_pANCOVA$bootstrap[,"(Intercept)"]
size_AEI_C_pANCOVA_fly_int_boot <- size_AEI_C_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANCOVA$bootstrap[,"ecology4flying"]
size_AEI_C_pANCOVA_ground_int_boot <- size_AEI_C_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANCOVA$bootstrap[,"ecology4ground"]
size_AEI_C_pANCOVA_tree_int_boot <- size_AEI_C_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANCOVA$bootstrap[,"ecology4tree"]

size_AEI_C_pANOVA_all_slope <- size_AEI_C_pANOVA$coefficients["lnsize"]
size_AEI_C_pANOVA_all_slope_boot <- size_AEI_C_pANOVA$bootstrap[,"lnsize"]

size_AEI_C_pANCOVA_chip_slope <- size_AEI_C_pANCOVA$coefficients["lnsize"]
size_AEI_C_pANCOVA_fly_slope <- size_AEI_C_pANCOVA$coefficients["lnsize"] + size_AEI_C_pANCOVA$coefficients["lnsize:ecology4flying"]
size_AEI_C_pANCOVA_ground_slope <- size_AEI_C_pANCOVA$coefficients["lnsize"] + size_AEI_C_pANCOVA$coefficients["lnsize:ecology4ground"]
size_AEI_C_pANCOVA_tree_slope <- size_AEI_C_pANCOVA$coefficients["lnsize"] + size_AEI_C_pANCOVA$coefficients["lnsize:ecology4tree"]
size_AEI_C_pANCOVA_chip_slope_boot <- size_AEI_C_pANCOVA$bootstrap[,"lnsize"]
size_AEI_C_pANCOVA_fly_slope_boot <- size_AEI_C_pANCOVA$bootstrap[,"lnsize"] + size_AEI_C_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_AEI_C_pANCOVA_ground_slope_boot <- size_AEI_C_pANCOVA$bootstrap[,"lnsize"] + size_AEI_C_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_AEI_C_pANCOVA_tree_slope_boot <- size_AEI_C_pANCOVA$bootstrap[,"lnsize"] + size_AEI_C_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_AEI_C_pANCOVA_int_cof <- rbind(
  size_AEI_C_pANOVA_all_int,
  size_AEI_C_pANCOVA_chip_int,
  size_AEI_C_pANCOVA_fly_int,
  size_AEI_C_pANCOVA_ground_int,
  size_AEI_C_pANCOVA_tree_int)
size_AEI_C_pANCOVA_int_95 <- rbind(
  quantile(size_AEI_C_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_AEI_C_pANCOVA_slope_cof <- rbind(
  size_AEI_C_pANOVA_all_slope,
  size_AEI_C_pANCOVA_chip_slope,
  size_AEI_C_pANCOVA_fly_slope,
  size_AEI_C_pANCOVA_ground_slope,
  size_AEI_C_pANCOVA_tree_slope)
size_AEI_C_pANCOVA_slope_95 <- rbind(
  quantile(size_AEI_C_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_AEI_C_pANCOVA_cof <- data.frame(size_AEI_C_pANCOVA_int_cof, size_AEI_C_pANCOVA_int_95, size_AEI_C_pANCOVA_slope_cof, size_AEI_C_pANCOVA_slope_95 )
colnames(size_AEI_C_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_AEI_C_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_AEI_C_pANCOVA_cof, digits = 2)

# plot lnAEI_C ~ lnsize*ecology
plot(lnAEI_C ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_C_pANCOVA_chip_int, b = size_AEI_C_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_AEI_C_pANCOVA_fly_int, b = size_AEI_C_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_C_pANCOVA_ground_int, b = size_AEI_C_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_C_pANCOVA_tree_int, b = size_AEI_C_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_C ANOVA ###
size_AEI_C_pANOVA <- phylolm(lnAEI_C~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_C_pANOVA)

# Make table of coefficients (intercept)
size_AEI_C_pANOVA_chip_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]
size_AEI_C_pANOVA_fly_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4flying"]
size_AEI_C_pANOVA_ground_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4ground"]
size_AEI_C_pANOVA_tree_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4tree"]
size_AEI_C_pANOVA_chip_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_C_pANOVA_fly_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_C_pANOVA_ground_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_C_pANOVA_tree_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_C_pANOVA_int_cof <- rbind(
  size_AEI_C_pANOVA_chip_int,
  size_AEI_C_pANOVA_fly_int,
  size_AEI_C_pANOVA_ground_int,
  size_AEI_C_pANOVA_tree_int)
size_AEI_C_pANOVA_int_95 <- rbind(
  quantile(size_AEI_C_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_C_pANOVA_cof <- data.frame(size_AEI_C_pANOVA_int_cof, size_AEI_C_pANOVA_int_95)
colnames(size_AEI_C_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_C_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_C_pANOVA_cof, digits = 2)

# plot lnAEI_C ~ lnsize*ecology
plot(lnAEI_C ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_C_pANCOVA_chip_int, b = size_AEI_C_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_AEI_C_pANCOVA_fly_int, b = size_AEI_C_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_C_pANCOVA_ground_int, b = size_AEI_C_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_C_pANCOVA_tree_int, b = size_AEI_C_pANCOVA_tree_slope, col = "#06d6a0")

### AEI_C ANOVA ###
size_AEI_C_pANOVA <- phylolm(lnAEI_C~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_C_pANOVA)

# Make table of coefficients (intercept)
size_AEI_C_pANOVA_chip_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]
size_AEI_C_pANOVA_fly_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4flying"]
size_AEI_C_pANOVA_ground_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4ground"]
size_AEI_C_pANOVA_tree_int <- size_AEI_C_pANOVA$coefficients["(Intercept)"]+size_AEI_C_pANOVA$coefficients["ecology4tree"]
size_AEI_C_pANOVA_chip_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_C_pANOVA_fly_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_C_pANOVA_ground_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_C_pANOVA_tree_int_boot <- size_AEI_C_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_C_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_C_pANOVA_int_cof <- rbind(
  size_AEI_C_pANOVA_chip_int,
  size_AEI_C_pANOVA_fly_int,
  size_AEI_C_pANOVA_ground_int,
  size_AEI_C_pANOVA_tree_int)
size_AEI_C_pANOVA_int_95 <- rbind(
  quantile(size_AEI_C_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_C_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_C_pANOVA_cof <- data.frame(size_AEI_C_pANOVA_int_cof, size_AEI_C_pANOVA_int_95)
colnames(size_AEI_C_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_C_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_C_pANOVA_cof, digits = 3)

# Plot AEI_C ~ ecotype
AEI_C_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnAEI_C, fill = ecology4)) + 
  scale_y_continuous(name="ln AEI_C (mm)", limits = c(0.9,1.4), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnAEI_C)) + 
  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") 
AEI_C_vp_ecology


##### AEI_T #####
### AEI_T ANCOVA ###

### AEI_T ANOVA ###
size_AEI_T_pANOVA <- phylolm(lnAEI_T~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_T_pANOVA)

### AEI_T ANCOVA ###
size_AEI_T_pANCOVA <- phylolm(lnAEI_T~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_T_pANCOVA)

# Make table of coefficients (intercept and slope)
size_AEI_T_pANOVA_all_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]
size_AEI_T_pANOVA_all_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"]

size_AEI_T_pANCOVA_chip_int <- size_AEI_T_pANCOVA$coefficients["(Intercept)"]
size_AEI_T_pANCOVA_fly_int <- size_AEI_T_pANCOVA$coefficients["(Intercept)"]+size_AEI_T_pANCOVA$coefficients["ecology4flying"]
size_AEI_T_pANCOVA_ground_int <- size_AEI_T_pANCOVA$coefficients["(Intercept)"]+size_AEI_T_pANCOVA$coefficients["ecology4ground"]
size_AEI_T_pANCOVA_tree_int <- size_AEI_T_pANCOVA$coefficients["(Intercept)"]+size_AEI_T_pANCOVA$coefficients["ecology4tree"]
size_AEI_T_pANCOVA_chip_int_boot <- size_AEI_T_pANCOVA$bootstrap[,"(Intercept)"]
size_AEI_T_pANCOVA_fly_int_boot <- size_AEI_T_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANCOVA$bootstrap[,"ecology4flying"]
size_AEI_T_pANCOVA_ground_int_boot <- size_AEI_T_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANCOVA$bootstrap[,"ecology4ground"]
size_AEI_T_pANCOVA_tree_int_boot <- size_AEI_T_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANCOVA$bootstrap[,"ecology4tree"]

size_AEI_T_pANOVA_all_slope <- size_AEI_T_pANOVA$coefficients["lnsize"]
size_AEI_T_pANOVA_all_slope_boot <- size_AEI_T_pANOVA$bootstrap[,"lnsize"]

size_AEI_T_pANCOVA_chip_slope <- size_AEI_T_pANCOVA$coefficients["lnsize"]
size_AEI_T_pANCOVA_fly_slope <- size_AEI_T_pANCOVA$coefficients["lnsize"] + size_AEI_T_pANCOVA$coefficients["lnsize:ecology4flying"]
size_AEI_T_pANCOVA_ground_slope <- size_AEI_T_pANCOVA$coefficients["lnsize"] + size_AEI_T_pANCOVA$coefficients["lnsize:ecology4ground"]
size_AEI_T_pANCOVA_tree_slope <- size_AEI_T_pANCOVA$coefficients["lnsize"] + size_AEI_T_pANCOVA$coefficients["lnsize:ecology4tree"]
size_AEI_T_pANCOVA_chip_slope_boot <- size_AEI_T_pANCOVA$bootstrap[,"lnsize"]
size_AEI_T_pANCOVA_fly_slope_boot <- size_AEI_T_pANCOVA$bootstrap[,"lnsize"] + size_AEI_T_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_AEI_T_pANCOVA_ground_slope_boot <- size_AEI_T_pANCOVA$bootstrap[,"lnsize"] + size_AEI_T_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_AEI_T_pANCOVA_tree_slope_boot <- size_AEI_T_pANCOVA$bootstrap[,"lnsize"] + size_AEI_T_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_AEI_T_pANCOVA_int_cof <- rbind(
  size_AEI_T_pANOVA_all_int,
  size_AEI_T_pANCOVA_chip_int,
  size_AEI_T_pANCOVA_fly_int,
  size_AEI_T_pANCOVA_ground_int,
  size_AEI_T_pANCOVA_tree_int)
size_AEI_T_pANCOVA_int_95 <- rbind(
  quantile(size_AEI_T_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_AEI_T_pANCOVA_slope_cof <- rbind(
  size_AEI_T_pANOVA_all_slope,
  size_AEI_T_pANCOVA_chip_slope,
  size_AEI_T_pANCOVA_fly_slope,
  size_AEI_T_pANCOVA_ground_slope,
  size_AEI_T_pANCOVA_tree_slope)
size_AEI_T_pANCOVA_slope_95 <- rbind(
  quantile(size_AEI_T_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_AEI_T_pANCOVA_cof <- data.frame(size_AEI_T_pANCOVA_int_cof, size_AEI_T_pANCOVA_int_95, size_AEI_T_pANCOVA_slope_cof, size_AEI_T_pANCOVA_slope_95 )
colnames(size_AEI_T_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_AEI_T_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_AEI_T_pANCOVA_cof, digits = 2)

# plot lnAEI_T ~ lnsize*ecology
plot(lnAEI_T ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_T_pANCOVA_chip_int, b = size_AEI_T_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_AEI_T_pANCOVA_fly_int, b = size_AEI_T_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_T_pANCOVA_ground_int, b = size_AEI_T_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_T_pANCOVA_tree_int, b = size_AEI_T_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_T ANOVA ###
size_AEI_T_pANOVA <- phylolm(lnAEI_T~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_T_pANOVA)

# Make table of coefficients (intercept)
size_AEI_T_pANOVA_chip_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]
size_AEI_T_pANOVA_fly_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4flying"]
size_AEI_T_pANOVA_ground_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4ground"]
size_AEI_T_pANOVA_tree_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4tree"]
size_AEI_T_pANOVA_chip_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_T_pANOVA_fly_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_T_pANOVA_ground_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_T_pANOVA_tree_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_T_pANOVA_int_cof <- rbind(
  size_AEI_T_pANOVA_chip_int,
  size_AEI_T_pANOVA_fly_int,
  size_AEI_T_pANOVA_ground_int,
  size_AEI_T_pANOVA_tree_int)
size_AEI_T_pANOVA_int_95 <- rbind(
  quantile(size_AEI_T_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_T_pANOVA_cof <- data.frame(size_AEI_T_pANOVA_int_cof, size_AEI_T_pANOVA_int_95)
colnames(size_AEI_T_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_T_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_T_pANOVA_cof, digits = 2)

# plot lnAEI_T ~ lnsize*ecology
plot(lnAEI_T ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_T_pANCOVA_chip_int, b = size_AEI_T_pANCOVA_chip_slope, col = "#073b4c", lty = 2)
abline(a = size_AEI_T_pANCOVA_fly_int, b = size_AEI_T_pANCOVA_fly_slope, col = "#118ab2", lty =2)
abline(a = size_AEI_T_pANCOVA_ground_int, b = size_AEI_T_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_T_pANCOVA_tree_int, b = size_AEI_T_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_T ANOVA ###
size_AEI_T_pANOVA <- phylolm(lnAEI_T~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_T_pANOVA)

# Make table of coefficients (intercept)
size_AEI_T_pANOVA_chip_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]
size_AEI_T_pANOVA_fly_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4flying"]
size_AEI_T_pANOVA_ground_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4ground"]
size_AEI_T_pANOVA_tree_int <- size_AEI_T_pANOVA$coefficients["(Intercept)"]+size_AEI_T_pANOVA$coefficients["ecology4tree"]
size_AEI_T_pANOVA_chip_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_T_pANOVA_fly_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_T_pANOVA_ground_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_T_pANOVA_tree_int_boot <- size_AEI_T_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_T_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_T_pANOVA_int_cof <- rbind(
  size_AEI_T_pANOVA_chip_int,
  size_AEI_T_pANOVA_fly_int,
  size_AEI_T_pANOVA_ground_int,
  size_AEI_T_pANOVA_tree_int)
size_AEI_T_pANOVA_int_95 <- rbind(
  quantile(size_AEI_T_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_T_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_T_pANOVA_cof <- data.frame(size_AEI_T_pANOVA_int_cof, size_AEI_T_pANOVA_int_95)
colnames(size_AEI_T_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_T_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_T_pANOVA_cof, digits = 3)

# Plot AEI_T ~ ecotype
AEI_T_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnAEI_T, fill = ecology4)) + 
  scale_y_continuous(name="ln AEI_T (mm)", limits = c(1.35,2.05), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnAEI_T)) + 
  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") 
AEI_T_vp_ecology


##### AEI_L #####

### AEI_L ANOVA ###
size_AEI_L_pANOVA <- phylolm(lnAEI_L~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_L_pANOVA)

### AEI_L ANCOVA ###
size_AEI_L_pANCOVA <- phylolm(lnAEI_L~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_L_pANCOVA)

# Make table of coefficients (intercept and slope)
size_AEI_L_pANOVA_all_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]
size_AEI_L_pANOVA_all_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"]

size_AEI_L_pANCOVA_chip_int <- size_AEI_L_pANCOVA$coefficients["(Intercept)"]
size_AEI_L_pANCOVA_fly_int <- size_AEI_L_pANCOVA$coefficients["(Intercept)"]+size_AEI_L_pANCOVA$coefficients["ecology4flying"]
size_AEI_L_pANCOVA_ground_int <- size_AEI_L_pANCOVA$coefficients["(Intercept)"]+size_AEI_L_pANCOVA$coefficients["ecology4ground"]
size_AEI_L_pANCOVA_tree_int <- size_AEI_L_pANCOVA$coefficients["(Intercept)"]+size_AEI_L_pANCOVA$coefficients["ecology4tree"]
size_AEI_L_pANCOVA_chip_int_boot <- size_AEI_L_pANCOVA$bootstrap[,"(Intercept)"]
size_AEI_L_pANCOVA_fly_int_boot <- size_AEI_L_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANCOVA$bootstrap[,"ecology4flying"]
size_AEI_L_pANCOVA_ground_int_boot <- size_AEI_L_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANCOVA$bootstrap[,"ecology4ground"]
size_AEI_L_pANCOVA_tree_int_boot <- size_AEI_L_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANCOVA$bootstrap[,"ecology4tree"]

size_AEI_L_pANOVA_all_slope <- size_AEI_L_pANOVA$coefficients["lnsize"]
size_AEI_L_pANOVA_all_slope_boot <- size_AEI_L_pANOVA$bootstrap[,"lnsize"]

size_AEI_L_pANCOVA_chip_slope <- size_AEI_L_pANCOVA$coefficients["lnsize"]
size_AEI_L_pANCOVA_fly_slope <- size_AEI_L_pANCOVA$coefficients["lnsize"] + size_AEI_L_pANCOVA$coefficients["lnsize:ecology4flying"]
size_AEI_L_pANCOVA_ground_slope <- size_AEI_L_pANCOVA$coefficients["lnsize"] + size_AEI_L_pANCOVA$coefficients["lnsize:ecology4ground"]
size_AEI_L_pANCOVA_tree_slope <- size_AEI_L_pANCOVA$coefficients["lnsize"] + size_AEI_L_pANCOVA$coefficients["lnsize:ecology4tree"]
size_AEI_L_pANCOVA_chip_slope_boot <- size_AEI_L_pANCOVA$bootstrap[,"lnsize"]
size_AEI_L_pANCOVA_fly_slope_boot <- size_AEI_L_pANCOVA$bootstrap[,"lnsize"] + size_AEI_L_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_AEI_L_pANCOVA_ground_slope_boot <- size_AEI_L_pANCOVA$bootstrap[,"lnsize"] + size_AEI_L_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_AEI_L_pANCOVA_tree_slope_boot <- size_AEI_L_pANCOVA$bootstrap[,"lnsize"] + size_AEI_L_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_AEI_L_pANCOVA_int_cof <- rbind(
  size_AEI_L_pANOVA_all_int,
  size_AEI_L_pANCOVA_chip_int,
  size_AEI_L_pANCOVA_fly_int,
  size_AEI_L_pANCOVA_ground_int,
  size_AEI_L_pANCOVA_tree_int)
size_AEI_L_pANCOVA_int_95 <- rbind(
  quantile(size_AEI_L_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_AEI_L_pANCOVA_slope_cof <- rbind(
  size_AEI_L_pANOVA_all_slope,
  size_AEI_L_pANCOVA_chip_slope,
  size_AEI_L_pANCOVA_fly_slope,
  size_AEI_L_pANCOVA_ground_slope,
  size_AEI_L_pANCOVA_tree_slope)
size_AEI_L_pANCOVA_slope_95 <- rbind(
  quantile(size_AEI_L_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_AEI_L_pANCOVA_cof <- data.frame(size_AEI_L_pANCOVA_int_cof, size_AEI_L_pANCOVA_int_95, size_AEI_L_pANCOVA_slope_cof, size_AEI_L_pANCOVA_slope_95 )
colnames(size_AEI_L_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_AEI_L_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_AEI_L_pANCOVA_cof, digits = 2)

# plot lnAEI_L ~ lnsize*ecology
plot(lnAEI_L ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_L_pANCOVA_chip_int, b = size_AEI_L_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_AEI_L_pANCOVA_fly_int, b = size_AEI_L_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_L_pANCOVA_ground_int, b = size_AEI_L_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_L_pANCOVA_tree_int, b = size_AEI_L_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_L ANOVA ###
size_AEI_L_pANOVA <- phylolm(lnAEI_L~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_L_pANOVA)

# Make table of coefficients (intercept)
size_AEI_L_pANOVA_chip_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]
size_AEI_L_pANOVA_fly_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4flying"]
size_AEI_L_pANOVA_ground_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4ground"]
size_AEI_L_pANOVA_tree_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4tree"]
size_AEI_L_pANOVA_chip_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_L_pANOVA_fly_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_L_pANOVA_ground_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_L_pANOVA_tree_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_L_pANOVA_int_cof <- rbind(
  size_AEI_L_pANOVA_chip_int,
  size_AEI_L_pANOVA_fly_int,
  size_AEI_L_pANOVA_ground_int,
  size_AEI_L_pANOVA_tree_int)
size_AEI_L_pANOVA_int_95 <- rbind(
  quantile(size_AEI_L_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_L_pANOVA_cof <- data.frame(size_AEI_L_pANOVA_int_cof, size_AEI_L_pANOVA_int_95)
colnames(size_AEI_L_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_L_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_L_pANOVA_cof, digits = 2)

# plot lnAEI_L ~ lnsize*ecology
plot(lnAEI_L ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_L_pANCOVA_chip_int, b = size_AEI_L_pANCOVA_chip_slope, col = "#073b4c", lty=2)
abline(a = size_AEI_L_pANCOVA_fly_int, b = size_AEI_L_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_L_pANCOVA_ground_int, b = size_AEI_L_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_L_pANCOVA_tree_int, b = size_AEI_L_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_L ANOVA ###
size_AEI_L_pANOVA <- phylolm(lnAEI_L~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_L_pANOVA)

# Make table of coefficients (intercept)
size_AEI_L_pANOVA_chip_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]
size_AEI_L_pANOVA_fly_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4flying"]
size_AEI_L_pANOVA_ground_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4ground"]
size_AEI_L_pANOVA_tree_int <- size_AEI_L_pANOVA$coefficients["(Intercept)"]+size_AEI_L_pANOVA$coefficients["ecology4tree"]
size_AEI_L_pANOVA_chip_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_L_pANOVA_fly_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_L_pANOVA_ground_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_L_pANOVA_tree_int_boot <- size_AEI_L_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_L_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_L_pANOVA_int_cof <- rbind(
  size_AEI_L_pANOVA_chip_int,
  size_AEI_L_pANOVA_fly_int,
  size_AEI_L_pANOVA_ground_int,
  size_AEI_L_pANOVA_tree_int)
size_AEI_L_pANOVA_int_95 <- rbind(
  quantile(size_AEI_L_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_L_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_L_pANOVA_cof <- data.frame(size_AEI_L_pANOVA_int_cof, size_AEI_L_pANOVA_int_95)
colnames(size_AEI_L_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_L_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_L_pANOVA_cof, digits = 3)

# Plot AEI_L ~ ecotype
AEI_L_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnAEI_L, fill = ecology4)) + 
  scale_y_continuous(name="ln AEI_L (mm)", limits = c(1.35,2.1), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnAEI_L)) + 
  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") 
AEI_L_vp_ecology


##### AEI_S #####

### AEI_S ANOVA ###
size_AEI_S_pANOVA <- phylolm(lnAEI_S~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_S_pANOVA)

### AEI_S ANCOVA ###
size_AEI_S_pANCOVA <- phylolm(lnAEI_S~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_S_pANCOVA)

# Make table of coefficients (intercept and slope)
size_AEI_S_pANOVA_all_int <- size_AEI_S_pANOVA$coefficients["(Intercept)"]
size_AEI_S_pANOVA_all_int_boot <- size_AEI_S_pANOVA$bootstrap[,"(Intercept)"]

size_AEI_S_pANCOVA_chip_int <- size_AEI_S_pANCOVA$coefficients["(Intercept)"]
size_AEI_S_pANCOVA_fly_int <- size_AEI_S_pANCOVA$coefficients["(Intercept)"]+size_AEI_S_pANCOVA$coefficients["ecology4flying"]
size_AEI_S_pANCOVA_ground_int <- size_AEI_S_pANCOVA$coefficients["(Intercept)"]+size_AEI_S_pANCOVA$coefficients["ecology4ground"]
size_AEI_S_pANCOVA_tree_int <- size_AEI_S_pANCOVA$coefficients["(Intercept)"]+size_AEI_S_pANCOVA$coefficients["ecology4tree"]
size_AEI_S_pANCOVA_chip_int_boot <- size_AEI_S_pANCOVA$bootstrap[,"(Intercept)"]
size_AEI_S_pANCOVA_fly_int_boot <- size_AEI_S_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANCOVA$bootstrap[,"ecology4flying"]
size_AEI_S_pANCOVA_ground_int_boot <- size_AEI_S_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANCOVA$bootstrap[,"ecology4ground"]
size_AEI_S_pANCOVA_tree_int_boot <- size_AEI_S_pANCOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANCOVA$bootstrap[,"ecology4tree"]

size_AEI_S_pANOVA_all_slope <- size_AEI_S_pANOVA$coefficients["lnsize"]
size_AEI_S_pANOVA_all_slope_boot <- size_AEI_S_pANOVA$bootstrap[,"lnsize"]

size_AEI_S_pANCOVA_chip_slope <- size_AEI_S_pANCOVA$coefficients["lnsize"]
size_AEI_S_pANCOVA_fly_slope <- size_AEI_S_pANCOVA$coefficients["lnsize"] + size_AEI_S_pANCOVA$coefficients["lnsize:ecology4flying"]
size_AEI_S_pANCOVA_ground_slope <- size_AEI_S_pANCOVA$coefficients["lnsize"] + size_AEI_S_pANCOVA$coefficients["lnsize:ecology4ground"]
size_AEI_S_pANCOVA_tree_slope <- size_AEI_S_pANCOVA$coefficients["lnsize"] + size_AEI_S_pANCOVA$coefficients["lnsize:ecology4tree"]
size_AEI_S_pANCOVA_chip_slope_boot <- size_AEI_S_pANCOVA$bootstrap[,"lnsize"]
size_AEI_S_pANCOVA_fly_slope_boot <- size_AEI_S_pANCOVA$bootstrap[,"lnsize"] + size_AEI_S_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_AEI_S_pANCOVA_ground_slope_boot <- size_AEI_S_pANCOVA$bootstrap[,"lnsize"] + size_AEI_S_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_AEI_S_pANCOVA_tree_slope_boot <- size_AEI_S_pANCOVA$bootstrap[,"lnsize"] + size_AEI_S_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_AEI_S_pANCOVA_int_cof <- rbind(
  size_AEI_S_pANOVA_all_int,
  size_AEI_S_pANCOVA_chip_int,
  size_AEI_S_pANCOVA_fly_int,
  size_AEI_S_pANCOVA_ground_int,
  size_AEI_S_pANCOVA_tree_int)
size_AEI_S_pANCOVA_int_95 <- rbind(
  quantile(size_AEI_S_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_AEI_S_pANCOVA_slope_cof <- rbind(
  size_AEI_S_pANOVA_all_slope,
  size_AEI_S_pANCOVA_chip_slope,
  size_AEI_S_pANCOVA_fly_slope,
  size_AEI_S_pANCOVA_ground_slope,
  size_AEI_S_pANCOVA_tree_slope)
size_AEI_S_pANCOVA_slope_95 <- rbind(
  quantile(size_AEI_S_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_AEI_S_pANCOVA_cof <- data.frame(size_AEI_S_pANCOVA_int_cof, size_AEI_S_pANCOVA_int_95, size_AEI_S_pANCOVA_slope_cof, size_AEI_S_pANCOVA_slope_95 )
colnames(size_AEI_S_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_AEI_S_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_AEI_S_pANCOVA_cof, digits = 2)

# plot lnAEI_S ~ lnsize*ecology
plot(lnAEI_S ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_AEI_S_pANCOVA_chip_int, b = size_AEI_S_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_AEI_S_pANCOVA_fly_int, b = size_AEI_S_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_AEI_S_pANCOVA_ground_int, b = size_AEI_S_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_AEI_S_pANCOVA_tree_int, b = size_AEI_S_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### AEI_S ANOVA ###
size_AEI_S_pANOVA <- phylolm(lnAEI_S~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_AEI_S_pANOVA)

# Make table of coefficients (intercept)
size_AEI_S_pANOVA_chip_int <- size_AEI_S_pANOVA$coefficients["(Intercept)"]
size_AEI_S_pANOVA_fly_int <- size_AEI_S_pANOVA$coefficients["(Intercept)"]+size_AEI_S_pANOVA$coefficients["ecology4flying"]
size_AEI_S_pANOVA_ground_int <- size_AEI_S_pANOVA$coefficients["(Intercept)"]+size_AEI_S_pANOVA$coefficients["ecology4ground"]
size_AEI_S_pANOVA_tree_int <- size_AEI_S_pANOVA$coefficients["(Intercept)"]+size_AEI_S_pANOVA$coefficients["ecology4tree"]
size_AEI_S_pANOVA_chip_int_boot <- size_AEI_S_pANOVA$bootstrap[,"(Intercept)"]
size_AEI_S_pANOVA_fly_int_boot <- size_AEI_S_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANOVA$bootstrap[,"ecology4flying"]
size_AEI_S_pANOVA_ground_int_boot <- size_AEI_S_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANOVA$bootstrap[,"ecology4ground"]
size_AEI_S_pANOVA_tree_int_boot <- size_AEI_S_pANOVA$bootstrap[,"(Intercept)"] + size_AEI_S_pANOVA$bootstrap[,"ecology4tree"]

size_AEI_S_pANOVA_int_cof <- rbind(
  size_AEI_S_pANOVA_chip_int,
  size_AEI_S_pANOVA_fly_int,
  size_AEI_S_pANOVA_ground_int,
  size_AEI_S_pANOVA_tree_int)
size_AEI_S_pANOVA_int_95 <- rbind(
  quantile(size_AEI_S_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_AEI_S_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_AEI_S_pANOVA_cof <- data.frame(size_AEI_S_pANOVA_int_cof, size_AEI_S_pANOVA_int_95)
colnames(size_AEI_S_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_AEI_S_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_AEI_S_pANOVA_cof, digits = 2)

# Plot AEI_S ~ ecotype
AEI_S_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnAEI_S, fill = ecology4)) + 
  scale_y_continuous(name="ln AEI_S (mm)", limits = c(.3,1), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnAEI_S)) + 
  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") 
AEI_S_vp_ecology



##### rib_sc #####

### rib_sc ANOVA ###
size_rib_sc_pANOVA <- phylolm(lnrib_sc~lnsize, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_rib_sc_pANOVA)

### rib_sc ANCOVA ###
size_rib_sc_pANCOVA <- phylolm(lnrib_sc~lnsize*ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_rib_sc_pANCOVA)

# Make table of coefficients (intercept and slope)
size_rib_sc_pANOVA_all_int <- size_rib_sc_pANOVA$coefficients["(Intercept)"]
size_rib_sc_pANOVA_all_int_boot <- size_rib_sc_pANOVA$bootstrap[,"(Intercept)"]

size_rib_sc_pANCOVA_chip_int <- size_rib_sc_pANCOVA$coefficients["(Intercept)"]
size_rib_sc_pANCOVA_fly_int <- size_rib_sc_pANCOVA$coefficients["(Intercept)"]+size_rib_sc_pANCOVA$coefficients["ecology4flying"]
size_rib_sc_pANCOVA_ground_int <- size_rib_sc_pANCOVA$coefficients["(Intercept)"]+size_rib_sc_pANCOVA$coefficients["ecology4ground"]
size_rib_sc_pANCOVA_tree_int <- size_rib_sc_pANCOVA$coefficients["(Intercept)"]+size_rib_sc_pANCOVA$coefficients["ecology4tree"]
size_rib_sc_pANCOVA_chip_int_boot <- size_rib_sc_pANCOVA$bootstrap[,"(Intercept)"]
size_rib_sc_pANCOVA_fly_int_boot <- size_rib_sc_pANCOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANCOVA$bootstrap[,"ecology4flying"]
size_rib_sc_pANCOVA_ground_int_boot <- size_rib_sc_pANCOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANCOVA$bootstrap[,"ecology4ground"]
size_rib_sc_pANCOVA_tree_int_boot <- size_rib_sc_pANCOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANCOVA$bootstrap[,"ecology4tree"]

size_rib_sc_pANOVA_all_slope <- size_rib_sc_pANOVA$coefficients["lnsize"]
size_rib_sc_pANOVA_all_slope_boot <- size_rib_sc_pANOVA$bootstrap[,"lnsize"]

size_rib_sc_pANCOVA_chip_slope <- size_rib_sc_pANCOVA$coefficients["lnsize"]
size_rib_sc_pANCOVA_fly_slope <- size_rib_sc_pANCOVA$coefficients["lnsize"] + size_rib_sc_pANCOVA$coefficients["lnsize:ecology4flying"]
size_rib_sc_pANCOVA_ground_slope <- size_rib_sc_pANCOVA$coefficients["lnsize"] + size_rib_sc_pANCOVA$coefficients["lnsize:ecology4ground"]
size_rib_sc_pANCOVA_tree_slope <- size_rib_sc_pANCOVA$coefficients["lnsize"] + size_rib_sc_pANCOVA$coefficients["lnsize:ecology4tree"]
size_rib_sc_pANCOVA_chip_slope_boot <- size_rib_sc_pANCOVA$bootstrap[,"lnsize"]
size_rib_sc_pANCOVA_fly_slope_boot <- size_rib_sc_pANCOVA$bootstrap[,"lnsize"] + size_rib_sc_pANCOVA$bootstrap[,"lnsize:ecology4flying"]
size_rib_sc_pANCOVA_ground_slope_boot <- size_rib_sc_pANCOVA$bootstrap[,"lnsize"] + size_rib_sc_pANCOVA$bootstrap[,"lnsize:ecology4ground"]
size_rib_sc_pANCOVA_tree_slope_boot <- size_rib_sc_pANCOVA$bootstrap[,"lnsize"] + size_rib_sc_pANCOVA$bootstrap[,"lnsize:ecology4tree"]

size_rib_sc_pANCOVA_int_cof <- rbind(
  size_rib_sc_pANOVA_all_int,
  size_rib_sc_pANCOVA_chip_int,
  size_rib_sc_pANCOVA_fly_int,
  size_rib_sc_pANCOVA_ground_int,
  size_rib_sc_pANCOVA_tree_int)
size_rib_sc_pANCOVA_int_95 <- rbind(
  quantile(size_rib_sc_pANOVA_all_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_tree_int_boot, prob = c(0.025, 0.975)))
size_rib_sc_pANCOVA_slope_cof <- rbind(
  size_rib_sc_pANOVA_all_slope,
  size_rib_sc_pANCOVA_chip_slope,
  size_rib_sc_pANCOVA_fly_slope,
  size_rib_sc_pANCOVA_ground_slope,
  size_rib_sc_pANCOVA_tree_slope)
size_rib_sc_pANCOVA_slope_95 <- rbind(
  quantile(size_rib_sc_pANOVA_all_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_chip_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_fly_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_ground_slope_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANCOVA_tree_slope_boot, prob = c(0.025, 0.975)))

size_rib_sc_pANCOVA_cof <- data.frame(size_rib_sc_pANCOVA_int_cof, size_rib_sc_pANCOVA_int_95, size_rib_sc_pANCOVA_slope_cof, size_rib_sc_pANCOVA_slope_95 )
colnames(size_rib_sc_pANCOVA_cof) <- c("intercept", "Int_L95", "Int_U95", "Slope", "Slope_L95", "Slope_U95")
rownames(size_rib_sc_pANCOVA_cof) <- c("all", "chip", "flying", "ground", "tree")
round(size_rib_sc_pANCOVA_cof, digits = 2)

# plot lnrib_sc ~ lnsize*ecology
plot(lnrib_sc ~ lnsize, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5)
abline(a = size_rib_sc_pANCOVA_chip_int, b = size_rib_sc_pANCOVA_chip_slope, col = "#073b4c")
abline(a = size_rib_sc_pANCOVA_fly_int, b = size_rib_sc_pANCOVA_fly_slope, col = "#118ab2")
abline(a = size_rib_sc_pANCOVA_ground_int, b = size_rib_sc_pANCOVA_ground_slope, col = "#ffd166")
abline(a = size_rib_sc_pANCOVA_tree_int, b = size_rib_sc_pANCOVA_tree_slope, col = "#06d6a0", lty = 2)

### rib_sc ANOVA ###
size_rib_sc_pANOVA <- phylolm(lnrib_sc~lnsize+ecology4, phy = tree_pr, model = "lambda", boot = 1000)
summary(size_rib_sc_pANOVA)

# Make table of coefficients (intercept)
size_rib_sc_pANOVA_chip_int <- size_rib_sc_pANOVA$coefficients["(Intercept)"]
size_rib_sc_pANOVA_fly_int <- size_rib_sc_pANOVA$coefficients["(Intercept)"]+size_rib_sc_pANOVA$coefficients["ecology4flying"]
size_rib_sc_pANOVA_ground_int <- size_rib_sc_pANOVA$coefficients["(Intercept)"]+size_rib_sc_pANOVA$coefficients["ecology4ground"]
size_rib_sc_pANOVA_tree_int <- size_rib_sc_pANOVA$coefficients["(Intercept)"]+size_rib_sc_pANOVA$coefficients["ecology4tree"]
size_rib_sc_pANOVA_chip_int_boot <- size_rib_sc_pANOVA$bootstrap[,"(Intercept)"]
size_rib_sc_pANOVA_fly_int_boot <- size_rib_sc_pANOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANOVA$bootstrap[,"ecology4flying"]
size_rib_sc_pANOVA_ground_int_boot <- size_rib_sc_pANOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANOVA$bootstrap[,"ecology4ground"]
size_rib_sc_pANOVA_tree_int_boot <- size_rib_sc_pANOVA$bootstrap[,"(Intercept)"] + size_rib_sc_pANOVA$bootstrap[,"ecology4tree"]

size_rib_sc_pANOVA_int_cof <- rbind(
  size_rib_sc_pANOVA_chip_int,
  size_rib_sc_pANOVA_fly_int,
  size_rib_sc_pANOVA_ground_int,
  size_rib_sc_pANOVA_tree_int)
size_rib_sc_pANOVA_int_95 <- rbind(
  quantile(size_rib_sc_pANOVA_chip_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANOVA_fly_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANOVA_ground_int_boot, prob = c(0.025, 0.975)),
  quantile(size_rib_sc_pANOVA_tree_int_boot, prob = c(0.025, 0.975)))

size_rib_sc_pANOVA_cof <- data.frame(size_rib_sc_pANOVA_int_cof, size_rib_sc_pANOVA_int_95)
colnames(size_rib_sc_pANOVA_cof) <- c("intercept", "Int_L95", "Int_U95")
rownames(size_rib_sc_pANOVA_cof) <- c("chip", "flying", "ground", "tree")
round(size_rib_sc_pANOVA_cof, digits = 2)

# Plot rib_sc ~ ecotype
rib_sc_vp_ecology <- ggplot(data, aes(x=fct_relevel(ecology4, "chip", "flying", "ground", "tree"), y= lnrib_sc, fill = ecology4)) + 
  scale_y_continuous(name="ln rib_sc (mm)", limits = c(-.1,.15), breaks = scales::pretty_breaks(n = 8)) +
  xlab("") +
  geom_violin(trim=TRUE) + 
  scale_fill_manual(values= col_ecology4) + 
  geom_hline(yintercept=mean(lnrib_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") 
rib_sc_vp_ecology