### =============================
###
### Name: Appendix 2 - Code
###
### Project: The process of island abandonment by humans.
###
### Purpose: Analyses
###
### Author: Dr. Fabio Mologni
###
### Date: 2024-08-20
###
### =============================
###
### Notes:
### 
### =============================

###Woring directory

setwd("G:/My Drive/2. Work/Research/island abandonment/abandonment")

### Print more records

options(max.print=10000)

### =============================

### packages

library(ggplot2)
library(tidyverse)
library(ggpubr)
library(sf)
library(spData)


### =====================================
### ================ Map ================
### =====================================

ARC = read.csv("./map.csv")  #data

(archipelagos <- data.frame(longitude = c(-9.91992, -2.5, -22.96216, -1.26594, 9.4252151, 14.5071969, -18.8886986), 
                            latitude = c(53.26164, 49.5, 64.59799, 60.52965, 41.2240592, 44.7176644,
                                         35.4283604)))

archipelagos2 <- archipelagos %>%
  add_column(name = c('West Ireland', 'Channel Islands', 'West Iceland',
                      'Shetland', 'Maddalena Islands', 'Quarnero', 'Macaronesian Islands'))

ggplot(data = world) +
  geom_sf() +
  geom_point(data = archipelagos2, aes(x = longitude, y = latitude), size = 10, 
             shape = 16, colour = 'darkgreen') +
  geom_text(data = archipelagos2, aes(x = longitude, y = latitude, label=name), 
            hjust=0, vjust=-1, size=6)+
  coord_sf(xlim = c(-30, 30), ylim = c(30, 70), expand = FALSE)+
  guides(fill = "none") +
  labs (x= "Longitude", y = "Latitude") +
  theme_bw()+
  theme(plot.title = element_blank(),
        axis.title.y = element_text(size = 30, colour='black'),
        axis.title.x = element_text(size = 30, colour='black'),
        axis.text.y = element_text(size = 15, colour='black'),
        axis.text.x = element_text(size = 15, colour='black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = 'none') +
  theme(legend.title=element_blank())



### =====================================
### ======= Rates of abandonment ========
### =====================================

ABAN = read_csv("./abandoment_rates.csv")

ABAN$inhabitation <- factor(ABAN$inhabitation, levels = c("uninhabited", "abandoned", 
                                                          'inhabited'))

ABAN$inhabitation <- factor(ABAN$inhabitation,
                             levels = c("uninhabited", "abandoned", 
                                        'inhabited'),
                             labels = c("Consistently uninhabited", "Abandoned", 
                                        "Consistently inhabited"))

ggplot(ABAN, aes(y=islands_num, x=archipelago, fill=inhabitation, alpha=archipelago)) + 
  geom_bar(position="fill", stat="identity")+
  scale_fill_manual(values=c("lightgrey", "#C14444", "#8DB987"))+
  scale_alpha_manual(values=c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 1))+
  labs (x= "Archipelago", y = "Islands (%)") +
  theme_bw() +
  theme(plot.title = element_blank(),
        axis.title.y = element_text(size = 30),
        axis.title.x = element_text(size = 30),
        axis.text.y = element_text(size = 15),
        axis.text.x = element_text(size = 15),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.text = element_text(size = 19),
        legend.position = 'top')+
  theme(legend.title=element_blank()) +
  guides(fill = guide_legend(title = NULL),  # Ensure fill legend is displayed
    alpha = guide_none())  # Remove alpha legend




### ========================================
### =============== LOGISTIC ===============
### ========================================

ABAN_log = read_csv("./abandoment_logistic.csv")

###+++++++++++
###ALL ISLANDS
###+++++++++++

model.all = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_log)
summary(model.all)

ggplot(ABAN_log, aes(log(area), occupancy, linetype = time, group = time)) +
  geom_jitter(aes(shape=time, size=2), height = 0.05, width=0) +
  stat_smooth(method = "glm", method.args = list(family = "binomial"), se = F, colour = 'black', size=2) +
  scale_shape_manual(values=c(1, 16))+
  labs (x= "ln (Area)", y = "Island inhabitation (presence/absence)") +
  theme_bw() +
  theme(plot.title = element_blank(),
        axis.title.y = element_text(size = 30, colour='black'),
        axis.title.x = element_text(size = 30, colour='black'),
        axis.text.y = element_text(size = 20, colour='black'),
        axis.text.x = element_text(size = 20, colour='black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = 'none') +
  theme(legend.title=element_blank())


###INFLECTION_POINT_ALL
#exp = transform back in km2

#1850
ABAN_all_past = subset(ABAN_log, time=="past")
model.all.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_all_past)
summary(model.all.past)

coefficients.all.past <- coef(model.all.past)
exp(inflection_point.all.past <- -coefficients.all.past[1] / coefficients.all.past[2])

#current
ABAN_all_current = subset(ABAN_log, time=="current")
model.all.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_all_current)
summary(model.all.current)

coefficients.all.current <- coef(model.all.current)
exp(inflection_point.all.current <- -coefficients.all.current[1] / coefficients.all.current[2])


###+++++++++++
###SHETLAND+++
###+++++++++++

ABAN_SHE = subset(ABAN_log, archipelago=="SHE")

model.SHE = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_SHE)
summary(model.SHE)  

###INFLECTION_POINT_SHE
#1850
ABAN_SHE_past = subset(ABAN_SHE, time=="past")
model.SHE.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_SHE_past)
summary(model.SHE.past)

coeff.SHE.past <- coef(model.SHE.past)
exp(inflection_point.SHE.past <- -coeff.SHE.past[1] / coeff.SHE.past[2])

#current
ABAN_SHE_current = subset(ABAN_SHE, time=="current")
model.SHE.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_SHE_current)
summary(model.SHE.current)

coefficients.SHE.current <- coef(model.SHE.current)
exp(inflection_point.SHE.current <- -coefficients.SHE.current[1] / coefficients.SHE.current[2])


###+++++++++++
###CHANNEL++++
###+++++++++++

ABAN_CHA = subset(ABAN_log, archipelago=="CHA")

model.CHA = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_CHA)
summary(model.CHA)  

###INFLECTION_POINT_CHA
#1850
ABAN_CHA_past = subset(ABAN_CHA, time=="past")
model.CHA.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_CHA_past)
summary(model.CHA.past)

coeff.CHA.past <- coef(model.CHA.past)
exp(inflection_point.CHA.past <- -coeff.CHA.past[1] / coeff.CHA.past[2])

#current
ABAN_CHA_current = subset(ABAN_CHA, time=="current")
model.CHA.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_CHA_current)
summary(model.CHA.current)

coefficients.CHA.current <- coef(model.CHA.current)
exp(inflection_point.CHA.current <- -coefficients.CHA.current[1] / coefficients.CHA.current[2])



###++++++++++++++++++
###QUARNERO/KVARNER++
###++++++++++++++++++

ABAN_QUA = subset(ABAN_log, archipelago=="QUA")

model.QUA = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_QUA)
summary(model.QUA)  

###INFLECTION_POINT_MAD
#1850
ABAN_QUA_past = subset(ABAN_QUA, time=="past")
model.QUA.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_QUA_past)
summary(model.QUA.past)

coeff.QUA.past <- coef(model.QUA.past)
exp(inflection_point.QUA.past <- -coeff.QUA.past[1] / coeff.QUA.past[2])

#current
ABAN_QUA_current = subset(ABAN_QUA, time=="current")
model.QUA.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_QUA_current)
summary(model.QUA.current)

coefficients.QUA.current <- coef(model.QUA.current)
exp(inflection_point.QUA.current <- -coefficients.QUA.current[1] / coefficients.QUA.current[2])



###+++++++++++++
###MACARONESIA++
###+++++++++++++

ABAN_MAC = subset(ABAN_log, archipelago=="MAC")

model.MAC = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_MAC)
summary(model.MAC)  

###INFLECTION_POINT_MAC
#1850
ABAN_MAC_past = subset(ABAN_MAC, time=="past")
model.MAC.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_MAC_past)
summary(model.MAC.past)

coeff.MAC.past <- coef(model.MAC.past)
exp(inflection_point.MAC.past <- -coeff.MAC.past[1] / coeff.MAC.past[2])

#current
ABAN_MAC_current = subset(ABAN_MAC, time=="current")
model.MAC.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_MAC_current)
summary(model.MAC.current)

coefficients.MAC.current <- coef(model.MAC.current)
exp(inflection_point.MAC.current <- -coefficients.MAC.current[1] / coefficients.MAC.current[2])



###+++++++++++++
###ICELAND++++++
###+++++++++++++

ABAN_ICE = subset(ABAN_log, archipelago=="ICE")

model.ICE = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_ICE)
summary(model.ICE)  

###INFLECTION_POINT_ICE
#1850
ABAN_ICE_past = subset(ABAN_ICE, time=="past")
model.ICE.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_ICE_past)
summary(model.ICE.past)

coeff.ICE.past <- coef(model.ICE.past)
exp(inflection_point.ICE.past <- -coeff.ICE.past[1] / coeff.ICE.past[2])

#current
ABAN_ICE_current = subset(ABAN_ICE, time=="current")
model.ICE.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_ICE_current)
summary(model.ICE.current)

coefficients.ICE.current <- coef(model.ICE.current)
exp(inflection_point.ICE.current <- -coefficients.ICE.current[1] / coefficients.ICE.current[2])



###+++++++++++++
###IRELAND++++++
###+++++++++++++

ABAN_IRE = subset(ABAN_log, archipelago=="IRE")

model.IRE = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_IRE)
summary(model.IRE)  

###INFLECTION_POINT_IRE
#1850
ABAN_IRE_past = subset(ABAN_IRE, time=="past")
model.IRE.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_IRE_past)
summary(model.IRE.past)

coeff.IRE.past <- coef(model.IRE.past)
exp(inflection_point.IRE.past <- -coeff.IRE.past[1] / coeff.IRE.past[2])

#current
ABAN_IRE_current = subset(ABAN_IRE, time=="current")
model.IRE.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_IRE_current)
summary(model.IRE.current)

coefficients.IRE.current <- coef(model.IRE.current)
exp(inflection_point.IRE.current <- -coefficients.IRE.current[1] / coefficients.IRE.current[2])



###+++++++++++++
###ISTRIA+++++++
###+++++++++++++

ABAN_IST = subset(ABAN_log, archipelago=="IST")

model.IST = glm(occupancy ~ log(area)+time, family = binomial(link = "logit"), ABAN_IST)
summary(model.IST)  

###INFLECTION_POINT_IST
#1850
ABAN_IST_past = subset(ABAN_IST, time=="past")
model.IST.past = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_IST_past)
summary(model.IST.past)

coeff.IST.past <- coef(model.IST.past)
exp(inflection_point.IST.past <- -coeff.IST.past[1] / coeff.IST.past[2])

#current
ABAN_IST_current = subset(ABAN_IST, time=="current")
model.IST.current = glm(occupancy ~ log(area), family = binomial(link = "logit"), ABAN_IST_current)
summary(model.IST.current)

coefficients.IST.current <- coef(model.IST.current)
exp(inflection_point.IST.current <- -coefficients.IST.current[1] / coefficients.IST.current[2])






### ========================================
### ======= Abandonment vs Latitude ========
### ========================================

ABAN2 = read_csv("./abandoment_latitude.csv")

model1 = lm(aban_per ~ lat, ABAN2)
summary(model1)

ggplot(ABAN2, aes(x=lat, y=aban_per, label = archipelago)) +  
  geom_point(size=4, aes(colour=archipelago)) +
  geom_text(hjust=0.5, vjust=-0.8, size=6, aes(colour=archipelago))+
  stat_smooth(method=lm, lwd=2, se=TRUE, colour='black') +
  scale_colour_manual(values=c('#D81B60', '#1E88E5', '#FFC107', 
                               '#004D40', '#DCA4CE', '#522763', '#A4A95F'))+
  stat_regline_equation(label.y = 65, size=8, aes(label = after_stat(eq.label))) +
  stat_regline_equation(label.y = 61 , size=8, aes(label = ..rr.label..)) +
  xlab("Latitude (°)") + 
  ylab("Abandoned islands (%)") +
  theme_bw() +
  theme(plot.title = element_blank(),
        axis.title.y = element_text(size = 25, colour='black'),
        axis.title.x = element_text(size = 25, colour='black'),
        axis.text.y = element_text(size = 20, colour='black'),
        axis.text.x = element_text(size = 20, colour='black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = 'none') +
  theme(legend.title=element_blank())


### ===============================================
### ======= Ratio vs Latitude =====================
### ===============================================

hist(log(ABAN2$infl_ratio))

model1 = lm(log(infl_ratio) ~ lat, ABAN2)
summary(model1)

ggplot(ABAN2, aes(x=lat, y=log(infl_ratio), label = archipelago)) +  
  geom_point(size=4, aes(colour=archipelago)) +
  geom_text(hjust=0.5, vjust=-0.8, size=6, aes(colour=archipelago))+
  stat_smooth(method=lm, lwd=2, se=TRUE, colour='black', linetype='dashed') +
  scale_colour_manual(values=c('#D81B60', '#1E88E5', '#FFC107', 
                               '#004D40', '#DCA4CE', '#522763', '#A4A95F'))+
  stat_regline_equation(label.y = 5.5, size=8, aes(label = after_stat(eq.label))) +
  stat_regline_equation(label.y = 5, size=8, aes(label = ..rr.label..)) +
  xlab("Latitude (°)") + 
  ylab("log(inflection points ratio)") +
  theme_bw() +
  theme(plot.title = element_blank(),
        axis.title.y = element_text(size = 25, colour='black'),
        axis.title.x = element_text(size = 25, colour='black'),
        axis.text.y = element_text(size = 20, colour='black'),
        axis.text.x = element_text(size = 20, colour='black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = 'none') +
  theme(legend.title=element_blank())