g=read.csv("island.csv",header=TRUE, na.strings=c("", "NA"),  sep = ";")

library(ape)
library(devtools)
library(ggplot2)
library(phyr)
library("gridExtra")
library("cowplot")

tree_vig <- read.nexus("island.nex") 
cons_vig = ape::consensus(tree_vig, p=0.5) # To get the majority-rule consensus tree, use p = 0.5


cons_vig2 = compute.brlen(cons_vig) # computes branch lengths, needed for next steps
vig_phylo <- vcv.phylo(cons_vig2, cor = T) 

g = data.frame(g)
island <- pglmm(logitsurvival ~ logbodymass + logclutchsize + hemisphere + status + abslat + hemisphere*abslat + breeding + family + (1 | species), data = g, 
     cov_ranef = list(sp = vig_phylo), family = "gaussian", REML = TRUE)
summary(island)
plot(island$H)

output=cbind(g, residuals=residuals(island), predicted=fitted(island))
output$hemisphere <- factor(output$hemisphere, levels = c("north", "south"), 
                  labels = c("Northern", "Southern"))
output$calculated=exp(output$predicted)/(1+(exp(output$predicted))
fig=ggplot(data=output, aes(x=abslat, y=calculated, color=status)) + 
    labs(x="Absolute latitude", y="Fitted annual adult apparent survival")+ 
    geom_smooth(method=lm, se = FALSE) + geom_point(size=1)+ theme_classic() +
                         scale_color_discrete(name="Location",
                         breaks=c("C", "I"), labels=c("Continent", "Island"), guide = guide_legend(reverse=TRUE))
figure=fig+facet_grid(. ~ hemisphere, scales="free")
figure
ggsave("survival.tiff", plot=figure, width = 17, height = 12, units = "cm", dpi = 300)
dev.off()