#####################data analysis##########################################
##Note:The analyses not included were conducted using SPSS or Origin.


##1.Random Forest model
####1.1Autocorrelation Analysis
library(corrplot)
library(pheatmap)
df.zyq<-read.csv("C:/Users/tsw/Desktop/BaiduSyncdisk/data1.csv")
attach(df.zyq)
cor(df.zyq)
pheatmap::pheatmap(cor(df.zyq))
corrplot(cor(df.zyq))
corrplot(cor(df.zyq),method = "pie")
corrplot(cor(df.zyq),method = "color",addCoef.col = "grey")
col=colorRampPalette(c("navy","white","firebrick3"))
corrplot(cor(df.zyq),type = "upper",col = col(10),tl.pos = "d")
corrplot(cor(df.zyq),add = TRUE,type = "lower",method="number"
         ,col ="black",tl.pos =  "d",cl.pos = "n")

#1.2Random Forest model

#bm
library(randomForest)
setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.train<-read.csv("bm.csv")
set.seed(6666)
df.train$X<-factor(df.train$X)
fit.forest<-randomForest(X~.,data=df.train,
                         na.action = na.pass,
                         importance=TRUE)
fit.forest
importance(fit.forest,type=2)
hist(treesize(fit.forest))
plot(fit.forest)
dev.new(
  width=16,
  height=8,
  noRStudioGD = TRUE
)
varImpPlot(fit.forest, main="variable importance")
partialPlot(fit.forest, pred.data=df.train, x.var=HA,which.class=1,
            plot = TRUE, add = FALSE,
            rug = FALSE, xlab=deparse(substitute(HA)))

df.zy<-importance(fit.forest,type=2)
write.csv(df.zy,file = "1.csv")

library(pdp)
library(randomForest)
df.train<-read.csv("bm.csv")
rf_model <- randomForest(X ~ ., data = df.train)
partial(rf_model, pred.var = "HA",plot = TRUE, plot.type = "prob")

pdp_res <- list()
for (var in colnames(df.train)) {
  pdp_res[[var]] <- partial(rf_model, pred.var = var)
}


write.csv(pdp_res[["HA"]], file = "my_file.csv")
write.csv(pdp_res[["T"]], file = "my_file.csv")
write.csv(pdp_res[["HE"]], file = "my_file.csv")
write.csv(pdp_res[["VC"]], file = "my_file.csv")
write.csv(pdp_res[["S"]], file = "my_file.csv")
write.csv(pdp_res[["DRO"]], file = "my_file.csv")
write.csv(pdp_res[["DW"]], file = "my_file.csv")
write.csv(pdp_res[["DRE"]], file = "my_file.csv")   



#bf
library(randomForest)
setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.train<-read.csv("bf.csv")
set.seed(6666)
df.train$X<-factor(df.train$X)
fit.forest<-randomForest(X~.,data=df.train,
                         na.action = na.pass,
                         importance=TRUE)
fit.forest
importance(fit.forest,type=2)
hist(treesize(fit.forest))
plot(fit.forest)
dev.new(
  width=16,
  height=8,
  noRStudioGD = TRUE
)
varImpPlot(fit.forest, main="variable importance")
partialPlot(fit.forest, pred.data=df.train, x.var=HA,which.class=1,
            plot = TRUE, add = FALSE,
            rug = FALSE, xlab=deparse(substitute(HA)))

df.zy<-importance(fit.forest,type=2)
write.csv(df.zy,file = "1.csv")

library(pdp)
library(randomForest)
df.train<-read.csv("bf.csv")
rf_model <- randomForest(X ~ ., data = df.train)
partial(rf_model, pred.var = "HA",plot = TRUE, plot.type = "prob")

pdp_res <- list()
for (var in colnames(df.train)) {
  pdp_res[[var]] <- partial(rf_model, pred.var = var)
}

write.csv(pdp_res[["HA"]], file = "my_file.csv")
write.csv(pdp_res[["T"]], file = "my_file.csv")
write.csv(pdp_res[["HE"]], file = "my_file.csv")
write.csv(pdp_res[["VC"]], file = "my_file.csv")
write.csv(pdp_res[["S"]], file = "my_file.csv")
write.csv(pdp_res[["DRO"]], file = "my_file.csv")
write.csv(pdp_res[["DW"]], file = "my_file.csv")
write.csv(pdp_res[["DRE"]], file = "my_file.csv")   


##nbm
library(randomForest)
setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.train<-read.csv("nbm.csv")
set.seed(6666)
df.train$X<-factor(df.train$X)
fit.forest<-randomForest(X~.,data=df.train,
                         na.action = na.pass,
                         importance=TRUE)
fit.forest
importance(fit.forest,type=2)
hist(treesize(fit.forest))
plot(fit.forest)
dev.new(
  width=16,
  height=8,
  noRStudioGD = TRUE
)
varImpPlot(fit.forest, main="variable importance")
partialPlot(fit.forest, pred.data=df.train, x.var=HA,which.class=1,
            plot = TRUE, add = FALSE,
            rug = FALSE, xlab=deparse(substitute(HA)))

df.zy<-importance(fit.forest,type=2)
write.csv(df.zy,file = "1.csv")

library(pdp)
library(randomForest)
df.train<-read.csv("nbm.csv")
rf_model <- randomForest(X ~ ., data = df.train)
partial(rf_model, pred.var = "HA",plot = TRUE, plot.type = "prob")

pdp_res <- list()
for (var in colnames(df.train)) {
  pdp_res[[var]] <- partial(rf_model, pred.var = var)
}


write.csv(pdp_res[["HA"]], file = "my_file.csv")
write.csv(pdp_res[["T"]], file = "my_file.csv")
write.csv(pdp_res[["HE"]], file = "my_file.csv")
write.csv(pdp_res[["VC"]], file = "my_file.csv")
write.csv(pdp_res[["S"]], file = "my_file.csv")
write.csv(pdp_res[["DRO"]], file = "my_file.csv")
write.csv(pdp_res[["DW"]], file = "my_file.csv")
write.csv(pdp_res[["DRE"]], file = "my_file.csv") 

#nbf
library(randomForest)
setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.train<-read.csv("nbf.csv")
set.seed(6666)
df.train$X<-factor(df.train$X)
fit.forest<-randomForest(X~.,data=df.train,
                         na.action = na.pass,
                         importance=TRUE)
fit.forest
importance(fit.forest,type=2)
hist(treesize(fit.forest))
plot(fit.forest)
dev.new(
  width=16,
  height=8,
  noRStudioGD = TRUE
)
varImpPlot(fit.forest, main="variable importance")
partialPlot(fit.forest, pred.data=df.train, x.var=HA,which.class=1,
            plot = TRUE, add = FALSE,
            rug = FALSE, xlab=deparse(substitute(HA)))

df.zy<-importance(fit.forest,type=2)
write.csv(df.zy,file = "1.csv")

library(pdp)
library(randomForest)
df.train<-read.csv("nbf.csv")
rf_model <- randomForest(X ~ ., data = df.train)
partial(rf_model, pred.var = "HA",plot = TRUE, plot.type = "prob")

pdp_res <- list()
for (var in colnames(df.train)) {
  pdp_res[[var]] <- partial(rf_model, pred.var = var)
}


write.csv(pdp_res[["HA"]], file = "my_file.csv")
write.csv(pdp_res[["T"]], file = "my_file.csv")
write.csv(pdp_res[["HE"]], file = "my_file.csv")
write.csv(pdp_res[["VC"]], file = "my_file.csv")
write.csv(pdp_res[["S"]], file = "my_file.csv")
write.csv(pdp_res[["DRO"]], file = "my_file.csv")
write.csv(pdp_res[["DW"]], file = "my_file.csv")
write.csv(pdp_res[["DRE"]], file = "my_file.csv")  


#2. niche overlap

##2.1  observed value
install.packages("spaa")
library(spaa)

setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.zyq<-read.csv("1.csv")


data2mat(data = df.zyq)
spmatrix <- data2mat(df.zyq)



df.mat<-read.csv("C:/Users/tsw/Desktop/BaiduSyncdisk/2.csv")
niche.width(df.mat,method = "levins")


niche.overlap(df.mat,method = "levins")

niche.overlap.boot(df.mat,method = "levins",times = 999, quant = c(0.025, 0.975))

##2.2 randomly resampled

x=1:92
y=sample(x=x,size=70)
write.csv(y,file = "C:/Users/tsw/Desktop/BaiduSyncdisk/t.CSV",)
Z=sample(1:2, 92, replace = TRUE )
write.csv(Z,file = "C:/Users/tsw/Desktop/BaiduSyncdisk/t2.CSV",)

library(spaa)

setwd("C:/Users/tsw/Desktop/BaiduSyncdisk/data analysis")
df.zyq<-read.csv("1.CSV")


data2mat(data = df.zyq)
spmatrix <- data2mat(df.zyq)



df.mat<-read.csv("2.csv")
niche.width(df.mat,method = "levins")


niche.overlap(df.mat,method = "levins")

niche.overlap.boot(df.mat,method = "levins",times = 999, quant = c(0.025, 0.975))