How to reproduce the MaxEnt Response Curves with R

package_list=c("ade4","adehabitatHR","ggplot2","maptools","readr","reshape2","reshape","spatstat",
               "spatialEco","spatstat.utils","raster","ggpubr","gridExtra","cowplot","ggThemeAssist",
               "ecospat","shapefiles","png","ggthemes")

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE, repos = "http://cran.us.r-project.org")
  sapply(pkg, require, character.only = TRUE)
}

ipak(package_list)

Maxent_responses_F2 <- list.files(path="D:/Luc/Articles/Mike/Results/MaxentGambiaCrossValidate_New_Road_Livestock/plots", 
                                  pattern=".dat$", full.names=TRUE)


Maxent_stack = lapply(Maxent_responses_F2, read.csv)

typeof(Maxent_stack)

Var_Num=length(Maxent_responses_F2)/11
Variables_order=seq(1,Var_Num,1)


dim_empty=dim(data.frame(Maxent_stack[[1]]))
mat_empty = matrix(ncol = 0, nrow = dim_empty[1])
df_empty=df=data.frame(mat_empty) 



selection=seq(1,1001, by=2)

VariableSelection=c(4,5,7,8,10,11,13,14,16,17,19,20,22,23,25,26,28,29,31,32)


for (i in Variables_order) {

  assign(paste0("Var_",i),cbind(data.frame(Maxent_stack[i]),data.frame(Maxent_stack[i+15]),
                                data.frame(Maxent_stack[i+15*2]),data.frame(Maxent_stack[i+15*3]),
                                data.frame(Maxent_stack[i+15*4]),data.frame(Maxent_stack[i+15*5]),
                                data.frame(Maxent_stack[i+15*6]),data.frame(Maxent_stack[i+15*7]),
                                data.frame(Maxent_stack[i+15*8]),data.frame(Maxent_stack[i+15*9]),
                                data.frame(Maxent_stack[i+15*10])[selection,])[,-VariableSelection])
}

df_Final=rbind(Var_1,Var_2, Var_3,Var_4, Var_5,Var_6,Var_7,Var_8, Var_9, Var_10,
               Var_11,Var_12,Var_13,Var_14,Var_15)

# converting the matrix to data  


R_colNames = c("Variable0","Avg_X","Var_Y1","Var_Y2","Var_Y3","Var_Y4","Var_Y5","Var_Y6","Var_Y7","Var_Y8","Var_Y9","Var_Y10","Avg_Y")
colnames(df_Final)=R_colNames

unique(df_Final$Variable0) # This is to see the order of variables in the dataframe.

Var_list_sorted=c("Mean temp. of coldest quarter","Precipitation of driest quarter", 
                  "Precipitation of warmest quarter","Precipitation of coldest quarter","Density of livestock",
                  "Elevation","Density of water","Density of trees","Density of flooded vegetation",
                  "Density of crops","Density of built area","Density of bare ground",
                  "Density of rangeland", "Distance to road","Slope")

Var_list_num=rep(1:Var_Num,1)
Variables11022_1 = rep(Var_list_sorted, each = 501)
Variables11022_2 = rep(Var_list_num, each = 501)
df_Final$Variable = Variables11022_1
df_Final$Variable1 = Variables11022_2

df_Final$Variable2=paste(df_Final$Variable1,"_",df_Final$Variable)
Head1=head(df_Final)

Maxent_Var_Resp_All3=df_Final

dimension=dim(Maxent_Var_Resp_All3)
Min_Y=c()

for (i in 1:dimension[1]) {
  Min_Y1=min(Maxent_Var_Resp_All3[i,c(3:12)])
  Min_Y=c(Min_Y,Min_Y1)
}

Max_Y=c()

for (i in 1:dimension[1]) {
  Max_Y1=max(Maxent_Var_Resp_All3[i,c(3:12)])
  Max_Y=c(Max_Y,Max_Y1)
}

Maxent_Var_Resp_All3$Max_Y=Max_Y
Maxent_Var_Resp_All3$Min_Y=Min_Y

Head2=rbind(head(Maxent_Var_Resp_All3),tail(Maxent_Var_Resp_All3))

Maxent_Var_Resp_All2=Maxent_Var_Resp_All3

Maxent_Var_Resp_All1=subset(Maxent_Var_Resp_All2, Variable1 %in% c(11,3,7,14,6,12))

plot1 = ggplot(data = Maxent_Var_Resp_All1, aes(x=Avg_X, y=Avg_Y)) +
  geom_ribbon(aes(ymin = Min_Y, ymax = Max_Y), alpha = 1.0,
              fill = "#0000FF", color = "transparent")+
  geom_line(linetype = "solid", col = "#DF0020", size = 0.6)+
  facet_wrap( ~factor(Variable, levels=c("Density of built area",
                                          "Precipitation of warmest quarter",
                                          "Density of water","Distance to road",
                                          "Elevation","Density of bare ground")),
              scales="free",ncol = 3)+
  ##facet_grid( ~ Variable)+
  theme_bw() + scale_colour_tableau()+
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5))+
  labs(y ="    ", x = NULL)+ylim(0,1)+
  labs(y = "Logistic output", x = "Variable values")