library(ape)
library(readr)
library(tidyverse)
library(cowplot)
library(nlme)
library(MuMIn)

#Set WD - to wherever data is downloaded
setwd() 
####Import .csv####
RangeFilling <- read_csv("RangeFillingData.csv")

####Reorder Factors####
RangeFilling$Lecty<-factor(RangeFilling$Lecty, levels = c("Polylectic", "Oligolectic","No Lectic Status"))
RangeFilling$Overwintering<-factor(RangeFilling$Overwintering, levels = c("Prepupa", "Adult within cocoon","Adult (female only)"))

#Create tree - Visit https://tree.opentreeoflife.org/opentree/argus/ottol@419526/Apoidea
#Select "Download subtree as Newick string"
Tree<-read.tree("subtree-node-ott419526-Apoidea.tre")
Tree <- compute.brlen(Tree, power=1)

#Format tip.labels to match RangeFilling$Species format
Tree$tip.label <- sub("_ott.*", "", Tree$tip.label)
Tree$tip.label <- gsub("_", " ", Tree$tip.label)

# Subset Tree to include species in the RangeFilling Dataset
species_in_dataframe <- RangeFilling$Species
row_numbers_to_include <- c(11715, 11302, 11389, 7287, 819, 930, 734, 639, 947, 10506, 7975)
subsetted_data_frame <- Tree$tip.label[row_numbers_to_include, drop = FALSE]
matching_tips <- which(Tree$tip.label %in% species_in_dataframe| Tree$tip.label %in% subsetted_data_frame)
Tree <- drop.tip(Tree, tip = setdiff(Tree$tip.label, Tree$tip.label[matching_tips]))

#View Tree
plot(Tree)

#Run PGLS Regression analysis
model<-gls(LogPercent ~ Habitat.Breadth + Lecty + Overwintering + ITD, data=RangeFilling, correlation=corPagel(1, phy = Tree, form = ~Species))
summary(model)
AICc(model)

####Plot Model####
a<-ggplot(RangeFilling, aes(x = Lecty, y = LogPercent))+
  geom_boxplot()+
  theme_bw()+
  labs(x = "Lecty", y = "log(Climate envelope\n filled (%))")

b<-ggplot(RangeFilling, aes(x = Overwintering, y = LogPercent))+
  geom_boxplot()+
  theme_bw()+
  labs(x = "Overwintering Stage", y = "log(Climate envelope\n filled (%))")

c<-ggplot(RangeFilling, aes(x = ITD, y = LogPercent))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_bw()+
  labs(x = "Mean Intertegular Distance (mm)", y = "log(Climate envelope\n filled (%))")


plot_grid(a,b,c, ncol = 1)