plot_missing(d.raw)

plot_bar(d.raw)

plot_histogram(d.raw)

ARRD_outcome<-ggplot(data=d.raw, aes(y=ARRD, x=outcome, col=outcome,size=precision))+
                     geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+ # line for y = 0 (no difference between males and females)
                     geom_jitter(alpha=0.5)+ # display data points with a jitter, so they don't overlap
                     scale_size_continuous(range = c(1, 7))+ # change scaling of data points
                     labs(x="Type of outcome variable",y="ARRD") + # rename axis labels
                     theme_classic() # classic theme

ARRD_habitat<-ggplot(data=d.raw, aes(y=ARRD, x=habitat, col=habitat,size=precision))+
                     geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
                     geom_jitter(alpha=0.5)+ 
                     scale_size_continuous(range = c(1, 7))+
                     labs(x="Species habitat",y="ARRD") +
                     theme_classic()

ARRD_taxa<-ggplot(data=d.raw, aes(y=ARRD, x=taxonomic_group, col=taxonomic_group,size=precision))+
                  geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
                  geom_jitter(alpha=0.5)+ 
                  scale_size_continuous(range = c(1, 7))+
                  labs(x="Taxonomic group",y="ARRD",color = "taxon") +
                  theme_classic()

ARRD_origin<-ggplot(data=d.raw, aes(y=ARRD, x=origin, col=origin,size=precision))+
                    geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
                    geom_jitter(alpha=0.5)+ 
                    scale_size_continuous(range = c(1, 7))+
                    labs(x="Animals origin",y="ARRD")+
                    theme_classic()

plot_grid(ARRD_outcome, ARRD_habitat, ARRD_taxa, ARRD_origin, nrow = 2,ncol=2, labels = "auto", rel_heights = c(1, 1, 1,1))

ARRD_latitude <- ggplot(data = d.raw, aes(y = ARRD, x = abs_latitude, size = precision)) +
    geom_hline(yintercept = 0, linetype = 2, colour = "black", alpha = 0.5) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Absolute latitude",
    y = "ARRD") + theme_classic()

ARRD_longitude <- ggplot(data = d.raw, aes(y = ARRD, x = abs_longitude, size = precision)) +
    geom_hline(yintercept = 0, linetype = 2, colour = "black", alpha = 0.5) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Absolute Longitude",
    y = "ARRD") + theme_classic()

ARRD_elevation <- ggplot(data = d.raw, aes(y = ARRD, x = elevation, size = precision)) +
    geom_hline(yintercept = 0, linetype = 2, colour = "black", alpha = 0.5) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Elevation",
    y = "ARRD") + theme_classic()

plot_grid(ARRD_latitude, ARRD_longitude, ARRD_elevation, nrow = 2, labels = "auto",
    rel_heights = c(1, 1, 1, 1, 1))

ARRD_SDI <- ggplot(data = d.raw, aes(y = ARRD, x = SDI, size = precision)) + geom_hline(yintercept = 0,
    linetype = 2, colour = "black", alpha = 0.5) + geom_vline(xintercept = 0, linetype = 2,
    colour = "black", alpha = 0.5) + geom_point(col = "darkcyan", alpha = 0.5) +
    scale_size_continuous(range = c(1, 9)) + labs(x = "Size dimorphism index", y = "ARRD") +
    theme_classic()

ARRD_MDI <- ggplot(data = d.raw, aes(y = ARRD, x = MDI, size = precision)) + geom_hline(yintercept = 0,
    linetype = 2, colour = "black", alpha = 0.5) + geom_vline(xintercept = 0, linetype = 2,
    colour = "black", alpha = 0.5) + geom_point(col = "darkcyan", alpha = 0.5) +
    scale_size_continuous(range = c(1, 9)) + labs(x = "Mass dimorphism index", y = "ARRD") +
    theme_classic()

plot_grid(ARRD_SDI, ARRD_MDI, nrow = 1, ncol = 2, labels = "auto", rel_heights = c(1,
    1, 1, 1, 1))

ARRD_duration <- ggplot(data = d.raw, aes(y = ARRD, x = acc_duration, size = precision)) +
    geom_hline(yintercept = 0, linetype = 2, colour = "black", alpha = 0.5) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Duration of acclimation",
    y = "ARRD") + theme_classic()

ARRD_ramping <- ggplot(data = d.raw, aes(y = ARRD, x = ramping.deg.min, size = precision)) +
    geom_hline(yintercept = 0, linetype = 2, colour = "black", alpha = 0.5) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Ramping rate",
    y = "ARRD") + theme_classic()

plot_grid(ARRD_duration, ARRD_ramping, nrow = 1, ncol = 2, labels = "auto", rel_heights = c(1,
    1, 1, 1, 1))

ARRD_pub_year <- ggplot(data = d.raw, aes(y = ARRD, x = pub_year, col = pub_year,
    size = precision)) + geom_hline(yintercept = 0, linetype = 2, colour = "black",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + geom_point(col = "darkcyan",
    alpha = 0.5) + scale_size_continuous(range = c(1, 9)) + labs(x = "Publication year",
    y = "ARRD") + theme_classic()

ARRD_data_type <- ggplot(data = d.raw, aes(y = ARRD, x = data_type, col = data_type,
    size = precision)) + geom_hline(yintercept = 0, linetype = 2, colour = "black",
    alpha = 0.5) + scale_size_continuous(range = c(1, 7)) + geom_jitter(alpha = 0.5) +
    labs(x = "Data type", y = "ARRD") + theme_classic()

ARRD_peer <- ggplot(data = d.raw, aes(y = ARRD, x = peer_reviewed, col = peer_reviewed,
    size = precision)) + geom_hline(yintercept = 0, linetype = 2, colour = "black",
    alpha = 0.5) + scale_size_continuous(range = c(1, 7)) + geom_jitter(alpha = 0.5) +
    labs(x = "Peer-review", y = "ARRD") + theme_classic()


plot_grid(ARRD_pub_year, ARRD_data_type, ARRD_peer, nrow = 2, ncol = 2, labels = "auto",
    rel_heights = c(1, 1, 1, 1, 1))

# Distribution of the data points

WorldData <- map_data("world")

ggplot(d.raw, aes(x = longitude, y = latitude, size = precision, col = taxonomic_group_sub)) +
    geom_map(data = WorldData, map = WorldData, aes(x = long, y = lat, group = group,
        map_id = region), fill = "white", colour = "#7f7f7f", size = 0.5) + geom_point(alpha = 0.6) +
    scale_colour_viridis(discrete = T) + xlim(min(WorldData$long), max(WorldData$long)) +
    ylim(min(WorldData$lat), max(WorldData$lat)) + theme_minimal() + labs(x = "Longitude",
    y = "Latitude") + scale_size_continuous(range = c(2, 8))

taxa <- tnrs_match_names(names = levels(d$genus_species), context_name = "Animals")
kable(taxa) %>%
    kable_styling("striped", position = "left")  # Check if all species names match with the Open Tree Taxonomy

d$genus_species <- str_replace(d$genus_species, "Gammarus marinus", "Marinogammarus marinus")  # Replace names to match with Open Tree Taxonomy
d$genus_species <- str_replace(d$genus_species, "Hyale hirtipalma", "Apohyale hirtipalma")  # Replace names to match with Open Tree Taxonomy
d$genus_species <- str_replace(d$genus_species, "Uca panacea", "Leptuca panacea")  # Replace names to match with Open Tree Taxonomy

phylo_tree <- tol_induced_subtree(ott_ids = taxa$ott_id, label_format = "name")  # Generate the phylogenetic tree

is.binary(phylo_tree)  # Check if the tree is binary (= are there some polytomies?)

plot(phylo_tree)  # check for polytomies visually

set.seed(1234)
binary_tree <- multi2di(phylo_tree, random = TRUE)  # resolve polytomy at random with a set seed

is.binary(binary_tree)  # the tree is now binary

set.seed(1234)
binary_tree <- multi2di(phylo_tree, random = TRUE)  # This part need to be run once again to generate the same tree in a different code chunk

png("output/phylogeny.png", bg = "transparent", height = 15, width = 8, units = "in",
    res = 1200)
plot(binary_tree)
dev.off()

d$genus_species <- as.factor(d$genus_species)  # convert genus_species (species name) as a factor 
binary_tree$tip.label <- as.factor(binary_tree$tip.label)  # convert the species names in the tree as a factor

levels(binary_tree$tip.label) <- levels(d$genus_species)  # make sure names match
binary_tree$tip.label <- as.character(binary_tree$tip.label)  # convert names back to characters

# Compute branch length
phylo_branch <- compute.brlen(binary_tree, method = "Grafen", power = 1)

is.ultrametric(phylo_branch)  # check if tree is ultrametric

# Create a phylogenetic matrix

phylo_matrix <- vcv(phylo_branch, cor = T)  # the vcv function returns a variance-covariance matrix

kable(phylo_matrix, "html") %>%
    kable_styling("striped", position = "left") %>%
    scroll_box(width = "100%", height = "500px")