# Supplementary R script for manuscript: A global comparison of stream diatom beta diversity on islands vs. continents across scales.

# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
# SESSION PREPARATION
# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .

## Load libraries

library(readxl)
library(betapart)
library(BAT)
library(rcompanion)
library(geodist)
library(psych)
library(minpack.lm)
library(vegan)
library(MuMIn)
library(lme4)
library(RLRsim)


# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
# DATA OF SELECTED SAMPLED SITES FOR ISLANDS AND CAEs 
# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .

### Read table with spatial coordinates, environmental and presence/absence data ###

## For the analyses at GENUS LEVEL

# Read the American region data
american_islands <- read_excel("Islands_American_region-genus.xlsx")
american_caes <-  read_excel("CAEs_American_region-genus.xlsx")

# Read the European region data
european_islands <- read_excel("Islands_European_region-genus.xlsx")
european_caes <-  read_excel("CAEs_European_region-genus.xlsx")

# Read the African region data
african_islands <- read_excel("Islands_African_region-genus.xlsx")
african_caes <-  read_excel("CAEs_African_region-genus.xlsx")

# Read the Pacific region data
pacific_islands <- read_excel("Islands_Pacific_region-genus.xlsx")


## Islands Characteristics
# Read the island characteristics table
island_char <- read_excel("Islands_characteristics.xlsx")



# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
# FUNCTIONS CREATED TO PERFORM THE SMALL SCALE ANALYSES 
# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
#######


### Function for compute the dissimilarity within CAEs and the comparison with dissimilarities within islands ###

beta_small_scale <- function(isl,isl_id,cae_pa,cae_id){ 
  
  cae_sor <- c(); cae_sim <- c(); cae_sne <- c()
  cae_total <- c(); cae_repl <- c(); cae_rich <- c() 
  
  for (i in 1:length(isl_id)) {
    
    CAE_selec <- cae_pa[which(cae_id$CAE_code %in% isl_id[i]),] 
    
    cae_beta_bas <- beta.pair(CAE_selec)
    vals0 <- as.vector(cae_beta_bas$beta.sor); vals1 <-as.vector(cae_beta_bas$beta.sim) ; vals2<-as.vector(cae_beta_bas$beta.sne)
    cae_sor <- c(cae_sor,vals0); cae_sim<- c(cae_sim,vals1); cae_sne<- c(cae_sne,vals2)
    
    cae_beta_pod <- beta(CAE_selec,func = "sorensen" ,abund = F)
    vals0 <- as.vector(cae_beta_pod$Btotal); vals1 <-as.vector(cae_beta_pod$Brepl) ; vals2<-as.vector(cae_beta_pod$Brich)
    cae_total <- c(cae_total,vals0); cae_repl <- c(cae_repl,vals1); cae_rich <- c(cae_rich,vals2)
    
    message(i)
  }
  
  small_scale_sor <- data.frame(values = c(isl$sor, cae_sor), system = c(rep("isl",length(isl$sor)),rep("con",length(cae_sor))))
  kw_sor <- kruskal.test(small_scale_sor$values, small_scale_sor$system)
  
  small_scale_sim <- data.frame(values = c(isl$sim, cae_sim), system = c(rep("isl",length(isl$sim)),rep("con",length(cae_sim))))
  kw_sim <- kruskal.test(small_scale_sim$values, small_scale_sim$system)
  
  small_scale_sne <- data.frame(values = c(isl$sne, cae_sne), system = c(rep("isl",length(isl$sne)),rep("con",length(cae_sne))))
  kw_sne <- kruskal.test(small_scale_sne$values, small_scale_sne$system)
  
  
  # Podani
  small_scale_total <- data.frame(values = c(isl$total, cae_total), system = c(rep("isl",length(isl$total)),rep("con",length(cae_total))))
  kw_total <- kruskal.test(small_scale_total$values, small_scale_total$system)
  
  small_scale_repl <- data.frame(values = c(isl$repl, cae_repl), system = c(rep("isl",length(isl$repl)),rep("con",length(cae_repl))))
  kw_repl <- kruskal.test(small_scale_repl$values, small_scale_repl$system)
  
  small_scale_rich <- data.frame(values = c(isl$rich, cae_rich), system = c(rep("isl",length(isl$rich)),rep("con",length(cae_rich))))
  kw_rich <- kruskal.test(small_scale_rich$values, small_scale_rich$system)
  
  
  results <- data.frame(mean_caes = c(mean(cae_sor), mean(cae_sim), mean(cae_sne), mean(cae_total), mean(cae_repl), mean(cae_rich)),
                        sd_caes = c(sd(cae_sor), sd(cae_sim), sd(cae_sne), sd(cae_total), sd(cae_repl), sd(cae_rich)),
                        mean_island = c(mean(isl$sor), mean(isl$sim), mean(isl$sne), mean(isl$total), mean(isl$repl), mean(isl$rich)),
                        sd_island = c(sd(isl$sor), sd(isl$sim), sd(isl$sne), sd(isl$total), sd(isl$repl), sd(isl$rich)),
                        p_value = c(kw_sor$p.value, kw_sim$p.value, kw_sne$p.value, kw_total$p.value, kw_repl$p.value, kw_rich$p.value))
  
  rownames(results) <- c("sor (Baselga)","sim (Baselga)","sne (Baselga)", "total (Podani)", "repl (Podani)", "rich (Podani)")
  
  return(results)
}  




### Function for distance-decay fitting ###

  fit <- function(sim,dist){
    
    # Model fitting
    linear <- lm(sim~dist)                                                       # Linear model
    power_law <- nlsLM(sim~a*(dist^b),start =list(a=0,b=0),trace=F)              # Power_law
    exponential <- nlsLM(sim~a*exp(-b*dist),start =list(a=0,b=0),trace=F)        # Negative exponential
    gompertz <- nlsLM(sim~1-(exp(-a*exp(-b*dist))),start =list(a=0,b=0),trace=F) # Gompertz
    
    # AIC
    aic_lm <- AIC(linear); aic_pow <- AIC(power_law); aic_exp <- AIC(exponential); aic_gom <- AIC(gompertz)
    
    # Pseudo r2
    null_model <- lm(sim ~ 1)
    pseudo_lm <- 1-deviance(linear)/deviance(null_model) 
    pseudo_pow <- 1-deviance(power_law)/deviance(null_model)
    pseudo_exp <- 1-deviance(exponential)/deviance(null_model)
    pseudo_gom <- 1-deviance(gompertz)/deviance(null_model)
    
    results <- data.frame(aic = c(aic_lm, aic_pow, aic_exp, aic_gom), pseudo_r2 = c(pseudo_lm, pseudo_pow, pseudo_exp, pseudo_gom))
    rownames(results) <- c("lm","pow","exp","gom")
    
    return(results)
    
  }
  

  
### Function for the comparison of distance-decay models between CAEs and islands at small scale ###
  
  dd_small_scale <- function(isl,isl_id,cae_pa,cae_coor,cae_env,cae_id){ 
    
    sim <- 1-isl$sim
    
    pseudo_geodist_island <- fit(sim = sim, dist = isl$geodist)$pseudo_r2[3] 
    pseudo_envdist_island <- fit(sim = sim, dist = isl$envdist)$pseudo_r2[3]
    
    
    pseudo_geodist_cae <- matrix(nrow = 4, ncol = length(isl_id))
    pseudo_envdist_cae <- matrix(nrow = 4, ncol = length(isl_id))
    
    for (i in 1:length(isl_id)) {
      
      CAE_selec_pa <- cae_pa[which(cae_id$CAE_code %in% isl_id[i]),]
      CAE_selec_coor <- cae_coor[which(cae_id$CAE_code %in% isl_id[i]),]
      CAE_selec_env <- cae_env[which(cae_id$CAE_code %in% isl_id[i]),]
      
      cae_geodist <- as.dist(geodist(CAE_selec_coor, measure = "geodesic")/1000)
      
      # For the environmental distances a small scale is necessary to include this conditional as is some CAEs
      # there is no variance on precipitation seasonality among sites, so the PCA cannot be computed with taht variable
      if (sd(CAE_selec_env$prec_sea)==0){
        cae_envdist <- dist(principal(CAE_selec_env[,-6], nfactors = 2)$scores)
      } else {
        cae_envdist <- dist(principal(CAE_selec_env, nfactors = 2)$scores)
      }
      
      cae_similarity <- 1 - beta.pair(CAE_selec_pa)$beta.sim
      
      pseudo_geodist_cae[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_geodist))$pseudo_r2
      pseudo_envdist_cae[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_envdist))$pseudo_r2
      
      message(i)
    }
    
    propgeo <- mean(pseudo_geodist_island > pseudo_geodist_cae[3,])
    propenv <- mean(pseudo_envdist_island > pseudo_envdist_cae[3,])
    
    geodist_test <- prop.test(sum(pseudo_geodist_island > pseudo_geodist_cae[3,]), length(isl_id))
    envdist_test <- prop.test(sum(pseudo_envdist_island > pseudo_envdist_cae[3,]), length(isl_id))
    
    results <- data.frame(island_pseudo_r2 = c(pseudo_geodist_island, pseudo_envdist_island),
                          proportion = c(propgeo, propenv),
                          p_val = c(geodist_test$p.value, envdist_test$p.value))
    
    rownames(results) <- c("Spatial", "Environmental")
    
    return(results)
    
  }

  
#####

  


# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
# Beta diversity analyses
# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .

  
####
# Beta diversity analyses for the AMERICAN REGION
################
  
### Data  
## Islands 
  american_isla_coor <- american_islands[,4:5]                            # Coordinates of the sample sites on islands
  american_isla_env <- american_islands[,6:13]                            # Environmental variables of each sample site on islands
  american_isla_pa <- american_islands[,14:dim(american_islands)[2]]      # Presence/absence data of each site on islands

## CAEs
  american_cae_id <- american_caes[,1:3]                                # Id of different CAEs and their correspondent islands
  american_cae_coor <- american_caes[,4:5]                              # Coordinates of the sample sites contained on each CAE
  american_cae_env <- american_caes[,6:13]                              # Environmental variables of each sample site on CAEs
  american_cae_pa <- american_caes[,14:dim(american_caes)[2]]           # Presence/absence data of each site on CAEs
  
  
  
#### Beta diversity comparisons for the AMERICAN REGION at LARGE SCALE ###
  
 ## Beta diversity for islands at large scale
  
   # With Baselga's (bas) approach 
   american_isla_beta_bas <- beta.pair(american_isla_pa)
   mean(american_isla_beta_bas$beta.sor); sd(american_isla_beta_bas$beta.sor)
   mean(american_isla_beta_bas$beta.sim); sd(american_isla_beta_bas$beta.sim)
   mean(american_isla_beta_bas$beta.sne); sd(american_isla_beta_bas$beta.sne)
   
   # With Podani's (pod) approach
   american_isla_beta_pod <- beta(american_isla_pa, func = "sorensen", abund = F)
   mean(american_isla_beta_pod$Btotal); sd(american_isla_beta_pod$Btotal)
   mean(american_isla_beta_pod$Brepl); sd(american_isla_beta_pod$Brepl)
   mean(american_isla_beta_pod$Brich); sd(american_isla_beta_pod$Brich)
   
   
 ## Beta diversity for CAEs
   
  # Id of the CAEs correspondent to the islands in the region 
   Gu_id <- unique(subset(american_cae_id, Island_correspondence=="Gu")$CAE_code)
   Ha_id <- unique(subset(american_cae_id, Island_correspondence=="Ha")$CAE_code)
   Ka_id <- unique(subset(american_cae_id, Island_correspondence=="Ka")$CAE_code)
   Ma_id <- unique(subset(american_cae_id, Island_correspondence=="Ma")$CAE_code)
   Oa_id <- unique(subset(american_cae_id, Island_correspondence=="Oa")$CAE_code)
     
   
  # Vectors to store the results of the CAEs beta diversity measures at large scale   
   american_cae_sor <- c(); american_cae_sim <- c(); american_cae_sne <- c()     # Baselga
   american_cae_total <- c(); american_cae_repl <- c(); american_cae_rich <- c() # Podani
   
  # Loop for the 100 resaples of beta diversity between CAEs
   
   for (i in 1:100) {
     
     CAE_id <- c(sample(Gu_id, 1),sample(Ha_id,1),sample(Ka_id,1),sample(Ma_id,1),sample(Oa_id,1))
     CAE_selec <- (american_cae_pa[which(american_cae_id$CAE_code %in% CAE_id),])
     
     cae_beta_bas <- beta.pair(CAE_selec)
     vals0 <- as.vector(cae_beta_bas$beta.sor); vals1 <-as.vector(cae_beta_bas$beta.sim) ; vals2<-as.vector(cae_beta_bas$beta.sne)
     american_cae_sor <- c(american_cae_sor,vals0); american_cae_sim<- c(american_cae_sim,vals1); american_cae_sne<- c(american_cae_sne,vals2)
     
     cae_beta_pod <- beta(CAE_selec,func = "sorensen" ,abund = F)
     vals0 <- as.vector(cae_beta_pod$Btotal); vals1 <-as.vector(cae_beta_pod$Brepl) ; vals2<-as.vector(cae_beta_pod$Brich)
     american_cae_total <- c(american_cae_total,vals0); american_cae_repl <- c(american_cae_repl,vals1); american_cae_rich <- c(american_cae_rich,vals2)
     
     message(i)
   }
   

   # Results with Baselga's approach
   mean(american_cae_sor); sd(american_cae_sor)
   mean(american_cae_sim); sd(american_cae_sim)
   mean(american_cae_sne); sd(american_cae_sne)
   
   # Results with Podani's approach
   mean(american_cae_total); sd(american_cae_total)
   mean(american_cae_repl); sd(american_cae_repl)
   mean(american_cae_rich); sd(american_cae_rich)
   
   
  ## Testing for differences in beta diversity between CAEs and islands
   # Baselga
   america_sor <- data.frame(values = c(american_isla_beta_bas$beta.sor, american_cae_sor), system = c(rep("isl",length(american_isla_beta_bas$beta.sor)),rep("con",length(american_cae_sor))))
   kruskal.test(america_sor$values, america_sor$system)
   
   america_sim <- data.frame(values = c(american_isla_beta_bas$beta.sim, american_cae_sim), system = c(rep("isl",length(american_isla_beta_bas$beta.sim)),rep("con",length(american_cae_sim))))
   kruskal.test(america_sim$values, america_sim$system)
   
   america_sne <- data.frame(values = c(american_isla_beta_bas$beta.sne, american_cae_sne), system = c(rep("isl",length(american_isla_beta_bas$beta.sne)),rep("con",length(american_cae_sne))))
   kruskal.test(america_sne$values, america_sne$system)
   
   # Podani
   america_total <- data.frame(values = c(american_isla_beta_pod$Btotal, american_cae_total), system = c(rep("isl",length(american_isla_beta_pod$Btotal)),rep("con",length(american_cae_total))))
   kruskal.test(america_total$values, america_total$system)
   
   america_repl <- data.frame(values = c(american_isla_beta_pod$Brepl, american_cae_repl), system = c(rep("isl",length(american_isla_beta_pod$Brepl)),rep("con",length(american_cae_repl))))
   kruskal.test(america_repl$values, america_repl$system)
   
   america_rich <- data.frame(values = c(american_isla_beta_pod$Brich, american_cae_rich), system = c(rep("isl",length(american_isla_beta_pod$Brich)),rep("con",length(american_cae_rich))))
   kruskal.test(america_rich$values, america_rich$system)
   
   
### Partial Mantel test and Distance-decay models
   
 ## For islands
  
  # Compute spatial distance between sample sites
   american_isla_geodist <- as.dist(geodist(american_isla_coor, measure = "geodesic")/1000)
     
  # Compute spatial distance between sample sites 
   american_isla_envpca <- principal(american_isla_env, nfactors = 2)
   american_isla_envdist <- dist(american_isla_envpca$scores,method = "euclidean")
   

  # Distance-decay
   american_isla_similarity <- 1 - american_isla_beta_bas$beta.sim
   
   fit_american_isla_geodist <- fit(sim = as.vector(american_isla_similarity), dist = as.vector(american_isla_geodist))
   fit_american_isla_envdist <- fit(sim = as.vector(american_isla_similarity), dist = as.vector(american_isla_envdist))
    
    
 ## For CAEs
   
  # Vectors to store the results 
   
   dist_pseudo <- matrix(nrow = 4, ncol=200)
   env_pseudo <- matrix(nrow = 4, ncol=200)
    
   
  # Loop for the 200 distance-decay models and pseudo_r2 for CAEs at large scale  
   
    for (i in 1:200) {
      
      CAE_id <- c(sample(Gu_id, 1),sample(Ha_id,1),sample(Ka_id,1),sample(Ma_id,1),sample(Oa_id,1))
      
      CAE_selec <- american_cae_pa[which(american_cae_id$CAE_code %in%  CAE_id),] 
      CAE_selec_coor <- american_cae_coor[which(american_cae_id$CAE_code %in% CAE_id),] 
      CAE_selec_env <- american_cae_env[which(american_cae_id$CAE_code %in% CAE_id),] 
      
      cae_geodist <- as.dist(geodist(CAE_selec_coor, measure = "geodesic")/1000)
      cae_envdist <- dist(principal(CAE_selec_env, nfactors = 2)$scores)
      
      # Distance-decay
      cae_similarity <- 1 - beta.pair(CAE_selec)$beta.sim
      
      dist_pseudo[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_geodist))$pseudo_r2
      env_pseudo[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_envdist))$pseudo_r2
      
      message(i)
    }
    
   
 ## Comparison between Islands and CAEs
    
   gtest <- sum (fit_american_isla_geodist$pseudo_r2[3] > dist_pseudo[3,]); prop.test(gtest, 200)
   etest <- sum (fit_american_isla_geodist$pseudo_r2[3] > env_pseudo[3,]); prop.test(etest, 200)
    
   
#### Beta diversity comparisons for the AMERICAN REGION at SMALL SCALE ###
    
  ## Beta diversity within islands 
   
  # General data frame with the geographical, environmental distances and dissimilarity measures (Baselga and Podani) for all pairwise comparisons between American islands sites   
   pairwise_american_islands <- as.data.frame(t(combn(american_islands$Island_id,2))); pairwise_IAM <- paste(pairwise_american_islands$V1,pairwise_american_islands$V2,sep="-")
   distances_american_islands <- data.frame(ID=pairwise_IAM, geodist = as.vector(american_isla_geodist), envdist = as.vector(american_isla_envdist),
                                            sor = as.vector(american_isla_beta_bas$beta.sor), sim = as.vector(american_isla_beta_bas$beta.sim), sne = as.vector(american_isla_beta_bas$beta.sne),
                                            total = as.vector(american_isla_beta_pod$Btotal), repl = as.vector(american_isla_beta_pod$Brepl), rich = as.vector(american_isla_beta_pod$Brich))
   
 ## Data frame with the within islands beta diversity values for each island
   # Guadeloupe (Gu)
   Gu <- subset(distances_american_islands, distances_american_islands$ID=="Gu-Gu")
   # Hawaii (Ha)
   Ha <- subset(distances_american_islands, distances_american_islands$ID=="Ha-Ha")
   # Kauai (Ka)
   Ka <- subset(distances_american_islands, distances_american_islands$ID=="Ka-Ka")
   # Martinique (Ma)
   Ma <- subset(distances_american_islands, distances_american_islands$ID=="Ma-Ma")
   # Oahu (Oa)
   Oa <- subset(distances_american_islands, distances_american_islands$ID=="Oa-Oa")

   
   
 ## Beta diversity comparison between CAEs and islands at small scale (within islands / within CAEs)
   # and beta ratio for CAEs
  
  # Guadeloupe
  beta_small_scale(isl = Gu, isl_id = Gu_id, cae_pa = american_cae_pa, cae_id = american_cae_id)
  # Hawaii
  beta_small_scale(isl = Ha, isl_id = Ha_id, cae_pa = american_cae_pa, cae_id = american_cae_id)
  # Kauai
  beta_small_scale(isl = Ka, isl_id = Ka_id, cae_pa = american_cae_pa, cae_id = american_cae_id)
  # Martinique
  beta_small_scale(isl = Ma, isl_id = Ma_id, cae_pa = american_cae_pa, cae_id = american_cae_id)
  # Oahu
  beta_small_scale(isl = Oa, isl_id = Oa_id, cae_pa = american_cae_pa, cae_id = american_cae_id)
  
  
  
 ## Distance-decay models for islands and their comparison with CAEs at small scale (within islands / within CAEs)
  
  # Guadeloupe
  dd_small_scale(isl = Gu, isl_id = Gu_id, cae_pa = american_cae_pa, cae_coor = american_cae_coor, cae_env = american_cae_env, cae_id = american_cae_id) 
  # Hawaii
  dd_small_scale(isl = Ha, isl_id = Ha_id, cae_pa = american_cae_pa, cae_coor = american_cae_coor, cae_env = american_cae_env, cae_id = american_cae_id)
  # Kauai
  dd_small_scale(isl = Ka, isl_id = Ka_id, cae_pa = american_cae_pa, cae_coor = american_cae_coor, cae_env = american_cae_env, cae_id = american_cae_id)
  # Martinique
  dd_small_scale(isl = Ma, isl_id = Ma_id, cae_pa = american_cae_pa, cae_coor = american_cae_coor, cae_env = american_cae_env, cae_id = american_cae_id)
  # Oahu
  dd_small_scale(isl = Oa, isl_id = Oa_id, cae_pa = american_cae_pa, cae_coor = american_cae_coor, cae_env = american_cae_env, cae_id = american_cae_id)
  
######  
  
  

####
# Beta diversity analyses for the EUROPEAN REGION
################
  
### Data  
  ## Islands 
  european_isla_coor <- european_islands[,4:5]                            # Coordinates of the sample sites on islands
  european_isla_env <- european_islands[,6:13]                            # Environmental variables of each sample site on islands
  european_isla_pa <- european_islands[,14:dim(european_islands)[2]]      # Presence/absence data of each site on islands
  
  ## CAEs  
  european_cae_id <- european_caes[,1:3]                              # Id of different CAEs and their correspondent islands
  european_cae_coor <- european_caes[,4:5]                            # Coordinates of the sample sites contained on each CAE
  european_cae_env <- european_caes[,6:13]                            # Environmental variables of each sample site on CAEs
  european_cae_pa <- european_caes[,14:dim(european_caes)[2]]         # Presence/absence data of each site on CAEs
  
  
  
  #### Beta diversity comparisons for the EUROPEAN REGION at LARGE SCALE ###
  
  ## Beta diversity for islands
  
  # With Baselga's (bas) approach 
  european_isla_beta_bas <- beta.pair(european_isla_pa)
  mean(european_isla_beta_bas$beta.sor); sd(european_isla_beta_bas$beta.sor)
  mean(european_isla_beta_bas$beta.sim); sd(european_isla_beta_bas$beta.sim)
  mean(european_isla_beta_bas$beta.sne); sd(european_isla_beta_bas$beta.sne)
  
  # With Podani's (pod) approach
  european_isla_beta_pod <- beta(european_isla_pa, func = "sorensen", abund = F)
  mean(european_isla_beta_pod$Btotal); sd(european_isla_beta_pod$Btotal)
  mean(european_isla_beta_pod$Brepl); sd(european_isla_beta_pod$Brepl)
  mean(european_isla_beta_pod$Brich); sd(european_isla_beta_pod$Brich)
  
  
  ## Beta diversity for CAEs
  
  # Id of the CAEs correspondent to the islands in the region 
  Co_id <- unique(subset(european_cae_id, Island_correspondence=="Co")$CAE_code) # Corsica
  Cy_id <- unique(subset(european_cae_id, Island_correspondence=="Cy")$CAE_code) # Cyprus
  Ic_id <- unique(subset(european_cae_id, Island_correspondence=="Ic")$CAE_code) # Iceland
  Ir_id <- unique(subset(european_cae_id, Island_correspondence=="Ir")$CAE_code) # Ireland
  Md_id <- unique(subset(european_cae_id, Island_correspondence=="Md")$CAE_code) # Madeira
  Mj_id <- unique(subset(european_cae_id, Island_correspondence=="Mj")$CAE_code) # Majorca
  Sa_id <- unique(subset(european_cae_id, Island_correspondence=="Sa")$CAE_code) # Sardinia
  Sm_id <- unique(subset(european_cae_id, Island_correspondence=="Sm")$CAE_code) # Sao Miguel
  
  
  # Vectors to store the results of the CAEs beta diversity measures at large scale   
  european_cae_sor <- c(); european_cae_sim <- c(); european_cae_sne <- c()     # Baselga
  european_cae_total <- c(); european_cae_repl <- c(); european_cae_rich <- c() # Podani
  
  
  # Loop for the 100 resaples of beta diversity between CAEs
  
  for (i in 1:100) {
    
    CAE_id <- c(sample(Co_id, 1),sample(Cy_id,1),sample(Ic_id,1),sample(Ir_id,1),sample(Md_id,1),sample(Mj_id, 1),sample(Sa_id,1),sample(Sm_id,1))
    CAE_selec <- (european_cae_pa[which(european_cae_id$CAE_code %in% CAE_id),])
    
    cae_beta_bas <- beta.pair(CAE_selec)
    vals0 <- as.vector(cae_beta_bas$beta.sor); vals1 <-as.vector(cae_beta_bas$beta.sim) ; vals2<-as.vector(cae_beta_bas$beta.sne)
    european_cae_sor <- c(european_cae_sor,vals0); european_cae_sim<- c(european_cae_sim,vals1); european_cae_sne<- c(european_cae_sne,vals2)
    
    cae_beta_pod <- beta(CAE_selec,func = "sorensen" ,abund = F)
    vals0 <- as.vector(cae_beta_pod$Btotal); vals1 <-as.vector(cae_beta_pod$Brepl) ; vals2<-as.vector(cae_beta_pod$Brich)
    european_cae_total <- c(european_cae_total,vals0); european_cae_repl <- c(european_cae_repl,vals1); european_cae_rich <- c(european_cae_rich,vals2)
    
    message(i)
  }
  
  
  # Results with Baselga's approach
  mean(european_cae_sor); sd(european_cae_sor)
  mean(european_cae_sim); sd(european_cae_sim)
  mean(european_cae_sne); sd(european_cae_sne)
  
  # Results with Podani's approach
  mean(european_cae_total); sd(european_cae_total)
  mean(european_cae_repl); sd(european_cae_repl)
  mean(european_cae_rich); sd(european_cae_rich)
  

  ## Test for differences 
  # Baselga
  europe_sor <- data.frame(values = c(european_isla_beta_bas$beta.sor, european_cae_sor), system = c(rep("isl",length(european_isla_beta_bas$beta.sor)),rep("con",length(european_cae_sor))))
  kruskal.test(europe_sor$values, europe_sor$system)
  
  europe_sim <- data.frame(values = c(european_isla_beta_bas$beta.sim, european_cae_sim), system = c(rep("isl",length(european_isla_beta_bas$beta.sim)),rep("con",length(european_cae_sim))))
  kruskal.test(europe_sor$values, europe_sor$system)
  
  europe_sne <- data.frame(values = c(european_isla_beta_bas$beta.sne, european_cae_sne), system = c(rep("isl",length(european_isla_beta_bas$beta.sne)),rep("con",length(european_cae_sne))))
  kruskal.test(europe_sor$values, europe_sor$system)
  
  
  # Podani
  europe_total <- data.frame(values = c(european_isla_beta_pod$Btotal, european_cae_total), system = c(rep("isl",length(european_isla_beta_pod$Btotal)),rep("con",length(european_cae_total))))
  kruskal.test(europe_total$values, europe_total$system)
  
  europe_repl <- data.frame(values = c(european_isla_beta_pod$Brepl, european_cae_repl), system = c(rep("isl",length(european_isla_beta_pod$Brepl)),rep("con",length(european_cae_repl))))
  kruskal.test(europe_repl$values, europe_repl$system)
  
  europe_rich <- data.frame(values = c(european_isla_beta_pod$Brich, european_cae_rich), system = c(rep("isl",length(european_isla_beta_pod$Brich)),rep("con",length(european_cae_rich))))
  kruskal.test(europe_rich$values, europe_rich$system)
  
  
  ### Partial Mantel test and Distance-decay models
  
  ## For islands
  
  # Compute spatial distance between sample sites
  european_isla_geodist <- as.dist(geodist(european_isla_coor, measure = "geodesic")/1000)
  
  # Compute spatial distance between sample sites 
  european_isla_envpca <- principal(european_isla_env, nfactors = 2)
  european_isla_envdist <- dist(european_isla_envpca$scores,method = "euclidean")
  
  
  # Distance-decay
  european_isla_similarity <- 1 - european_isla_beta_bas$beta.sim
  
  fit_european_isla_geodist <- fit(sim = as.vector(european_isla_similarity), dist = as.vector(european_isla_geodist))
  fit_european_isla_envdist <- fit(sim = as.vector(european_isla_similarity), dist = as.vector(european_isla_envdist))
  
  
  ## For CAEs
  
  # Vectors to store the results 
  
  dist_pseudo <- matrix(nrow = 4, ncol=200)
  env_pseudo <- matrix(nrow = 4, ncol=200)
  
  # Loop for the 200 distance-decay models and pseudo_r2 for CAEs at large scale  
  
  for (i in 1:200) {
    
    CAE_id <- c(sample(Gu_id, 1),sample(Ha_id,1),sample(Ka_id,1),sample(Ma_id,1),sample(Oa_id,1))
    
    CAE_selec <- european_cae_pa[which(european_cae_id$CAE_code %in%  CAE_id),] 
    CAE_selec_coor <- european_cae_coor[which(european_cae_id$CAE_code %in% CAE_id),] 
    CAE_selec_env <- european_cae_env[which(european_cae_id$CAE_code %in% CAE_id),] 
    
    cae_geodist <- as.dist(geodist(CAE_selec_coor, measure = "geodesic")/1000)
    cae_envdist <- dist(principal(CAE_selec_env, nfactors = 2)$scores)
    
    # Distance-decay
    cae_similarity <- 1 - beta.pair(CAE_selec)$beta.sim
    
    dist_pseudo[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_geodist))$pseudo_r2
    env_pseudo[,i] <- fit(sim = as.vector(cae_similarity), dist = as.vector(cae_envdist))$pseudo_r2
    
    message(i)
  }
  
  
  ## Comparison between Islands and CAEs
  
  gtest <- sum (fit_european_isla_geodist$pseudo_r2[3] > dist_pseudo[3,]); prop.test(gtest, 200)
  etest <- sum (fit_european_isla_geodist$pseudo_r2[3] > env_pseudo[3,]); prop.test(etest, 200)
  
  
  
#### Beta diversity comparisons for the EUROPEAN REGION at SMALL SCALE ###
  
  ## Beta diversity within islands 
  
  # General data frame with the geographical, environmental distances and dissimilarity measures (Baselga and Podani) for all pairwise comparisons on European islands   
  pairwise_european_islands <- as.data.frame(t(combn(european_islands$Island_id,2))); pairwise_IAM <- paste(pairwise_european_islands$V1,pairwise_european_islands$V2,sep="-")
  distances_european_islands <- data.frame(ID=pairwise_IAM, geodist = as.vector(european_isla_geodist), envdist = as.vector(european_isla_envdist),
                                           sor = as.vector(european_isla_beta_bas$beta.sor), sim = as.vector(european_isla_beta_bas$beta.sim), sne = as.vector(european_isla_beta_bas$beta.sne),
                                           total = as.vector(european_isla_beta_pod$Btotal), repl = as.vector(european_isla_beta_pod$Brepl), rich = as.vector(european_isla_beta_pod$Brich))
  
  
  ## Data frame with the within islands beta diversity values for each island
  # Corsica (Co)
  Co <- subset(distances_european_islands, distances_european_islands$ID=="Co-Co")
  # Cyprus (Cy)
  Cy <- subset(distances_european_islands, distances_european_islands$ID=="Cy-Cy")
  # Iceland (Ic)
  Ic <- subset(distances_european_islands, distances_european_islands$ID=="Ic-Ic")
  # Ireland (Ma)
  Ir <- subset(distances_european_islands, distances_european_islands$ID=="Ir-Ir")
  # Madeira (Md)
  Md <- subset(distances_european_islands, distances_european_islands$ID=="Md-Md")
  # Majorca (Mj)
  Mj <- subset(distances_european_islands, distances_european_islands$ID=="Mj-Mj")
  # Sardinia (Sa)
  Sa <- subset(distances_european_islands, distances_european_islands$ID=="Sa-Sa")
  # Sao Miguel (Sm)
  Sm <- subset(distances_european_islands, distances_european_islands$ID=="Sm-Sm")
  
  
  ## Beta diversity comparison between CAEs and islands at small scale (within islands / within CAEs)
  # and beta ratio for CAEs
  
  # Corsica
  beta_small_scale(isl = Co, isl_id = Co_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Cyprus
  beta_small_scale(isl = Cy, isl_id = Cy_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Iceland
  beta_small_scale(isl = Ic, isl_id = Ic_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Ireland
  beta_small_scale(isl = Ir, isl_id = Ir_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Madeira
  beta_small_scale(isl = Md, isl_id = Md_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Majorca
  beta_small_scale(isl = Mj, isl_id = Mj_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Sardinia
  beta_small_scale(isl = Sa, isl_id = Sa_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  # Sao Miguel
  beta_small_scale(isl = Sm, isl_id = Sm_id, cae_pa = european_cae_pa, cae_id = european_cae_id)
  
  
  ## Distance-decay models for islands and their comparison with CAEs at small scale (within islands / within CAEs)
  
  # Corsica
  dd_small_scale(isl = Co, isl_id = Co_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id) 
  # Cyprus
  dd_small_scale(isl = Cy, isl_id = Cy_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Iceland
  dd_small_scale(isl = Ic, isl_id = Ic_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Ireland
  dd_small_scale(isl = Ir, isl_id = Ir_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Madeira
  dd_small_scale(isl = Md, isl_id = Md_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Majorca
  dd_small_scale(isl = Mj, isl_id = Mj_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Sardinia
  dd_small_scale(isl = Sa, isl_id = Sa_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  # Sao Muguel
  dd_small_scale(isl = Sm, isl_id = Sm_id, cae_pa = european_cae_pa, cae_coor = european_cae_coor, cae_env = european_cae_env, cae_id = european_cae_id)
  
#####
  
  
  
####
# Beta diversity analyses for the AFRICAN REGION
################
  
## Data  
  ## Islands 
  african_isla_coor <- african_islands[,4:5]                           # Coordinates of the sample sites on islands
  african_isla_env <- african_islands[,6:13]                           # Environmental variables of each sample site on islands
  african_isla_pa <- african_islands[,14:dim(african_islands)[2]]      # Presence/absence data of each site on islands
  
  ## CAEs
  african_cae_id <- african_caes[,1:3]                            # Id of different CAEs and their correspondent islands
  african_cae_coor <- african_caes[,4:5]                          # Coordinates of the sample sites contained on each CAE
  african_cae_env <- african_caes[,6:13]                          # Environmental variables of each sample site on CAEs
  african_cae_pa <- african_caes[,14:dim(african_caes)[2]]        # Presence/absence data of each site on CAEs
  
  
  
  ## Beta diversity for islands at LARGE SCALE
  
  # With Baselga's (bas) approach 
  african_isla_beta_bas <- beta.pair(african_isla_pa)
  mean(african_isla_beta_bas$beta.sor); sd(african_isla_beta_bas$beta.sor)
  mean(african_isla_beta_bas$beta.sim); sd(african_isla_beta_bas$beta.sim)
  mean(african_isla_beta_bas$beta.sne); sd(african_isla_beta_bas$beta.sne)
  
  # With Podani's (pod) approach
  african_isla_beta_pod <- beta(african_isla_pa, func = "sorensen", abund = F)
  mean(african_isla_beta_pod$Btotal); sd(african_isla_beta_pod$Btotal)
  mean(african_isla_beta_pod$Brepl); sd(african_isla_beta_pod$Brepl)
  mean(african_isla_beta_pod$Brich); sd(african_isla_beta_pod$Brich)
  
  
  ### Partial Mantel test
  
  ## For islands
  # Compute spatial distance between sample sites
  african_isla_geodist <- as.dist(geodist(african_isla_coor, measure = "geodesic")/1000)
  
  # Compute spatial distance between sample sites 
  african_isla_envpca <- principal(african_isla_env, nfactors = 2)
  african_isla_envdist <- dist(african_isla_envpca$scores,method = "euclidean")
  
  
  # Id of the CAEs correspondent to the islands on the region 
  My_id <- unique(subset(african_cae_id, Island_correspondence=="My")$CAE_code) # Mayotte
  Po_id <- unique(subset(african_cae_id, Island_correspondence=="Po")$CAE_code) # Posession
  Re_id <- unique(subset(african_cae_id, Island_correspondence=="Re")$CAE_code) # Reunion
  
  
  
#### Beta diversity comparisons for the AFRICAN REGION at SMALL SCALE ###
  
  ## Beta diversity within islands 
  
  # General data frame with the geographical, environmental distances and dissimilarity measures (Baselga and Podani) for all pairwise comparisons on AFRICAN islands   
  pairwise_african_islands <- as.data.frame(t(combn(african_islands$Island_id,2))); pairwise_IAM <- paste(pairwise_african_islands$V1,pairwise_african_islands$V2,sep="-")
  distances_african_islands <- data.frame(ID=pairwise_IAM, geodist = as.vector(african_isla_geodist), envdist = as.vector(african_isla_envdist),
                                           sor = as.vector(african_isla_beta_bas$beta.sor), sim = as.vector(african_isla_beta_bas$beta.sim), sne = as.vector(african_isla_beta_bas$beta.sne),
                                           total = as.vector(african_isla_beta_pod$Btotal), repl = as.vector(african_isla_beta_pod$Brepl), rich = as.vector(african_isla_beta_pod$Brich))
  
  
  ## Mean beta diversity within islands
  # Mayotte (My)
  My <- subset(distances_african_islands, distances_african_islands$ID=="My-My")
  # Posession (Po)
  Po <- subset(distances_african_islands, distances_african_islands$ID=="Po-Po")
  # Reunion (Re)
  Re <- subset(distances_african_islands, distances_african_islands$ID=="Re-Re")
  
  
  
  ## Beta diversity comparison between CAEs and islands at small scale (within islands / within CAEs)
  # and beta ratio for CAEs
  
  # Mayotte
  beta_small_scale(isl = My, isl_id = My_id, cae_pa = african_cae_pa, cae_id = african_cae_id)
  # Posession
  beta_small_scale(isl = Po, isl_id = Po_id, cae_pa = african_cae_pa, cae_id = african_cae_id)
  # Reunion
  beta_small_scale(isl = Re, isl_id = Re_id, cae_pa = african_cae_pa, cae_id = african_cae_id)

  
  
  ## Distance-decay models for islands and their comparison with CAEs at small scale (within islands / within CAEs)
  
  # Mayotte
  dd_small_scale(isl = My, isl_id = My_id, cae_pa = african_cae_pa, cae_coor = african_cae_coor, cae_env = african_cae_env, cae_id = african_cae_id) 
  # Posession
  dd_small_scale(isl = Po, isl_id = Po_id, cae_pa = african_cae_pa, cae_coor = african_cae_coor, cae_env = african_cae_env, cae_id = african_cae_id)
  # Reunion
  dd_small_scale(isl = Re, isl_id = Re_id, cae_pa = african_cae_pa, cae_coor = african_cae_coor, cae_env = african_cae_env, cae_id = african_cae_id)
  
  
#####
  
  
  
####
# Beta diversity analyses for the PACIFIC REGION
################  
  
  
  ## Islands 
  pacific_isla_coor <- pacific_islands[,4:5]                           # Coordinates of the sample sites on islands
  pacific_isla_env <- pacific_islands[,6:13]                           # Environmental variables of each sample site on islands
  pacific_isla_pa <- pacific_islands[,14:dim(pacific_islands)[2]]      # Presence/absence data of each site on islands
  
  
  #### Beta diversity comparisons for the PACIFIC REGION at LARGE SCALE ###
  
  ## Beta diversity for islands
  
  # With Baselga's (bas) approach 
  pacific_isla_beta_bas <- beta.pair(pacific_isla_pa)
  mean(pacific_isla_beta_bas$beta.sor); sd(pacific_isla_beta_bas$beta.sor)
  mean(pacific_isla_beta_bas$beta.sim); sd(pacific_isla_beta_bas$beta.sim)
  mean(pacific_isla_beta_bas$beta.sne); sd(pacific_isla_beta_bas$beta.sne)
  
  # With Podani's (pod) approach
  pacific_isla_beta_pod <- beta(pacific_isla_pa, func = "sorensen", abund = F)
  mean(pacific_isla_beta_pod$Btotal); sd(pacific_isla_beta_pod$Btotal)
  mean(pacific_isla_beta_pod$Brepl); sd(pacific_isla_beta_pod$Brepl)
  mean(pacific_isla_beta_pod$Brich); sd(pacific_isla_beta_pod$Brich)
  
  
  ### Partial Mantel test and Distance-decay models
  
  ## For islands
  
  # Compute spatial distance between sample sites
  pacific_isla_geodist <- as.dist(geodist(pacific_isla_coor, measure = "geodesic")/1000)
  
  # Compute spatial distance between sample sites 
  pacific_isla_envpca <- principal(pacific_isla_env, nfactors = 2)
  pacific_isla_envdist <- dist(pacific_isla_envpca$scores,method = "euclidean")
  
  # Partial Mantel test
  mantel_europe_isl_spatial <- mantel.partial(pacific_isla_beta_bas$beta.sim, pacific_isla_geodist, pacific_isla_envdist, permutations = 1000)
  mantel_europe_isl_spatial$statistic; mantel_europe_isl_spatial$signif
  mantel_europe_isl_environmental <- mantel.partial(pacific_isla_beta_bas$beta.sim, pacific_isla_envdist, pacific_isla_geodist, permutations = 1000)
  mantel_europe_isl_environmental$statistic; mantel_europe_isl_environmental$signif
  
  # Distance-decay
  pacific_isla_similarity <- 1 - pacific_isla_beta_bas$beta.sim
  
  fit_pacific_isla_geodist <- fit(sim = as.vector(pacific_isla_similarity), dist = as.vector(pacific_isla_geodist))
  fit_pacific_isla_geodist
  fit_pacific_isla_envdist <- fit(sim = as.vector(pacific_isla_similarity), dist = as.vector(pacific_isla_envdist))
  fit_pacific_isla_envdist
  
  
  
  #### Beta diversity comparisons for the pacific REGION at SMALL SCALE ###
  
  ## Beta diversity within islands 
  
  # General data frame with the geographical, environmental distances and dissimilarity measures (Baselga and Podani) for all pairwise comparisons on European islands   
  pairwise_pacific_islands <- as.data.frame(t(combn(pacific_islands$Island_id,2))); pairwise_IAM <- paste(pairwise_pacific_islands$V1,pairwise_pacific_islands$V2,sep="-")
  distances_pacific_islands <- data.frame(ID=pairwise_IAM, geodist = as.vector(pacific_isla_geodist), envdist = as.vector(pacific_isla_envdist),
                                           sor = as.vector(pacific_isla_beta_bas$beta.sor), sim = as.vector(pacific_isla_beta_bas$beta.sim), sne = as.vector(pacific_isla_beta_bas$beta.sne),
                                           total = as.vector(pacific_isla_beta_pod$Btotal), repl = as.vector(pacific_isla_beta_pod$Brepl), rich = as.vector(pacific_isla_beta_pod$Brich))
  
  
  ## Mean beta diversity within islands
  # New Caledonia (Nc)
  Nc <- subset(distances_pacific_islands, distances_pacific_islands$ID=="Nc-Nc")
  # North New Zealand (Nz)
  Nz <- subset(distances_pacific_islands, distances_pacific_islands$ID=="Nz-Nz")
  # South New Zealand (Sz)
  Sz <- subset(distances_pacific_islands, distances_pacific_islands$ID=="Sz-Sz")
  
  
  # Beta diversity for the Pacific islands at small scale (within islands)
  
  # New Caledonia
  Nc_beta <- data.frame(mean_island = c(mean(Nc$sor), mean(Nc$sim), mean(Nc$sne), mean(Nc$total), mean(Nc$repl), mean(Nc$rich)),
             sd_island = c(sd(Nc$sor), sd(Nc$sim), sd(Nc$sne), sd(Nc$total), sd(Nc$repl), sd(Nc$rich))); rownames(Nc_beta) <- c("sor (Baselga)","sim (Baselga)","sne (Baselga)", "total (Podani)", "repl (Podani)", "rich (Podani)")
  Nc_beta
  # North New Zealand
  Nz_beta <- data.frame(mean_island = c(mean(Nz$sor), mean(Nz$sim), mean(Nz$sne), mean(Nz$total), mean(Nz$repl), mean(Nz$rich)),
                        sd_island = c(sd(Nz$sor), sd(Nz$sim), sd(Nz$sne), sd(Nz$total), sd(Nz$repl), sd(Nz$rich))); rownames(Nz_beta) <- c("sor (Baselga)","sim (Baselga)","sne (Baselga)", "total (Podani)", "repl (Podani)", "rich (Podani)")
  Nz_beta
  # South New Zealand
  Sz_beta <- data.frame(mean_island = c(mean(Sz$sor), mean(Sz$sim), mean(Sz$sne), mean(Sz$total), mean(Sz$repl), mean(Sz$rich)),
                        sd_island = c(sd(Sz$sor), sd(Sz$sim), sd(Sz$sne), sd(Sz$total), sd(Sz$repl), sd(Sz$rich))); rownames(Sz_beta) <- c("sor (Baselga)","sim (Baselga)","sne (Baselga)", "total (Podani)", "repl (Podani)", "rich (Podani)")
  Sz_beta
  
  
  ### Distance-decay models for islands at small scale (within islands)
  
  # New Caledonia
  fit(sim = 1-Nc$sor, dist = Nc$geodist) # Spatial 
  fit(sim = 1-Nc$sor, dist = Nc$envdist) # Environmental
  # North New Zealand
  fit(sim = 1-Nz$sor, dist = Nz$geodist) # Spatial 
  fit(sim = 1-Nz$sor, dist = Nz$envdist) # Environmental
  # South New Zealand
  fit(sim = 1-Sz$sor, dist = Sz$geodist) # Spatial 
  fit(sim = 1-Sz$sor, dist = Sz$envdist) # Environmental
  
  
#####
  
  
  
  

# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .
# RELATIONSHIP BETWEEN ISLANDS BETA DIVERSITY AND THEIR ATTRIBUTES
# ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ . ~ .  
  
  
####
# Multiple beta diversity and islands characteristics analyses
################ 
  
 detach("package:BAT", unload = TRUE) # Unload the package BAT because it overlaps with the beta.multi function of
                                      # the package betapart, used in these analyses 
  
## Computing multiple beta diversity for each island 
 # American region
 Gu_multi <- beta.multi(american_isla_pa[which(american_islands$Island_id=="Gu"),])
 Ha_multi <- beta.multi(american_isla_pa[which(american_islands$Island_id=="Ha"),])
 Ka_multi <- beta.multi(american_isla_pa[which(american_islands$Island_id=="Ka"),])
 Ma_multi <- beta.multi(american_isla_pa[which(american_islands$Island_id=="Ma"),])
 Oa_multi <- beta.multi(american_isla_pa[which(american_islands$Island_id=="Oa"),])
 # European region
 Co_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Co"),])
 Cy_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Cy"),])
 Ic_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Ic"),])
 Ir_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Ir"),])
 Md_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Md"),])
 Mj_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Mj"),])
 Sa_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Sa"),])
 Sm_multi <- beta.multi(european_isla_pa[which(european_islands$Island_id=="Sm"),])
 # African region
 My_multi <- beta.multi(african_isla_pa[which(african_islands$Island_id=="My"),])
 Po_multi <- beta.multi(african_isla_pa[which(african_islands$Island_id=="Po"),])
 Re_multi <- beta.multi(african_isla_pa[which(african_islands$Island_id=="Re"),])
 # Pacific region
 Nc_multi <- beta.multi(pacific_isla_pa[which(pacific_islands$Island_id=="Nc"),])
 Nz_multi <- beta.multi(pacific_isla_pa[which(pacific_islands$Island_id=="Nz"),])
 Sz_multi <- beta.multi(pacific_isla_pa[which(pacific_islands$Island_id=="Sz"),])
 
 
 
 # Building the data frame with the multiple beta diversity results and islands characteristics
 
 beta_multiple_values <- data.frame(Island_id = c("Gu","Ha","Ka","Ma","Oa","Co","Cy","Ic","Ir","Md","Mj", "Sa","Sm","My","Po","Re", "Nc","Nz","Sz"),
   
                                    SOR = c(Gu_multi$beta.SOR, Ha_multi$beta.SOR, Ka_multi$beta.SOR, Ma_multi$beta.SOR, Oa_multi$beta.SOR, Co_multi$beta.SOR, Cy_multi$beta.SOR, Ic_multi$beta.SOR,
                                            Ir_multi$beta.SOR, Md_multi$beta.SOR, Mj_multi$beta.SOR, Sa_multi$beta.SOR, Sm_multi$beta.SOR, My_multi$beta.SOR, Po_multi$beta.SOR, Re_multi$beta.SOR,
                                            Nc_multi$beta.SOR, Nz_multi$beta.SOR, Sz_multi$beta.SOR),
                                    
                                    SIM = c(Gu_multi$beta.SIM, Ha_multi$beta.SIM, Ka_multi$beta.SIM, Ma_multi$beta.SIM, Oa_multi$beta.SIM, Co_multi$beta.SIM, Cy_multi$beta.SIM, Ic_multi$beta.SIM,
                                            Ir_multi$beta.SIM, Md_multi$beta.SIM, Mj_multi$beta.SIM, Sa_multi$beta.SIM, Sm_multi$beta.SIM, My_multi$beta.SIM, Po_multi$beta.SIM, Re_multi$beta.SIM,
                                            Nc_multi$beta.SIM, Nz_multi$beta.SIM, Sz_multi$beta.SIM),
                                    
                                    SNE = c(Gu_multi$beta.SNE, Ha_multi$beta.SNE, Ka_multi$beta.SNE, Ma_multi$beta.SNE, Oa_multi$beta.SNE, Co_multi$beta.SNE, Cy_multi$beta.SNE, Ic_multi$beta.SNE,
                                            Ir_multi$beta.SNE, Md_multi$beta.SNE, Mj_multi$beta.SNE, Sa_multi$beta.SNE, Sm_multi$beta.SNE, My_multi$beta.SNE, Po_multi$beta.SNE, Re_multi$beta.SNE,
                                            Nc_multi$beta.SNE, Nz_multi$beta.SNE, Sz_multi$beta.SNE))
                                    
 islands_df <- merge(island_char, beta_multiple_values, by = "Island_id", sort = F)                                    
 
 
 ##  Computing the linear models
 
 # Models for variable selection
 multi_sor_select <- lm(SOR ~ mean_lat + size_sample + age + isolation, data = islands_df, na.action = "na.fail")
 dredge(multi_sor_select)
 
 multi_sim_select <- lm(SIM ~ mean_lat + size_sample + age + isolation, data = islands_df, na.action = "na.fail")
 dredge(multi_sim_select)
 
 multi_sne_select <- lm(SNE ~ mean_lat + size_sample + age + isolation, data = islands_df, na.action = "na.fail")
 dredge(multi_sne_select)
 
 
 # Mix linear models for testing the effect of the regions on islands beta diversity taking into account 
 # the variables selected previously
 mix_sor_model <- lmer(SOR ~ mean_lat + age + isolation + (1 | Region), data = islands_df)
 exactRLRT(mix_sor_model)
 
 mix_sim_model <- lmer(SIM ~ mean_lat + age + isolation + (1 | Region), data = islands_df)
 exactRLRT(mix_sim_model)
 
 mix_sne_model <- lmer(SIM ~ mean_lat + age + isolation + (1 | Region), data = islands_df)
 exactRLRT(mix_sne_model)
 
 
 # Final models
 final_sor_model <- lm(SOR ~ mean_lat + age + isolation, data = islands_df)
 summary(final_sor_model)
 
 final_sim_model <- lm(SIM ~ mean_lat + age + isolation, data = islands_df)
 summary(final_sim_model)
 
 final_sne_model <- lm(SNE ~ mean_lat + age + isolation, data = islands_df)
 summary(final_sne_model)
 
 