# Analysis of MalNAF version 1

# Load libraries

library(tidyverse)
library(readr)
library(magrittr)
library(geojsonio)
library(broom)
library(dplyr)
library(viridis)
library(RColorBrewer)
library(sf)
library(readr)
library(classInt)
library(ggpubr)
library(lcvplants)
library(LCVP)
library(hrbrthemes)
extrafont::loadfonts()
# Load in dataset

MalNAF <- read_csv("MalNAF_Naturalized.csv")
Native_range <- read_csv("MalNAF_Native_range.csv")
Malesia_environmental <- read_csv("Malesia_environmental_data.csv")
hab <- read_csv("MalNAF_SynHab.csv")

# total unique naturalized taxa and species in dataset after harmonization

MalNAF %$%
  n_distinct(Accepted_Taxon) # taxa

MalNAF %$%
  n_distinct(Accepted_binomial) # binomial

# dataset has multiple entries for some taxa due to multiple references
# simplify data into single entries to get p/a as binary

naturalize_taxa <-  MalNAF %>%
  arrange(region_id, Accepted_Taxon)%>%
  group_by(region_id, Accepted_Taxon) %>%
  slice_head()# 1 entry per species x region

head(naturalize_taxa)

# but to get a dataset for species only need to use accepted_binomial

naturalize_sp <-  MalNAF %>%
  arrange(region_id, Accepted_binomial)%>%
  group_by(region_id, Accepted_binomial) %>%
  slice_head() # 1 entry per species x region 

head(naturalize_sp)

# Get species richness per island group

naturalized_rich_sp <- naturalize_sp %>%
  group_by (region_id) %>%
  dplyr::count(region_id)%>%
  arrange(-n) %>% 
  rename (richness = n)

naturalized_rich_sp

summary(naturalized_rich_sp)
sd(naturalized_rich_sp$richness)

# add relative and perc naturalized to environmental data

Malesia_env <- Malesia_environmental %>%
  mutate(Relative_naturalized_richness = Naturalize/Native_ric) %>%
  mutate(Perc_total_flora_naturalized = 
           (Naturalize/(Naturalize + Native_ric))*100)

# get a species list for the naturalized taxa

mal_sp <- naturalize_sp %>%
  group_by(Accepted_binomial) %>%
  slice_head() %>%
  select(Accepted_Family, Accepted_Taxon, Accepted_binomial) 

## How frequent are naturalized species across the island groups?

naturalized_freq <- naturalize_sp %>%
  group_by (Accepted_binomial) %>%
  dplyr::count(Accepted_binomial) %>%
  arrange(-n) %>% 
  rename (frequency = n)

naturalized_freq

summary(naturalized_freq)

sd(naturalized_freq$frequency)

# which species are common to every biogeographic unit?

freq_allBU <- naturalized_freq %>%
  filter(frequency == "10")

# add family name to table

most_freq_sp_fam <- freq_allBU %>%
  left_join(naturalize_sp, by = "Accepted_binomial") %>%
  group_by(Accepted_binomial) %>% 
  slice_head() %>%
  select (Accepted_binomial, Accepted_Taxon, Accepted_Family)

# how many families represented by the most frequent species?

most_freq_sp_fam %$%
  n_distinct(Accepted_Family)

# which families have the most widespread species?

top_fam_freq <- most_freq_sp_fam %>%
  group_by(Accepted_Family) %>%
  mutate (count = n()) %>%
  select(Accepted_Family, count) %>%
  group_by (Accepted_Family) %>%
  slice_head () %>%
  arrange(-count) %>%
  rename(No_sp = count)

top_fam_freq

# plot frequency

p_freq <- naturalized_freq %>%
  ggplot(aes(x = frequency)) + 
  geom_histogram(color="darkgray", fill="darkgray") +
  # geom_vline(aes(xintercept=mean(frequency)),
  # color="black", linetype="dashed", size=1) +
  xlab("No. Island Groups") + 
  ylab("No. Species") +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5))

p_freq

# how many are widespread and occur in 5 island groups or more

widespread <- naturalized_freq %>%
  filter(frequency >= 5) # 5 or over

widespread %$%
  n_distinct(Accepted_binomial) 

rare <- naturalized_freq %>%
  filter(frequency == 1) # only 1

rare %$%
  n_distinct(rare)
  
## How many families and genera in the dataset?

# first one will produce a table of all families 

mal_fam <- naturalize_sp %>%
  group_by(Accepted_Family) %>%
  slice_head() %>%
  select(Accepted_Family) %>%
  arrange(Accepted_Family)

head(mal_fam)

# all genera by families

mal_gen <- naturalize_sp %>%
  group_by(Accepted_genus) %>%
  slice_head() %>%
  select(Accepted_Family, Accepted_genus) %>%
  arrange(Accepted_Family)

# all species in each families

mal_sp <-  naturalize_sp %>%
  group_by(Accepted_binomial) %>%
  slice_head() %>%
  select(Accepted_Family, Accepted_genus, Accepted_Taxon, Accepted_binomial) %>%
  arrange(Accepted_Family, Accepted_genus, Accepted_Taxon)

# count no genus in family and gen in family and make new table

mal_gen_count <- mal_gen %>%
  group_by (Accepted_Family) %>% 
  mutate (count = n()) %>%
  select(Accepted_Family, count) %>%
  group_by (Accepted_Family) %>%
  slice_head () %>%
  arrange(-count) %>%
  rename(No_genus = count)

# count no. species in family 

mal_sp_count <- mal_sp %>%
  group_by (Accepted_Family) %>% 
  mutate (count = n()) %>%
  select(Accepted_Family, count) %>%
  group_by (Accepted_Family) %>%
  slice_head () %>%
  arrange(-count) %>%
  rename(No_species = count)

## combine the two together

mal_fam_count <- mal_gen_count %>%
  left_join(mal_sp_count) %>%
  arrange(-No_genus, -No_species)

head(mal_fam_count)

# average 

mal_fam_count %>%
  summary()

## Species richness plots

#import spatial data 

mal <- st_read(
  "Malesia_environmental_data.shp")

# check data

plot(mal)

# need to combine the .shp file with species richness vars.

mal_env_filt <- Malesia_env %>%
  select("Island_nam", "Relative_naturalized_richness", "Perc_total_flora_naturalized")

mal_rich <- mal %>%
  left_join (mal_env_filt, by = ("Island_nam"))

# fig 3
p1 <- mal_rich %>%
  ggplot() +
  geom_sf(aes(fill= Native_ric))+
  scale_fill_continuous(
    low = "#fed8b1", high = "#820000")+
  labs(x='Longitude',y='Latitude') +
  theme_bw()+
  #theme(legend.position = 'bottom')+
  labs(fill = "Native richness")

p1

p2 <- mal_rich %>%
  ggplot() +
  geom_sf(aes(fill= Naturalize))+
  scale_fill_continuous(
    low = "#fed8b1", high = "#820000")+
  labs(x='Longitude',y='Latitude') +
  theme_bw()+
  #theme(legend.position = 'bottom')+
  labs(fill = "Naturalized richness")

p2 # raw richness



p3 <- mal_rich %>%
  ggplot() +
  geom_sf(aes(fill= Relative_naturalized_richness))+
  scale_fill_continuous(
    low = "#fed8b1", high = "#820000")+
  labs(x='Longitude',y='Latitude') +
  theme_bw()+
  #theme(legend.position = 'bottom')+
  labs(fill = "Relative richness") # not sure how to caption

p3

ggarrange(p1,p2,p3, nrow=3, align = "v")

## Family Naturalization index

lcvp_fam <- lcvp_group_search(mal_fam$Accepted_Family, 
                              search_by = "Family") # search for fam in LCVP

lcvp_fam_count <- lcvp_fam %>%
  filter(Status == "accepted") %>%
  group_by(Family) %>%
  summarise(fam_species_global = n_distinct(Output.Taxon)) # get spp. counts

# count no. species in family 

mal_fam_bor <- naturalize_sp %>%
  filter(region_id == "Borneo")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Borneo = n_distinct(Accepted_binomial)
  )


mal_fam_pen <- naturalize_sp %>%
  filter(region_id == "Peninsular Malaysia")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Peninsular_Malaysia = n_distinct(Accepted_binomial)
  )

mal_fam_sgp <- naturalize_sp %>%
  filter(region_id == "Singapore")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Singapore = n_distinct(Accepted_binomial)
  )

mal_fam_sum <- naturalize_sp %>%
  filter(region_id == "Sumatra")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Sumatra = n_distinct(Accepted_binomial)
  )

mal_fam_jav <- naturalize_sp %>%
  filter(region_id == "Java")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Java = n_distinct(Accepted_binomial)
  )

mal_fam_LSI <- naturalize_sp %>%
  filter(region_id == "LSI")%>%
  group_by(Accepted_Family) %>%
  summarise(
    LSI = n_distinct(Accepted_binomial)
  )

mal_fam_sul <- naturalize_sp %>%
  filter(region_id == "Sulawesi")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Sulawesi = n_distinct(Accepted_binomial)
  )

mal_fam_maluku <- naturalize_sp %>%
  filter(region_id == "Maluku")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Maluku = n_distinct(Accepted_binomial)
  )

mal_fam_phi <- naturalize_sp %>%
  filter(region_id == "Philippines")%>%
  group_by(Accepted_Family) %>%
  summarise(
    Philippines = n_distinct(Accepted_binomial)
  )

mal_fam_nwg <- naturalize_sp %>%
  filter(region_id == "New Guinea")%>%
  group_by(Accepted_Family) %>%
  summarise(
    NWG = n_distinct(Accepted_binomial)
  )

# now combine

fam_island <- mal_fam_bor %>%
  full_join(mal_fam_pen) %>%
  full_join(mal_fam_sgp) %>%
  full_join(mal_fam_sum) %>%
  full_join(mal_fam_jav) %>%
  full_join(mal_fam_LSI) %>%
  full_join(mal_fam_sul) %>%
  full_join(mal_fam_phi) %>%
  full_join(mal_fam_maluku) %>%
  full_join(mal_fam_nwg) %>%
  rename(Family = Accepted_Family) %>%
  replace(is.na(.), 0)


# now to calculate 


# need to pivot longer then add species richness per island group 

FI_long <- fam_island %>%
  pivot_longer(cols = -Family, names_to = "Island", values_to = "species")

# finally add back in the totals for world and Malesia to allow for calc
# also need total richness per island

family_index <- FI_long %>%
  mutate(Island = replace(Island, Island == "Peninsular_Malaysia", "Peninsular Malaysia")) %>%
  mutate(Island = replace(Island, Island == "NWG", "New Guinea")) %>%
  left_join(lcvp_fam_count) %>%
  left_join(naturalized_rich_sp, by = c("Island" = "region_id")) %>%
  mutate(FI = (species/(fam_species_global * richness))*100) %>% # calculate FI 
  mutate(Island = replace(Island, Island == "LSI", "Lesser Sunda Islands")) %>%
  mutate(F_ratio = (species/fam_species_global)*100)

# only include families which have at least 10 species naturalized in an individual island 
# as these are important in the flora and cannot plot all families

family_10 <- family_index %>%
  filter(species > 9) %>% # above 10 in a region
  group_by(Family) %>%
  slice_head() %>%
  select(Family)

family_index_10 <- family_10 %>%
  left_join(family_index)

family_index_10 %$%
  n_distinct(Family)

FI_1 <- family_index_10 %>%
  ggplot(aes(y= FI, x= fct_reorder(Family, FI))) +
  geom_boxplot(outlier.shape = NA) +
  scale_fill_viridis(discrete = TRUE, alpha=0.6) +
  geom_jitter(aes(color= Island), 
              size=2.5, alpha=0.9) +
  theme_ipsum(axis_text_size = 14, hrbrthemes.loadfonts = TRUE) +
  theme(
    legend.title = element_text(size=14, color = "dimgrey"), 
    legend.text = element_text(size=14, color = "dimgrey"),
    legend.position="right"
  ) +
  xlab("")+
  ylab("") +
  coord_flip()

FI_1 <- FI_1 + labs(fill = "Island/Island group")

FI_1

FI_1 <- print(FI_1 + labs(colour = "Island/Island group"))

FI_2 <- family_index_10 %>%
  ggplot(aes(y= FI, x= fct_reorder(Island, FI))) +
  geom_boxplot(outlier.shape = NA) +
  scale_fill_viridis(discrete = TRUE, alpha=0.6) +
  geom_jitter(aes(color= Family), 
              size=2.5, alpha=0.9) +
  theme_ipsum(axis_text_size = 14) +
  theme(
    legend.title = element_text(size=14, color = "dimgrey"), 
    legend.text = element_text(size=14, color = "dimgrey"),
  ) +
  xlab("")+
  ylab("") +
  coord_flip()

ggarrange(FI_1, FI_2, nrow = 2, align = "hv", legend = "right") # fig.4


## Native Range of the naturalized taxa

range <- Native_range %>%
  select (binomial, Europe, Africa, Asia.Temperate, 
          Asia.Tropical, Australasia, Pacific, Northern.America, 
          Southern.America, Antarctica, Hybrid, Unknown)

head(range)

# need it in long format for analysis using dplyr

range_long <- range %>%
  pivot_longer(cols = -binomial, names_to = "Continent", values_to = "Present") %>%
  filter(Present >= 1) %>%
  select(binomial, Continent)

head(range_long)

# which continent has donated the most naturalized plant species?

range_rich <- range_long %>%
  group_by (Continent) %>%
  dplyr::count(Continent)%>%
  arrange(-n) %>% 
  rename (richness = n)
#  knitr::kable()

range_rich # all good!

## plot it

p_range <- range_long %>%
  group_by (Continent) %>%
  dplyr::count(Continent)%>%
  ggplot(aes(
    x = n,
    y = reorder(Continent, n))) +
  geom_bar(stat = "identity",
           colour = "darkGray", fill = "darkGray", 
           width = 0.5) +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  xlab("Naturalized Species") + 
  ylab("Native Range")


p_range

range_rich %>%
  summary()

sd(range_rich$richness)

range_rich <- range_rich %>%
  mutate (perc_origin = round(((richness/3494) * 100), 2)) 

## need a new table combining the taxon x region data with taxon x native origin data 

mal_range <- naturalize_sp %>% 
  left_join(range_long, by = c("Accepted_binomial" = "binomial")) %>%
  select(Accepted_binomial, region_id, Continent)

mal_range  

mal_r_1 <- mal_range %>% 
  group_by(region_id) %>%
  count(Continent) %>%
  rename (no_spp = n)

mal_r_2 <- mal_r_1 %>%
  group_by(region_id) %>%
  mutate (total_orig = sum(no_spp))

mal_r_2

cont_origin <- mal_r_2 %>% 
  mutate (perc_origin = (no_spp / total_orig) * 100)

p_orig_stack <- cont_origin %>%
  mutate(Continent = replace(Continent, Continent == "Asia.Tropical", "Asia Tropical")) %>%
  mutate(Continent = replace(Continent, Continent == "Asia.Temperate", "Asia Temperate")) %>%
  mutate(Continent = replace(Continent, Continent == "Southern.America", "Southern America")) %>%
  mutate(Continent = replace(Continent, Continent == "Northern.America", "Northern America")) %>%
  mutate(region_id = replace(region_id, region_id == "LSI", "Lesser Sunda Islands")) %>%
  mutate(Continent = fct_relevel(Continent, "Unknown","Hybrid","Antarctica",
                                 "Pacific","Australasia","Asia Tropical", 
                                 "Asia Temperate", "Europe","Africa",
                                 "Southern America", "Northern America")) %>%
  ggplot(aes(fill = Continent, y= perc_origin, 
             x= fct_relevel(region_id, "New Guinea", "Maluku", "Philippines",
                            "Sulawesi","Lesser Sunda Islands", "Java","Borneo","Singapore", 
                            "Peninsular Malaysia","Sumatra"))) +                                                   
  geom_bar(position="stack", stat="identity") +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  scale_fill_brewer(palette = "Spectral") +
  ylab("% Total naturalized species") + 
  xlab("Islands/Island groups") +
  guides(fill=guide_legend(title="Native range", reverse = TRUE)) +
  coord_flip()

p_orig_stack # figure 5

## Habitat Data

# how many species in data?

hab %$%
  n_distinct(binomial)

# how many species in each country?

sp_no <- hab %>%
  group_by(Region_ID) %>%
  summarise(no_species = n_distinct(binomial))
  
sp_no

# how many records by country total

sp_records <- hab %>%
  group_by(Region_ID) %>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records)

sp_records

# What proportion of naturalized species occur in each habitat type across Malesia?

all_hab <- hab %>%
  group_by(binomial, Habitat) %>%
  slice_head()%>% # for each species no double habitat entries
  ungroup() %>%
  group_by(Habitat)%>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records)

all_hab

# What proportion of records reported from each habitat type across Malesia?

all_hab2 <- hab %>%
  group_by(binomial, Habitat) %>%
  ungroup() %>%
  group_by(Habitat)%>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records)

all_hab2

# plot

Malesia_hab <- hab %>%
  group_by(binomial, Habitat) %>%
  slice_head()%>% # for each species no double habitat entries
  ungroup() %>%
  group_by(Habitat)%>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records) %>%
  mutate(Habitat = replace(Habitat, Habitat == "1", "Forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "2", "Open forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "3", "Scrub")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4", "Grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4a", "Natural grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4b", "Anthropogenic grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "5", "Sandy")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "6", "Rocky")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "7", "Dryland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "8", "Saline")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "9", "Riparian")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "10", "Wetland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "11", "Aquatic")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12", "Anthropogenic")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12a", "Ruderal/Urban")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12b", "Agricultural")) %>%
  mutate(Region = "Malesia")

# set up colour

## You need to expand palette size

cols <- 15
mycolors <- colorRampPalette(brewer.pal(11, "Spectral"))(cols)
mycolours <- rev(mycolors)

# plot

mal_hab_stack <- Malesia_hab %>%
  ggplot(aes(fill = fct_relevel(Habitat, "Forest", "Open forest", "Scrub", "Grassland", 
                                "Natural grassland", "Anthropogenic grassland",
                                "Sandy", "Rocky", "Dryland", "Saline", 
                                "Riparian", "Wetland", "Aquatic", 
                                "Ruderal/Urban","Agricultural"),
             y= perc_hab_records, x = Region)) +                                                   
  geom_bar(position="stack", stat="identity") +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  scale_fill_manual(values = mycolours) +
  guides(fill=guide_legend(title="Habitat type", reverse = TRUE)) +
  ylab("") + 
  xlab("Region") + 
  coord_flip()

mal_hab_stack # use as fig. 6a


# What is the proportion by country?

country_hab <- hab %>%
  group_by(Region_ID, Habitat)%>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records) %>%
  mutate(Region_ID = replace(Region_ID, Region_ID == "PH", "Philippines")) %>%
  mutate(Region_ID = replace(Region_ID, Region_ID == "ID", "Indonesia")) %>%
  mutate(Region_ID = replace(Region_ID, Region_ID == "TL", "Timor Leste")) %>%
  mutate(Region_ID = replace(Region_ID, Region_ID == "PG", "Papua New Guinea")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "1", "Forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "2", "Open forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "3", "Scrub")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4", "Grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4a", "Natural grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4b", "Anthropogenic grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "5", "Sandy")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "6", "Rocky")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "7", "Dryland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "8", "Saline")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "9", "Riparian")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "10", "Wetland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "11", "Aquatic")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12", "Anthropogenic")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12a", "Ruderal/Urban")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12b", "Agricultural"))

# set up colour

## You need to expand palette size

cols <- 15
mycolors <- colorRampPalette(brewer.pal(11, "Spectral"))(cols)
mycolours <- rev(mycolors)

# plot

p_hab_stack <- country_hab %>%
  ggplot(aes(fill = fct_relevel(Habitat, "Forest", "Open forest", "Scrub", "Grassland", 
                                "Natural grassland", "Anthropogenic grassland",
                                "Sandy", "Rocky", "Dryland", "Saline", 
                                "Riparian", "Wetland", "Aquatic", 
                                "Ruderal/Urban","Agricultural"),
             y= perc_hab_records, 
             x= fct_relevel(Region_ID, "Papua New Guinea", "Philippines",
                            "Timor Leste", "Indonesia"))) +                                                   
  geom_bar(position="stack", stat="identity") +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  scale_fill_manual(values = mycolours) +
  guides(fill=guide_legend(title="Habitat type", reverse = TRUE)) +
  ylab("") + 
  xlab("Country") +
  coord_flip()

p_hab_stack # use as fig. 6b

# what proportion do widespread species occur in vs. non-widespread

widespread_sp <- widespread %>%
  filter (frequency > 4) %>%
  select (Accepted_binomial) 

habitat_filt <- hab %>%
  group_by(binomial, Habitat) %>%
  slice_head()%>% # for each species no double habitat entries
  ungroup() 

hab_wide <- habitat_filt %>%
  semi_join(widespread_sp, by = c("binomial" = "Accepted_binomial")) %>%
  mutate(Distribution = "Widespread")

hab_rare <- habitat_filt %>%
  anti_join(widespread_sp, by = c("binomial" = "Accepted_binomial")) %>%
  mutate(Distribution = "Rare")

hab_frequency <- hab_wide %>%
  full_join(hab_rare)

freq_hab <- hab_frequency %>%
  ungroup() %>%
  group_by(Distribution, Habitat)%>%
  summarise(hab_records = (count = n())) %>%
  mutate(perc_hab_records = (hab_records/sum(hab_records))*100) %>%
  arrange(-perc_hab_records) %>%
  mutate(Habitat = replace(Habitat, Habitat == "1", "Forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "2", "Open forest")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "3", "Scrub")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4", "Grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4a", "Natural grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "4b", "Anthropogenic grassland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "5", "Sandy")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "6", "Rocky")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "7", "Dryland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "8", "Saline")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "9", "Riparian")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "10", "Wetland")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "11", "Aquatic")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12a", "Ruderal/Urban")) %>%
  mutate(Habitat = replace(Habitat, Habitat == "12b", "Agricultural"))

## You need to expand palette size

cols <- 15
mycolors <- colorRampPalette(brewer.pal(11, "Spectral"))(cols)
mycolours <- rev(mycolors)

# plot

p_hab_stack_freq <- freq_hab %>%
  ggplot(aes(fill = fct_relevel(Habitat, "Forest", "Open forest", "Scrub", "Grassland", 
                                "Natural grassland", "Anthropogenic grassland",
                                "Sandy", "Rocky", "Dryland", "Saline", 
                                "Riparian", "Wetland", "Aquatic", 
                                "Ruderal/Urban","Agricultural"),
             y= perc_hab_records, 
             x= Distribution)) +                                                   
  geom_bar(position="stack", stat="identity") +
  theme (text = element_text(size = 20), 
         # Hide panel borders and remove grid lines
         panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         panel.background = element_blank(),
         # Change axis line
         axis.line = element_line(colour = "black"),
         axis.title = element_text(size = 15)) +
  scale_fill_manual(values = mycolours) +
  guides(fill=guide_legend(title="Habitat type", reverse = TRUE)) +
  ylab("% Species") + 
  xlab("Distribution") +
  coord_flip()

p_hab_stack_freq # use as fig 6c.

ggarrange(mal_hab_stack, p_hab_stack , p_hab_stack_freq, nrow = 3, align = "hv",  
          common.legend = TRUE, legend = "right")

## Multi-Linear Regression 

# libraries

library(ggbiplot)    # plotting PCA
library(tidymodels)
library(psych)

### Step 1 PCA of bioclim 

mal_bioclim <- mal_rich %>%
  as.data.frame() %>% # to drop geometry
  select(Island, bio1_mean, bio2_mean, bio3_mean, bio4_mean, bio5_mean, bio6_mean,
         bio7_mean, bio8_mean, bio9_mean, bio10_mean, bio11_mean, bio12_mean, 
         bio13_mean, bio14_mean, bio15_mean, bio16_mean, bio17_mean, bio18_mean,
         bio19_mean) # select bioclim var

mal_bioclim.pca <- mal_bioclim %>%
  select(-Island) %>%
  prcomp(center = TRUE,scale. = TRUE)

mal_bioclim.pca %>%
  summary()

mal_bioclim.pca %>%
  str ()

mal_bioclim.pca %>%
  ggbiplot(labels= mal_bioclim$Island)


bioclim_vars <- as.data.frame(mal_bioclim.pca$rotation) # extract loadings


bioclim_scores <- as.data.frame(mal_bioclim.pca$x)

bioclim_scores <- bioclim_scores %>%
  add_column(Island = mal_bioclim$Island) # add labels back in

# add scores back into the env var data frame for regression tree analysis and
# subset only var needed for the analysis

mal_rich_filt <- mal_rich %>%
  left_join(bioclim_scores) %>%
  select(Island, Phyto, x_coord, y_coord, area, alt_range, 
         perc_urb, perc_agri, pop_dens, GDP_sum, Naturalize, 
         PC1, PC2, PC3) %>%
  mutate(log10_perc_urb = log10(perc_urb)) %>%
  mutate(log10_pop_dens = log10(pop_dens)) %>%
  mutate(sqrt_GDP_sum = sqrt(GDP_sum))%>%
  mutate(sqrt_Area = sqrt(area))%>%
  st_drop_geometry()

# Normality 
# Check out normality of the distributions

# First taking a look at dependent variables.

mal_rich_filt %>% 
  pull(Naturalize) %>% stats::shapiro.test() # normal

# Geographic drivers

# Area 

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = area
    )
  ) +
  ggplot2::geom_histogram(bins=5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(area) %>% stats::shapiro.test() # not normal, skewed left with long right tail

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = sqrt_Area
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()

mal_rich_filt %>%
  pull(sqrt_Area) %>% stats::shapiro.test() # now normal

# Altitude 


mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = alt_range
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()

mal_rich_filt %>% 
  pull(alt_range) %>% stats::shapiro.test() 

# Climate 

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = PC1
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(PC1) %>% stats::shapiro.test()

#pc2

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = PC2
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(PC2) %>% stats::shapiro.test()


# pc3

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = PC3
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(PC3) %>% stats::shapiro.test()

# Anthropogenic

# Urban area

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = perc_urb
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(perc_urb) %>% stats::shapiro.test() # not normal 

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = log10_perc_urb
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()

mal_rich_filt %>% 
  pull(log10_perc_urb) %>% stats::shapiro.test() # now normal

# Agricultural area

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = perc_agri
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(perc_agri) %>% stats::shapiro.test()  

# Population density

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = pop_dens
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(pop_dens) %>% stats::shapiro.test()  # not normal

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = log10_pop_dens
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>% 
  pull(log10_pop_dens) %>% stats::shapiro.test()  #  normal

# GDP

mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = GDP_sum
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>%
  ggplot2::ggplot(
    aes(
      x = sqrt_GDP_sum
    )
  ) +
  ggplot2::geom_histogram(bins = 5) +
  ggplot2::theme_bw()


mal_rich_filt %>%
  pull(GDP_sum) %>% stats::shapiro.test() # not normal

mal_rich_filt %>%
  pull(sqrt_GDP_sum) %>% stats::shapiro.test() # normal

## Correlation tests 

# Use Spearmans because small dataset even though all variables are 
# normally distributed after transformation

mal_rich_filt %>%
  select(
    Naturalize,
    PC1, PC2, PC3, alt_range, log10_perc_urb, perc_agri,
    log10_pop_dens, sqrt_GDP_sum,sqrt_Area
  ) %>% 
  pairs.panels(method = "spearman")

# MLR

library(lmtest)
library(lm.beta)
library(reshape2)
library(car)

# Model 1

mal_rich_mod1 <- lm(Naturalize ~ log10_pop_dens + PC3, data = mal_rich_filt)

mal_rich_mod1 %>%
  summary()

lm.beta(mal_rich_mod1)

# Test the model meets all the assumptions

#1) Normality (Shapiro Wilke)

mal_rich_mod1 %>% 
  stats::rstandard() %>% 
  stats::shapiro.test()

#Test is not significant so the residuals are normally distributed. The model is robust.

#2) Homoscedasticity (Breusch-Pagan)

mal_rich_mod1 %>% 
  lmtest::bptest()

# Not significant so the model is robust. 

# 3) Independance (Durbin-Watson)

mal_rich_mod1 %>%
  lmtest::dwtest()

# 4) multicolinnearity

mal_rich_mod1 %>%
  vif()

#other checks

mal_rich_mod1 %>%
  plot(which = c(1))

mal_rich_mod1 %>%
  plot(which = c(2))

# cooks distance for outliers

mal_rich_mod1 %>%
  plot(which = c(5))

# Model is significant and robust meeting all assumptions. 

## Influence of outliers on important drivers

mal_rich2 <- mal_rich_filt %>%
  filter(Island != "SGP") %>%
  filter(Island != "NWG")

mal_rich_filt %$%
  cor.test(Naturalize, perc_agri, method = "spearman")

mal_rich2 %$%
  cor.test(Naturalize, perc_agri, method = "spearman")