#Previously to run the code please UPLOAD the workspace "Dataset1.RData" with the needed data for obtaining the results displayed in the manuscript. Visit the readme archive to understand each variable 
# --------------------step0. Preload the needed R packages. The following code will install and load the required packages if they are not preinstalled in your computer--------------#
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = FALSE)
sapply(pkg, require, character.only = TRUE)
}
packages<-c("tidyverse","ggplot2","ggrepel","gridExtra","V.PhyloMaker","MCMCglmm","bayestestR")

ipak(packages)
#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#Hypothesis 1. Facilitation enables the establishment of rare, sensitive species
#----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#----------------------step1. Generate phylogenetic relationship among the species presented in gyp_recruitment1 in order to phylogenetically inform the MCMCglmm--------#
tree<-phylo.maker(sp.list=phylo,tree=GBOTB.extended,nodes=nodes.info.1,scenario="S2",r=1)
mytree<-as.phylo(tree$scenario.2[[1]])
tree<-keep.tip(mytree,gyp_recruitment1$animal)#We only retain those species found in gyp_recruitment1
tree$node.label<-c("a","b","c","d","e","f","g","h","i","j","k","l","m","n","ñ","p","q","r","s")#We avoid node label with the same name that would be  problematic further

# --------------------- Step2. Check hypothesis. MCMCglmm and results-----------------------------------------------------------------------------------------------------------------------------------------#
#1 Model assessing the relationship between facilitation (positives RII) and stress sensitivity ( g)
RII_G<-(MCMCglmm(RII~g, random =~animal +Localidad , pedigree = tree ,family="gaussian",data=gyp_recruitment1,nitt=600000,burnin=100000,thin=100,pr=TRUE))
summary(RII_G)


# 2 Model assessing the relationship between Relative abundance (Ra) and stress sensitivity (g)

RII_Ra<-(MCMCglmm(Ra~g, random =~animal +Localidad , pedigree = tree ,family="gaussian",data=gyp_recruitment1,nitt=600000,burnin=100000,thin=100,pr=TRUE))
summary(RII_Ra)




#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#Hypothesis 2. Facilitation balances species performances
#----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#--------------------------------Step1. Calculation of the NPD with real performance data-------------------------------------------#
PD<-list()#list saving the localities
for (i in 1:length(levels(performance$site_ID)))
{

mysite<-droplevels(performance[performance$site_ID==levels(performance$site_ID)[i],])
mysp<-levels(mysite[,"sp"])
mysite_ID<-paste(unique(mysite$site_ID))

myinter<-combn(mysp,2)#possible pairs per locality and type of soils
par<-list()#new list where every possible interaction is saved

for( j in 1:dim(myinter)[2])
{
par[[j]]<-c(myinter[1,j],myinter[2,j])
}#close the j


PD_site<-list()
for(k in 1:length(par))
{
mypar<-droplevels(mysite[mysite$sp%in%par[[k]],])
mypar_solas<-mypar[mypar$Posicion=="sola",]
mypar_asociadas<-mypar[mypar$Posicion=="asociada",]
#Non-associated#
sp1s<-levels(mypar_solas$sp)[1]
sp2s<-levels(mypar_solas$sp)[2]
n_sp1s<-sum(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[1],"abundancia"])
n_sp2s<-sum(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[2],"abundancia"])

mean_sp1s<-median(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[1],"cover"])
mean_sp2s<-median(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[2],"cover"])
#Associated
sp1a<-levels(mypar_asociadas$sp)[1]
sp2a<-levels(mypar_asociadas$sp)[2]
n_sp1a<-sum(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[1],"abundancia"])
n_sp2a<-sum(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[2],"abundancia"])


mean_sp1a<-median(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[1],"cover"])
mean_sp2a<-median(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[2],"cover"])

PD_site[[k]]<-cbind(mysite_ID,sp1s,sp2s,mean_sp1s,mean_sp2s,n_sp1s,n_sp2s,sp1a,sp2a,mean_sp1a,mean_sp2a,n_sp1a,n_sp2a)


}#close the k 

PD[[i]]<-data.frame(do.call("rbind",PD_site))

}#close the i

PD_db<-data.frame(do.call("rbind",PD))
PD_db<-(PD_db[,setdiff(colnames(PD_db),c("sp1a","sp2a"))])
pp<-PD_db[,setdiff(colnames(PD_db),c("mysite_ID","sp1s","sp2s"))]%>%
mutate_if(is.character,as.numeric)
PD_db<-as.data.frame(cbind(PD_db[,c("mysite_ID","sp1s","sp2s")],pp))
PD_db$site_ID<-as.factor(PD_db$mysite_ID)
PD_db$sp1s<-as.factor(PD_db$sp1s)
PD_db$sp2s<-as.factor(PD_db$sp2s)
#for calculated PD (PERFORMANCE DIFFERENCES)#
PD_db$s1<-pmax(PD_db$mean_sp1s,PD_db$mean_sp2s)
PD_db$s2<-pmin(PD_db$mean_sp1s,PD_db$mean_sp2s)
PD_db$a1<-pmax(PD_db$mean_sp1a,PD_db$mean_sp2a)
PD_db$a2<-pmin(PD_db$mean_sp1a,PD_db$mean_sp2a)
#################################################
PD_db$diff_solas<-(PD_db$s1/PD_db$s2)
PD_db$diff_asociadas<-(PD_db$a1/PD_db$a2)


PD_db$NPD<-(PD_db$diff_asociadas-PD_db$diff_solas)# Effec of Neighborhood in performance differences
PD_db<-separate(PD_db,"site_ID",c("Localidad","Suelo"),sep="-")
mean(PD_db$NPD)# Overall Effec of Neighborhood in performance differences

#------------Step3. Null model of NPDs produced by reshuffling the performance dataset-----------------

nul<-list()
for (h in 1:10000)#number of times that we reshuffled the dataset

{


mysp_loc_all<-list()
for (z in 1:length(levels(performance$ID)))
{
mysp_loc<-droplevels(performance[performance$ID==levels(performance$ID)[z],])
mysp_loc$Posicion<-sample(mysp_loc$Posicion,replace=F)
mysp_loc_all[[z]]<-mysp_loc
}

performance_rand<-do.call(rbind,mysp_loc_all)#

PD<-list()
for (i in 1:length(levels(performance_rand$site_ID)))
{

mysite<-droplevels(performance_rand[performance_rand$site_ID==levels(performance_rand$site_ID)[i],])
mysp<-levels(mysite[,"sp"])
mysite_ID<-paste(unique(mysite$site_ID))

myinter<-combn(mysp,2)
par<-list()

for( j in 1:dim(myinter)[2])
{
par[[j]]<-c(myinter[1,j],myinter[2,j])
}#cierra la j

PD_site<-list()
for(k in 1:length(par))
{
mypar<-droplevels(mysite[mysite$sp%in%par[[k]],])
mypar_solas<-mypar[mypar$Posicion=="sola",]
mypar_asociadas<-mypar[mypar$Posicion=="asociada",]
#Non-associated#
sp1s<-levels(mypar_solas$sp)[1]
sp2s<-levels(mypar_solas$sp)[2]
n_sp1s<-sum(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[1],"abundancia"])
n_sp2s<-sum(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[2],"abundancia"])


mean_sp1s<-median(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[1],"cover"])
mean_sp2s<-median(mypar_solas[mypar_solas$sp==levels(mypar_solas$sp)[2],"cover"])


#Associated·
sp1a<-levels(mypar_asociadas$sp)[1]
sp2a<-levels(mypar_asociadas$sp)[2]
n_sp1a<-sum(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[1],"abundancia"])
n_sp2a<-sum(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[2],"abundancia"])

mean_sp1a<-median(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[1],"cover"])
mean_sp2a<-median(mypar_asociadas[mypar_asociadas$sp==levels(mypar_asociadas$sp)[2],"cover"])

PD_site[[k]]<-cbind(mysite_ID,sp1s,sp2s,mean_sp1s,mean_sp2s,n_sp1s,n_sp2s,sp1a,sp2a,mean_sp1a,mean_sp2a,n_sp1a,n_sp2a)


}#cierra la k

PD[[i]]<-data.frame(do.call("rbind",PD_site))

}#cierra la i

PD_db_null<-data.frame(do.call("rbind",PD))
PD_db_null<-(PD_db_null[,setdiff(colnames(PD_db_null),c("sp1a","sp2a"))])
pp<-PD_db_null[,setdiff(colnames(PD_db_null),c("mysite_ID","sp1s","sp2s"))]%>%
mutate_if(is.character,as.numeric)
PD_db_null<-as.data.frame(cbind(PD_db_null[,c("mysite_ID","sp1s","sp2s")],pp))
PD_db_null$site_ID<-as.factor(PD_db_null$mysite_ID)
PD_db_null$sp1s<-as.factor(PD_db_null$sp1s)
PD_db_null$sp2s<-as.factor(PD_db_null$sp2s)
PD_db_null$s1<-pmax(PD_db_null$mean_sp1s,PD_db_null$mean_sp2s)
PD_db_null$s2<-pmin(PD_db_null$mean_sp1s,PD_db_null$mean_sp2s)
PD_db_null$a1<-pmax(PD_db_null$mean_sp1a,PD_db_null$mean_sp2a)
PD_db_null$a2<-pmin(PD_db_null$mean_sp1a,PD_db_null$mean_sp2a)
#################################################
PD_db_null$diff_solas<-(PD_db_null$s1/PD_db_null$s2)
PD_db_null$diff_asociadas<-(PD_db_null$a1/PD_db_null$a2)

PD_db_null$NPD<-(PD_db_null$diff_asociadas-PD_db_null$diff_solas)
PD_db_null<-separate(PD_db_null,"site_ID",c("Localidad","Suelo"),sep="-")

nul[[h]]<-PD_db_null$NPD
}
null_model<-NULL
for (t in 1 :10000)
{
null_model[t]<-mean(nul[[t]])
}

#----------------Step.4 Comparing the observed NPD with the NPDa derived from the null model------------------#

ci(null_model,ci=0.95,method="HDI")# null model average and 95% credible interval
mean(PD_db$NPD)# oberved association effect

#-----------------------Figures---------------------------------#
# We only show the code for figure 1 which is the only one displaying results
# Fig.1 R code for generating Fig.1

all<-ggplot(gyp_recruitment1, aes(x=g, y=RII)) + 
  geom_point()+
  geom_smooth(method=lm, se=FALSE)

Crevillente<-ggplot(gyp_recruitment1[gyp_recruitment1$Localidad=="Crevillente",], aes(x=g, y=RII)) + 
  geom_point()+geom_smooth(method=lm, se=FALSE)+geom_text(aes(label=animal))

Perrio<-ggplot(gyp_recruitment1[gyp_recruitment1$Localidad=="Perrio",], aes(x=g, y=RII)) + 
  geom_point()+geom_smooth(method=lm, se=FALSE)+geom_text(aes(label=animal))

La_mola<-ggplot(gyp_recruitment1[gyp_recruitment1$Localidad=="La_mola",], aes(x=g, y=RII)) + 
  geom_point()+geom_smooth(method=lm, se=FALSE)+geom_text(aes(label=animal))

Salinas<-ggplot(gyp_recruitment1[gyp_recruitment1$Localidad=="Salinas",], aes(x=g, y=RII)) + 
  geom_point()+geom_smooth(method=lm, se=FALSE)+geom_text(aes(label=animal))
	
grid.arrange(all, Crevillente, Perrio, La_mola, Salinas,ncol = 3, nrow = 2)