##################################################################################################################################################
### Disentangling fine- and large-scale colonization processes in metapopulation dynamics: a case study on a threatenedepiphytic bryophyte" 
### Script R - Description of metapopulation dynamic at stand scale 
### G. PERCEL & F. LAROCHE (June 2023)
##################################################################################################################################################


#*****************************************
# Installing & Loading R-packages
#*****************************************

#install.packages("")
library(maptools)
library(ggplot2)


#************************************
# Loading stand data (local-scale)
#*************************************

root = "C:/Users/HP050/Documents/Metapop_colo_Octobre2017_4/09_ArticleBryo_SoumisFL.170223/Revision1/"
data_DV17 <- read.table(paste(root,"Data_Dviride.2017_AVCHESS_d14.txt",sep=""),dec=",",header=T)
View(data_DV17)
data_DV17$PARCELLE_ID <- as.factor(data_DV17$PARCELLE_ID)
data_DV17$SUBSTRAT_ID <- as.factor(data_DV17$SUBSTRAT_ID)
data_DV17$DIAMETRE <- as.numeric(data_DV17$DIAMETRE)
data_DV17$DETECT14 <- as.numeric(data_DV17$DETECT14)
data_DV17$DETECT17 <- as.numeric(data_DV17$DETECT17)
str(data_DV17)

# Dataset simplification and cleaning
data_DV17 <- data_DV17[,c(4,6,7,12,13,15,16)] 
data_DV17 <- data_DV17[data_DV17$PARCELLE_ID!="5272",]
data_DV17NA <- data_DV17[complete.cases(data_DV17), ]
str(data_DV17NA) #1404 obs. of 7 variables

data_DV17_1 <- data_DV17NA[data_DV17NA$PARCELLE_ID=="1570",]
str(data_DV17_1) #701 obs
data_DV17_2 <- data_DV17NA[data_DV17NA$PARCELLE_ID=="3450",]
str(data_DV17_2) #131 obs
data_DV17_3 <- data_DV17NA[data_DV17NA$PARCELLE_ID=="5042",]
str(data_DV17_3) #512 obs

#lisData_DV17 <- list(data_DV17_1,data_DV17_2,data_DV17_3)



#******************************************************
# Stand data description (Table 1 ; Figure 2)
#******************************************************

### Characterization of substrates
### Description of metapopulation dynamic within each stand (Current & Past occupancy, colonization, extinction)
### Plot "potential unoccupied" (0-0), "survival"(1-1), "potential colonization" (0-1), "extinction"(1-0) within each stand
### Regional connectivity of each stand 

#--------------------------------------
# Characterization of substrates 
#--------------------------------------

### Number of substrates
nb.sub.p1570 <- nrow(data_DV17_1) 
nb.sub.p3450 <- nrow(data_DV17_2) 
nb.sub.p5042 <- nrow(data_DV17_3) 

### Diameter of substrat (mean+-sd)
diam.mean.p1570 <- mean(data_DV17_1$DIAMETRE)
diam.sd.p1570 <- sd(data_DV17_1$DIAMETRE)
diam.mean.p3450 <- mean(data_DV17_2$DIAMETRE)
diam.sd.p3450 <- sd(data_DV17_2$DIAMETRE)
diam.mean.p5042 <- mean(data_DV17_3$DIAMETRE)
diam.sd.p5042 <- sd(data_DV17_3$DIAMETRE)

### Distribution of tree species
plot_ESS.CATEG_p1570 <- ggplot(data_DV17_1, aes(ESS_CATEG)) + geom_bar(fill = "#0073C2FF") 
plot_ESS.CATEG_p1570
nb.ess.aut.p1570 <- data.frame(table(data_DV17_1$ESS_CATEG))[1,2]
nb.ess.het.p1570 <- data.frame(table(data_DV17_1$ESS_CATEG))[2,2]
plot_ESS.CATEG_p3450 <- ggplot(data_DV17_2, aes(ESS_CATEG)) + geom_bar(fill = "#0073C2FF") 
plot_ESS.CATEG_p3450
nb.ess.aut.p3450 <- data.frame(table(data_DV17_2$ESS_CATEG))[1,2]
nb.ess.het.p3450 <- data.frame(table(data_DV17_2$ESS_CATEG))[2,2]
plot_ESS.CATEG_p5042 <- ggplot(data_DV17_3, aes(ESS_CATEG)) + geom_bar(fill = "#0073C2FF") 
plot_ESS.CATEG_p5042
nb.ess.aut.p5042 <- data.frame(table(data_DV17_3$ESS_CATEG))[1,2]
nb.ess.het.p5042 <- data.frame(table(data_DV17_3$ESS_CATEG))[2,2]

####### Compilation table 
sub.p1570 <- rbind(nb.sub.p1570,diam.mean.p1570,diam.sd.p1570,nb.ess.aut.p1570,nb.ess.het.p1570)
sub.p3450 <- rbind(nb.sub.p3450,diam.mean.p3450,diam.sd.p3450,nb.ess.aut.p3450,nb.ess.het.p3450)
sub.p5042 <- rbind(nb.sub.p5042,diam.mean.p5042,diam.sd.p5042,nb.ess.aut.p5042,nb.ess.het.p5042)
Caract.Sub <- data.frame(cbind(sub.p1570,sub.p3450,sub.p5042))
colnames(Caract.Sub) <- c("p1570","p3450","p5042")
rownames(Caract.Sub) <- c("nb.sub","diam.mean","diam.sd","nb.ess.aut","nb.ess.het")
View(Caract.Sub)


#------------------------------------------------------------------------------------------------------------------
# 1.2 Description of metapopulation dynamic within each stand (Current & Past occupancy, colonization, extinction)
#------------------------------------------------------------------------------------------------------------------
# /!\ We're talking about "potential colonization" here, given the uncertainty about past unoccupied substrates !

### Current & Past occupancy
nb.occ14.p1570 <- nrow(data_DV17_1[data_DV17_1$DETECT14==1,])
tx.occ14.p1570 <- nb.occ14.p1570/nb.sub.p1570
nb.occ17.p1570 <- nrow(data_DV17_1[data_DV17_1$DETECT17==1,])
tx.occ17.p1570 <- nb.occ17.p1570/nb.sub.p1570

nb.occ14.p3450 <- nrow(data_DV17_2[data_DV17_2$DETECT14==1,])
tx.occ14.p3450 <- nb.occ14.p3450/nb.sub.p3450
nb.occ17.p3450 <- nrow(data_DV17_2[data_DV17_2$DETECT17==1,])
tx.occ17.p3450 <- nb.occ17.p3450/nb.sub.p3450

nb.occ14.p5042 <- nrow(data_DV17_3[data_DV17_3$DETECT14==1,])
tx.occ14.p5042 <- nb.occ14.p5042/nb.sub.p5042
nb.occ17.p5042 <- nrow(data_DV17_3[data_DV17_3$DETECT17==1,])
tx.occ17.p5042 <- nb.occ17.p5042/nb.sub.p5042

### Potential colonization
nb.sub.unocc.p1570 <- nrow(data_DV17_1[data_DV17_1$DETECT14==0,])
nb.colo.p1570 <- (nrow(data_DV17_1[data_DV17_1$DETECT14==0 & data_DV17_1$DETECT17==1,]))
tx.colo.p1570 <- nb.colo.p1570/nb.sub.nocc.p1570 

nb.sub.unocc.p3450 <- nrow(data_DV17_2[data_DV17_2$DETECT14==0,])
nb.colo.p3450 <- (nrow(data_DV17_2[data_DV17_2$DETECT14==0 & data_DV17_2$DETECT17==1,]))
tx.colo.p3450 <- nb.colo.p3450/nb.sub.nocc.p3450 

nb.sub.unocc.p5042 <- nrow(data_DV17_3[data_DV17_3$DETECT14==0,])
nb.colo.p5042 <- (nrow(data_DV17_3[data_DV17_3$DETECT14==0 & data_DV17_3$DETECT17==1,]))
tx.colo.p5042 <- nb.colo.p5042/nb.sub.nocc.p5042 

### Extinction
nb.sub.occ.p1570 <- nrow(data_DV17_1[data_DV17_1$DETECT14==1,])
nb.ext.p1570 <- (nrow(data_DV17_1[data_DV17_1$DETECT14==1 & data_DV17_1$DETECT17==0,]))
tx.ext.p1570 <- nb.ext.p1570/nb.sub.occ.p1570 

nb.sub.occ.p3450 <- nrow(data_DV17_2[data_DV17_2$DETECT14==1,])
nb.ext.p3450 <- (nrow(data_DV17_2[data_DV17_2$DETECT14==1 & data_DV17_2$DETECT17==0,]))
tx.ext.p3450 <- nb.ext.p3450/nb.sub.occ.p3450 

nb.sub.occ.p5042 <- nrow(data_DV17_3[data_DV17_3$DETECT14==1,])
nb.ext.p5042 <- (nrow(data_DV17_3[data_DV17_3$DETECT14==1 & data_DV17_3$DETECT17==0,]))
tx.ext.p5042 <- nb.ext.p5042/nb.sub.occ.p5042 


####### Compilation table 
Caract.Metapop.p1570 <- rbind(nb.occ14.p1570,tx.occ14.p1570,nb.occ17.p1570,tx.occ17.p1570,nb.sub.unocc.p1570,nb.colo.p1570,tx.colo.p1570,nb.sub.occ.p1570,nb.ext.p1570,tx.ext.p1570)
Caract.Metapop.p3450 <- rbind(nb.occ14.p3450,tx.occ14.p3450,nb.occ17.p3450,tx.occ17.p3450,nb.sub.unocc.p3450,nb.colo.p3450,tx.colo.p3450,nb.sub.occ.p3450,nb.ext.p3450,tx.ext.p3450)
Caract.Metapop.p5042 <- rbind(nb.occ14.p5042,tx.occ14.p5042,nb.occ17.p5042,tx.occ17.p5042,nb.sub.unocc.p5042,nb.colo.p5042,tx.colo.p5042,nb.sub.occ.p5042,nb.ext.p5042,tx.ext.p5042)
Caract.Metapop <- data.frame(cbind(Caract.Metapop.p1570, Caract.Metapop.p3450, Caract.Metapop.p5042))
colnames(Caract.Metapop) <- c("p1570","p3450","p5042")
rownames(Caract.Metapop) <- c("nb.occ14","tx.occ14","nb.occ17","tx.occ17","nb.sub.unocc","nb.colo","tx.colo","nb.sub.occ","nb.ext","tx.ext")
View(Caract.Metapop)

#___________________________________________________________________
# Quantify extinction events (1-0) in relation to survival (1-1) 

dataext17  <- data_DV17NA[data_DV17NA$DETECT14==1 & data_DV17NA$DETECT17==0,]
dataocc1417  <- data_DV17NA[data_DV17NA$DETECT14==1 & data_DV17NA$DETECT17==1,]
prop.ext <- (nrow(dataext17)/nrow(dataocc1417))*100
prop.ext
#___________________________________________________________________


#------------------------------------------------------------------------------------------------------------------------------
# 1.3 Plot "potential unoccupied" (0-0), "survival"(1-1), "potential colonization" (0-1), "extinction"(1-0) within each stand
#------------------------------------------------------------------------------------------------------------------------------

#p1570
data_DV17_1$ETAT <-
  ifelse(data_DV17_1$DETECT14==0 & data_DV17_1$DETECT17==1 ,"COL",
         ifelse(data_DV17_1$DETECT14==0 & data_DV17_1$DETECT17==0,"NOC",
                ifelse(data_DV17_1$DETECT14==1 & data_DV17_1$DETECT17==1,"OCC","EXT"
                )))
str(data_DV17_1)

setwd("C:/Users/HP050/Documents/Metapop_colo_Octobre2017_4/SIG/BRYO_metapop2017/Vecteur/UGS_ech_Oct17/Parcelles_ech/polygones_L93/polygone_modif")	
polygon1570mod <- readShapeSpatial("p1570u3_echDviride_Oct2017_modif.shp")
datapolyg1570mod <- fortify(polygon1570mod)

plot.DV.p1570.ptpg.mod <- ggplot(data=datapolyg1570mod,aes(long,lat)) + 
  geom_polygon(colour="black", fill="white") + #grey90 
  geom_point(data=data_DV17_1,aes(x=LONGITUDE,y=LATITUDE, shape=ETAT, colour=ETAT),size=2)+
  scale_shape_manual(values = c(15,7,0,15))+ 
  scale_colour_manual(values = c("Red","Red","Black","Black"))+ 
  xlim(692500,NA)+
  ylim(NA,6922700)+
  ggtitle("Distribution Dicranum viride - peuplement 1570")+
  xlab("longitude")+
  ylab("latitude")+
  theme_classic() +
  theme(axis.line=element_blank(),panel.border = element_rect(colour = "black", fill=NA, size=2)) #theme(panel.background = element_rect(colour = "black", size=1))#multiplot(plot.DV.p1570.ptpg,plot.DV.p1570.ptpg.mod)
plot.DV.p1570.ptpg.mod

#p3450
data_DV17_2$ETAT <-
  ifelse(data_DV17_2$DETECT14==0 & data_DV17_2$DETECT17==1 ,"COL",
         ifelse(data_DV17_2$DETECT14==0 & data_DV17_2$DETECT17==0,"NOC",
                ifelse(data_DV17_2$DETECT14==1 & data_DV17_2$DETECT17==1,"OCC","EXT"
                )))
str(data_DV17_2)

setwd("C:/Users/HP050/Documents/Metapop_colo_Octobre2017_4/SIG/BRYO_metapop2017/Vecteur/UGS_ech_Oct17/Parcelles_ech/polygones_L93/polygone_modif")	
polygon3450mod <- readShapeSpatial("p3450c_echDviride_Oct2017_modif.shp")
datapolyg3450mod <- fortify(polygon3450mod)

plot.DV.p3450.ptpg.mod <- ggplot(data=datapolyg3450mod,aes(long,lat)) + 
  geom_polygon(colour="black", fill="white") + #grey90 
  geom_point(data=data_DV17_2,aes(x=LONGITUDE,y=LATITUDE, shape=ETAT, colour=ETAT),size=2)+
  scale_shape_manual(values = c(15,0,15))+ 
  scale_colour_manual(values = c("Red","Black","Black"))+ 
  xlim(NA,691500)+
  ylim(NA,6919400)+
  ggtitle("Distribution Dicranum viride - peuplement 3450")+
  xlab("longitude")+
  ylab("latitude")+
  theme_classic() +
  theme(axis.line=element_blank(),panel.border = element_rect(colour = "black", fill=NA, size=2)) #theme(panel.background = element_rect(colour = "black", size=1))#multiplot(plot.DV.p1570.ptpg,plot.DV.p1570.ptpg.mod)
plot.DV.p3450.ptpg.mod

#p5042
data_DV17_3$ETAT <-
  ifelse(data_DV17_3$DETECT14==0 & data_DV17_3$DETECT17==1 ,"COL",
         ifelse(data_DV17_3$DETECT14==0 & data_DV17_3$DETECT17==0,"NOC",
                ifelse(data_DV17_3$DETECT14==1 & data_DV17_3$DETECT17==1,"OCC","EXT"
                )))
str(data_DV17_3)

setwd("C:/Users/HP050/Documents/Metapop_colo_Octobre2017_4/SIG/BRYO_metapop2017/Vecteur/UGS_ech_Oct17/Parcelles_ech/polygones_L93/polygone_modif")	
polygon5042mod <- readShapeSpatial("p5042d_echDviride_Oct2017_modif.shp")
datapolyg5042mod <- fortify(polygon5042mod)

plot.DV.p5042.ptpg.mod <- ggplot(data=datapolyg5042mod,aes(long,lat)) + 
  geom_polygon(colour="black", fill="white") + #grey90 
  geom_point(data=data_DV17_3,aes(x=LONGITUDE,y=LATITUDE, shape=ETAT, colour=ETAT),size=2)+
  scale_shape_manual(values = c(15,7,0,15))+ 
  scale_colour_manual(values = c("Red","Red","Black","Black"))+ 
  #xlim(NA,NA)+
  #ylim(NA,6920400)+
  ggtitle("Distribution Dicranum viride - peuplement 5042")+
  xlab("longitude")+
  ylab("latitude")+
  theme_classic() +
  theme(axis.line=element_blank(),panel.border = element_rect(colour = "black", fill=NA, size=2)) #theme(panel.background = element_rect(colour = "black", size=1))#multiplot(plot.DV.p1570.ptpg,plot.DV.p1570.ptpg.mod)
plot.DV.p5042.ptpg.mod


#complilation plot
source("http://peterhaschke.com/Code/multiplot.R")
multiplot(plot.DV.p1570.ptpg.mod,plot.DV.p3450.ptpg.mod,plot.DV.p5042.ptpg.mod,cols=3)


#---------------------------------------------
# 1.4 Regional connectivity of each stand 
#---------------------------------------------

### See the R-script "Supplementary information and code for the study "Disentangling fine- and large-scale colonization processes in metapopulation dynamics: a case study on a threatened epiphytic bryophyte"
### Section "Scale parameter of regional colonization process"


###############################################################
### END
###############################################################