##tittle: The driver of dispersal model #########################
##date: 20240102#################################################
##Authors Zhiwen Gao et al, 2024##############################################

###package been used in this analysis#####

library(PerformanceAnalytics)
library(lme4)
library(lmerTest)
library(performance)
library(dplyr)
library(rsq)
library(MuMIn)
require(ggplot2)
require(pscl)
require(boot)
library(sjPlot)
library(glmmTMB)
library(GLMMadaptive)
library(partR2)
library(DHARMa)

##data preparing##############################

dis.patch <- read.csv("231219patch scale.csv")

str(dis.patch)

hist(dis.patch$City_size)
hist(log(dis.patch$City_size))
dis.patch$City_size<-log(dis.patch$City_size)

hist(dis.patch$Area)
hist(log(dis.patch$Area))
dis.patch$Area<-log(dis.patch$Area)


####avoid over dispersion, so area and city size been loged in this case
dis.patch$Area <- log(dis.patch$Area)
dis.patch$City_size <- log(dis.patch$City_size)


str(dis.patch)

dis.patch[,c(8:ncol(dis.patch))]<-lapply(dis.patch[,c(8:ncol(dis.patch))],
                                         FUN = function(y){as.numeric(y)})

str(dis.patch)


##standardize independent variables##########################################

dis.patch[,-c(1:15)] <- scale(dis.patch[,-c(1:15)])


##check multicollinearity#######################################################
#####To avoid collinearity between variables, if the correlation coefficient between 
#####independent variables was bigger than 0.7, only one of the variables was retained in the model.

chart.Correlation(dis.patch[,c(12:ncol(dis.patch))],method= "spearman")###removed perimeter, population.city, elevation and CMT


##Build the full models by 3 types of random effect#########################


options(na.action = na.omit)

global.Auto1<-glmer.nb(Auto.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       na.action = na.omit,
                       data = dis.patch,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+5)))
summary(global.Auto1)
check_model(global.Auto1)
plot(check_outliers(global.Auto1))
check_autocorrelation(global.Auto1)

simres.nb = simulateResiduals(global.Auto1)
plot(simres.nb)


dis.patch1 <- dis.patch[-c(773,127,294,794,300),]#remove outliers

global.Auto2<-glmer.nb(Auto.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       #na.action = na.omit,
                       data = dis.patch1,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+6)))##after remove outliers, global.Auto2 is much better

AICc(global.Auto1,global.Auto2)

check_model(global.Auto2)
plot(check_outliers(global.Auto2))

simres.nb = simulateResiduals(global.Auto2)
plot(simres.nb)

dis.patch2 <- dis.patch1[-c(480,589,162,320),]#remove outliers again

global.Auto3<-glmer.nb(Auto.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       #na.action = na.omit,
                       data = dis.patch2,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+6)))

AICc(global.Auto1,global.Auto2,global.Auto3)

check_model(global.Auto3)
plot(check_outliers(global.Auto3))

simres.nb = simulateResiduals(global.Auto3)
plot(simres.nb)

VarCorr(global.Auto3)




##model selection################################
##################################################
summary(global.Auto3)


options(na.action = na.fail)

se_Auto<- dredge(global.Auto3, evaluate = TRUE,  rank = "AICc",REML=F)



##model average#####################################

avg_Auto<-model.avg(se_Auto,subset = delta<2,revised.var = TRUE,fit = T)

summary(avg_Auto)
summary(global.Auto3)
best.Auto<-global.Auto3


plot_model(best.Auto,type = "eff", terms = c("Area","City_size"),show.data = T)
plot_model(best.Auto,type = "eff", terms = c("UR.city"),show.data = T)
plot_model(best.Auto,type = "eff", terms = c("MAT"),show.data = T)
plot_model(best.Auto,type = "eff", terms = c("MAP"),show.data = T)

plot_model(best.Auto,type = "re") ####for random effect

plot_model(best.Auto,type = "std2")####for standardized coefficient
plot_model(best.Auto,type = "pred", show.data = T) #####for prediction
plot_model(best.Auto,type = "resid")


##Variance composition####################################################
#################################
r.Auto=r.squaredGLMM(best.Auto)


R2m.Auto = r.Auto[1,1] 
R2c.Auto = r.Auto[1,2]
R2.Auto <- c(R2m.Auto,R2c.Auto, r.Auto[1,2]-r.Auto[1,1],1-r.Auto[1,2])
names(R2.Auto) <- c( "R2m.Auto","R2c.Auto", "Random", "Unexplained")

R2.Auto

d1.Auto=abs(summary(avg_Auto)$coefmat.full[-1,1])


####calculate the relative contribution of every fixed effects
####estimate/(sum(estimate))*marginal R2#######################


d2.Auto=d1.Auto/sum(d1.Auto)*R2m.Auto

Auto.R2=as.data.frame(c(R2.Auto,d2.Auto))

Auto.coef=as.data.frame(summary(avg_Auto)$coefmat.full)
Auto.R2$x=rownames(Auto.R2)
Auto.coef$x <- rownames(Auto.coef)

Auto.result=full_join(Auto.coef,Auto.R2)



##Anemochory#####################################################################

options(na.action = na.fail)

global.Wind<-glmer(Wind.s~Area+LSI+BD+Sealed25+
                     UR.city+
                     City_size+
                     MAP+MAT+(1|SiteID),
                   data = dis.patch,
                   na.action = na.omit,
                   control = glmerControl(optimizer = "bobyqa"),
                   poisson(link ="log"))


global.Wind1<-glmer(Wind.s~Area+LSI+BD+Sealed25+
                      UR.city+
                      City_size+
                      MAP+MAT+(1|CityID),
                    na.action = na.omit,
                    data = dis.patch,
                    control = glmerControl(optimizer = "bobyqa"),
                    poisson(link ="log"))

global.Wind2<-glmer(Wind.s~Area+LSI+BD+Sealed25+
                      UR.city+
                      City_size+
                      MAP+MAT+(1|CityID/SiteID),
                    data = dis.patch,
                    na.action = na.omit,
                    control = glmerControl(optimizer = "bobyqa"),
                    poisson(link ="log"))


global.Wind3<-glmer.nb(Wind.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       na.action = na.omit,
                       data = dis.patch,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+5)))


AICc(global.Wind3,global.Wind,global.Wind1,global.Wind2)

dis.wind1 <- dis.patch[-c(773,863,849,462,211),]

global.Wind4<-glmer.nb(Wind.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       na.action = na.omit,
                       data = dis.wind1,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+5)))

dis.wind2<-dis.wind1[-c(398,1,163,825,453),]

global.Wind5<-glmer.nb(Wind.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       na.action = na.omit,
                       data = dis.wind2,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+5)))

dis.wind3<-dis.wind2[-c(32,331,305,424,1),]


options(na.action = na.omit)

global.Wind6<-glmer.nb(Wind.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       #na.action = na.omit,
                       data = dis.wind3,
                       control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                              optCtrl=list(maxfun=1e+6)))



AICc(global.Wind3,global.Wind4,global.Wind5,
     global.Wind6)

options(na.action = na.fail)

se_Wind<- dredge(global.Wind6, evaluate = TRUE,  rank = "AICc",REML=F)

avg_Wind<-model.avg(se_Wind,subset = delta < 2, revised.var = TRUE,fit = T)

summary(avg_Wind)

best.wind<-glmer.nb(Wind.s~Area+LSI+BD+Sealed25+
                      City_size+
                      MAP+MAT+(1|CityID/SiteID),
                    na.action = na.omit,
                    data = dis.wind3,
                    control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                           optCtrl=list(maxfun=1e+6)))



r.Wind=r.squaredGLMM(best.wind)


R2m.Wind = r.Wind[1,1] 
R2c.Wind = r.Wind[1,2]
R2.Wind <- c(R2m.Wind,R2c.Wind, r.Wind[1,2]-r.Wind[1,1],1-r.Wind[1,2])
names(R2.Wind) <- c( "R2m.Wind","R2c.Wind", "Random", "Unexplained")
R2.Wind

d1.Wind=abs(summary(avg_Wind)$coefmat.full[-1,1])

d2.Wind=d1.Wind/sum(d1.Wind)*R2m.Wind

Wind.R2=as.data.frame(c(R2.Wind,d2.Wind))

Wind.coef=as.data.frame(summary(avg_Wind)$coefmat.full)
Wind.R2$x=rownames(Wind.R2)
Wind.coef$x <- rownames(Wind.coef)


Wind.result=full_join(Wind.coef,Wind.R2)


##Zoochory##################################################
##############################################################

options(na.action = na.omit)

global.Animal<-glmer(Animal.s~Area+LSI+BD+Sealed25+
                       UR.city+
                       City_size+
                       MAP+MAT+(1|SiteID),
                     data = dis.patch,
                     control = glmerControl(optimizer = "bobyqa"),
                     poisson(link ="log"))


global.Animal1<-glmer(Animal.s~Area+LSI+BD+Sealed25+
                        UR.city+
                        City_size+
                        MAP+MAT+(1|CityID),
                      data = dis.patch,
                      control = glmerControl(optimizer = "bobyqa"),
                      poisson(link ="log"))



global.Animal2<-glmer(Animal.s~Area+LSI+BD+Sealed25+
                        UR.city+
                        City_size+
                        MAP+MAT+(1|CityID/SiteID),
                      data = dis.patch,
                      control = glmerControl(optimizer = "bobyqa"),
                      poisson(link ="log"))

global.Animal3<-glmer.nb(Animal.s~Area+LSI+BD+Sealed25+
                           UR.city+
                           City_size+
                           MAP+MAT+(1|CityID/SiteID),
                         #na.action = na.omit,
                         data = dis.patch,
                         control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                                optCtrl=list(maxfun=1e+6)))

anova(global.Animal,global.Animal1,global.Animal2)

AICc(global.Animal3)

summary(global.Animal3)


VarCorr(global.Animal3)


simres.nb = simulateResiduals(global.Animal3)
plot(simres.nb)

options(na.action = na.fail)

se_Animal<- dredge(global.Animal3, evaluate = TRUE,  rank = "AICc",REML=F)


avg_Animal<-model.avg(se_Animal,subset = delta<2,revised.var = TRUE,fit = T)

summary(avg_Animal)
best.Animal<-global.Animal3





r.Animal=r.squaredGLMM(best.Animal)


R2m.Animal = r.Animal[1,1] 
R2c.Animal = r.Animal[1,2]
R2_Animal <- c(R2m.Animal,R2c.Animal, r.Animal[1,2]-r.Animal[1,1],1-r.Animal[1,2])
names(R2_Animal) <- c( "R2m.Animal","R2c.Animal", "Random", "Unexplained")
R2_Animal

d1.Animal=abs(summary(avg_Animal)$coefmat.full[-1,1])

d2.Animal=d1.Animal/sum(d1.Animal)*R2m.Animal

Animal.R2=as.data.frame(c(R2_Animal,d2.Animal))

Animal.coef=as.data.frame(summary(avg_Animal)$coefmat.full)

Animal.R2$x=rownames(Animal.R2)

Animal.coef$x <- rownames(Animal.coef)

Animal.result=full_join(Animal.coef,Animal.R2)

##Hydrochory#################################################
########################################################################


options(na.action = na.omit)

global.Water<-glmmTMB(Water.s~Area+LSI+BD+Sealed25+
                        UR.city+
                        City_size+
                        MAP+MAT+(1|SiteID),
                      data = dis.patch,
                      ziformula=~1,
                      family = poisson(link = "log"))



global.Water1<-glmmTMB(Water.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID),
                       data = dis.patch,
                       ziformula=~1,
                       family = poisson(link = "log"))



global.Water2<-glmmTMB(Water.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       data = dis.patch,
                       ziformula=~1,
                       family = poisson(link = "log"))


global.Water3<-glmmTMB(Water.s~Area+LSI+BD+Sealed25+
                         UR.city+
                         City_size+
                         MAP+MAT+(1|CityID/SiteID),
                       data = dis.patch,
                       ziformula=~1,
                       family = nbinom2(link = "log"))




global.Water4<-glmer.nb(Water.s~Area+LSI+BD+Sealed25+
                          UR.city+
                          City_size+
                          MAP+MAT+(1|CityID/SiteID),
                        #na.action = na.omit,
                        data = dis.patch,
                        control = glmerControl(optimizer = "bobyqa", tol = 1e-5,
                                               optCtrl=list(maxfun=1e+5)))




AIC(global.Water,global.Water1,global.Water2,global.Water3,global.Water4)

summary(global.Water4)


simres.nb = simulateResiduals(global.Water4)
plot(simres.nb)

options(na.action = na.fail)

se_Water<- dredge(global.Water4, evaluate = TRUE,  rank = "AICc",REML=F)

avg_Water<-model.avg(se_Water,subset = delta<2,revised.var = TRUE,fit = T)

summary(avg_Water)


best.Water<-global.Water4



r.Water=r.squaredGLMM(best.Water)


R2m.Water = r.Water[1,1] 
R2c.Water = r.Water[1,2]
R2_Water <- c(R2m.Water,R2c.Water, r.Water[1,2]-r.Water[1,1],1-r.Water[1,2])
names(R2_Water) <- c( "R2m.Water","R2c.Water", "Random", "Unexplained")
R2_Water

d1.Water=abs(summary(avg_Water)$coefficients[1,-1])

d2.Water=d1.Water/sum(d1.Water)*R2m.Water

Water.R2=as.data.frame(c(R2_Water,d2.Water))

Water.coef=as.data.frame(summary(avg_Water)$coefmat.full)
Water.R2$x=rownames(Water.R2)
Water.coef$x <- rownames(Water.coef)


Water.sw<-as.data.frame(summary(avg_Water)$sw)

Water.sw$x <- rownames(Water.sw)

Water.result=full_join(Water.coef,Water.R2)

Water.result=full_join(Water.result,Water.sw)


Auto.result$dis<-c("Autochory")
Animal.result$dis<-c("Zoochory")

Wind.result$dis<-c("Anemochory")

Water.result$dis<-c("Hydrochory")



##results############################

Auto.result
Animal.result
Wind.result
Water.result