#Manuscript "No general support of functional diversity enhancing resilience across terrestrial plant communities" ### Executing the following code performs the main analysis and figures presented in the manuscript. ##The code is structured as follows: #1. Required packages #2. Principal Coordinates Analysis (PCoA) #2.1. Description and objectives of the analysis #2.2. Analysis for the complete data set of Resilience and Species richness #2.2.1. Analysis per ecosystem type (Including both biomass and floristic composition data) #2.2.1.1. Woody ecosystems #2.2.1.2. Herbaceous ecosystems #2.2.2. Analysis per ecosystem type and response category #2.2.2.1. Resilience of biomass in woody ecosystems #2.2.2.2. Resilience of floristic composition in woody ecosystems #2.2.2.3. Resilience of biomass in herbaceous ecosystems #2.3. Analysis for the data set of Resilience and Indices of Functional Diversity #2.3.1. Analysis per ecosystem type (Including both biomass and floristic composition data) #2.3.1.1. Woody ecosystems #2.3.1.2. Herbaceous ecosystems #2.3.2. Analysis per ecosystem type and response category #2.3.2.1. Resilience of biomass in woody ecosystems #2.3.2.2. Resilience of floristic composition in woody ecosystems #2.3.2.3. Resilience of biomass in herbaceous ecosystems #3. Meta-analysis of the correlation between Resilience and Richness and Functional Diversity #3.1. Description and objectives of the analysis #3.2. Analysis for the complete data set of Resilience and Species richness #3.2.1. Analysis per ecosystem type and response category #3.2.1.1. Resilience of biomass in herbaceous ecosystems and Species richness #3.2.1.2. Resilience of biomass in woody ecosystems and Species richness #3.2.1.3. Resilience of floristic composition in woody ecosystems and Species richness #3.3. Analysis for the complete data set of Resilience and Functional Richness #3.3.1. Analysis per ecosystem type and response category #3.3.1.1. Resilience of biomass in herbaceous ecosystems and Functional Richness #3.3.1.2. Resilience of biomass in woody ecosystems and Functional Richness #3.3.1.3. Resilience of floristic composition in woody ecosystems and Functional Richness #3.4. Analysis for the complete data set of Resilience and Functional Divergence #3.4.1. Analysis per ecosystem type and response category #3.4.1.1. Resilience of biomass in herbaceous ecosystems and Functional Divergence #3.4.1.2. Resilience of biomass in woody ecosystems and Functional Divergence #3.4.1.3. Resilience of floristic composition in woody ecosystems and Functional Divergence #3.5. Analysis for the complete data set of Resilience and Functional Dispersion #3.5.1. Analysis per ecosystem type and response category #3.5.1.1. Resilience of biomass in herbaceous ecosystems and Functional Dispersion #3.5.1.2. Resilience of biomass in woody ecosystems and Functional Dispersion #3.5.1.3. Resilience of floristic composition in woody ecosystems and Functional Dispersion #3.6. Analysis for the complete data set of Resilience and Functional Evenness #3.6.1. Analysis per ecosystem type and response category #3.6.1.1. Resilience of biomass in herbaceous ecosystems and Functional Evenness #3.6.1.2. Resilience of biomass in woody ecosystems and Functional Evenness #3.6.1.3. Resilience of floristic composition in woody ecosystems and Functional Evenness #3.7. Analysis for the complete data set of Resilience and Small plant Dominance #3.7.1. Analysis per ecosystem type and response category #3.7.1.1. Resilience of biomass in herbaceous ecosystems and Small plant Dominance #3.7.1.2. Resilience of biomass in woody ecosystems and Small plant Dominance #3.7.1.3. Resilience of floristic composition in woody ecosystems and Small plant Dominance #3.8. Analysis for the complete data set of Resilience and Acquisitive plant Dominance #3.8.1. Analysis per ecosystem type and response category #3.8.1.1. Resilience of biomass in herbaceous ecosystems and Acquisitive plant Dominance #3.8.1.2. Resilience of biomass in woody ecosystems and Acquisitive plant Dominance #3.8.1.3. Resilience of floristic composition in woody ecosystems and Acquisitive plant Dominance #4. Figures #4.1. Figure 1 - MAP #4.2. Figure 2 - Grand mean effect sizes of the different meta-analyses #4.3. Figure 3 - Grand mean effect sizes of the different meta-analyses grouping by Ecosystem type and Response variable #4.3.1. Figure 3 - Resilience of biomass in herbaceous ecosystems #4.3.2. Figure 3 - Resilience of biomass in woody ecosystems #4.3.3. Figure 3 - Resilience of floristic composition in woody ecosystems #4.4. Figure 4 - Relationship between the axes of moderator variation and Zscores Herbaceous ecosystems #4.4.1. Figure 4 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and FDis) #4.4.2. Figure 4 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and Acquisitive plant dominance) #4.5. Figure 5 - Relationship between the axes of moderator variation and Zscores Woody ecosystems #4.5.1. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and FDiv) #4.5.2. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and Acquisitive plant dominance) #4.5.3. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and FRic) #4.5.4. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and FDiv) #4.5.5. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and Acquisitive plant dominance) #4.6. Appendix 4 - PCoA Data set on resilience and Species Richness #4.6.1. Appendix 4 - Complete dataset #4.6.2. Appendix 4 - Woody ecosystems #4.6.3. Appendix 4 - Herbaceous ecosystems #4.6.4. Appendix 4 - Resilience of biomass in woody ecosystems #4.6.5. Appendix 4 - Resilience of floristic compoitionin woody ecosystems #4.6.6. Appendix 4 - Resilience of biomass herbaceous ecosystems #4.7. Appendix 5 - PCoA Data set on resilience and Functional Diverstity indices #4.7.1. Appendix 5 - Complete dataset #4.7.2. Appendix 5 - Woody ecosystems #4.7.3. Appendix 5 - Herbaceous ecosystems #4.7.4. Appendix 5 - Resilience of biomass in woody ecosystems #4.7.5. Appendix 5 - Resilience of floristic compoitionin woody ecosystems #4.7.6. Appendix 5 - Resilience of biomass herbaceous ecosystems #1. Required packages #library(tidyr) library(tidyverse) library(dplyr) library(metafor) #library(maps) library(ggplot2) #library(maptools) library(FD) library(glue) library(patchwork) library(egg) ############################################### #2. Principal Coordinate Analysis (PCoA) #2.1. We used Principal Coordinate Analyses (PCoAs), based on Gower’s dissimilarities, across all moderators # to divide our data sets into studies with similar ecological and methodological study-characteristic and # delineate common axes of moderators’ variation (see Appendix 4 and 5). By subdividing the data sets, the # number of moderators and, as a consequence their correlation, was reduced. #2.2. Analysis for the complete data set of Resilience and Species richness (Results from this analysis are presented in Appendix 4). #Get data Resilience.SpeciesRichness.Correlation<-read.table("ResilienceSR.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) ##Select variables COMPLETE<-Resilience.SpeciesRichness.Correlation %>% select(Ecosystem.type,Response.category,Biome, Response.component,Disturbance.type, LogDisturbance.days,LogTemporal.scale,Community.origin,Study.design,Study.type, Resilience.category) COMPLETE %>%names #Reduce names for simplicity COMPLETE<-COMPLETE %>% rename(ET=Ecosystem.type, RESP=Response.category, COMP=Response.component, DIST.TYPE=Disturbance.type, DIST.DAYS=LogDisturbance.days, TEMP.SCALE=LogTemporal.scale, ORIG=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST.TYPE=case_match(DIST.TYPE, "Wet conditions"~"Wet", "Dry conditions"~ "Dry", "Direct biomass removal"~"Biomass removal")) #characters need to be transformed to factors in order to use Gower distance function COMPLETE<-COMPLETE%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixCOMPLETE <- gowdis(COMPLETE) #PCoA analysis test.pcoaCOMPLETE = capscale(test.dist.matrixCOMPLETE ~ 1) SCORESCOMPLETE<-scores(test.pcoaCOMPLETE) #Scores for each site PCOA1complete<-SCORESCOMPLETE$sites[,1] #Scores in PCoA1 PCOA2complete<-SCORESCOMPLETE$sites[,2] #Scores in PCoA2 #build a table to plot later PCOAcomplete<-tibble(PCOA1complete,PCOA2complete,Resilience.SpeciesRichness.Correlation$ID_Si)%>% rename(PCoA1=PCOA1complete, PCoA2=PCOA2complete, ID_Si="Resilience.SpeciesRichness.Correlation$ID_Si") #Eigenvalues for each axis eigCOMPLETE<-eigenvals(test.pcoaCOMPLETE) eigCOMPLETEpor<-(eigCOMPLETE/sum(eigCOMPLETE))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANcomplete <- envfit(test.pcoaCOMPLETE, COMPLETE, na.rm = TRUE) efitVEGANcomplete.vector<-data.frame(efitVEGANcomplete$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANcomplete.factor<-data.frame(efitVEGANcomplete$factors$centroids) #returns the centroids of the factor #2.2.1. Analysis per ecosystem type (Including both biomass and floristic composition data) #2.2.1.1. Woody ecosystems #Get data Resilience.SpeciesRichness.Woody.Correlation<-read.table("ResilienceSR-woody.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) WOODY<-Resilience.SpeciesRichness.Woody.Correlation%>% ungroup()%>% select(Response.category,Biome,Response.component,LogDisturbance.days,Study.design,Study.type, Resilience.category,LogTemporal.scale) WOODY%>%names #Reduce names for simplicity WOODY<-WOODY %>% rename(RESP=Response.category, COMP=Response.component, DIST.DAYS=LogDisturbance.days, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function WOODY<-WOODY%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixWOODY <- gowdis(WOODY) #PCoA analysis test.pcoaWOODY = capscale(test.dist.matrixWOODY ~ 1) SCORESWOODY<-scores(test.pcoaWOODY) #Scores for each site PCOA1WOODY<-SCORESWOODY$sites[,1] #Scores for PCoA1 PCOA2WOODY<-SCORESWOODY$sites[,2] #Scores for PCoA2 #build a table to plot later PCOAWOODY<-tibble(PCOA1WOODY,PCOA2WOODY,Resilience.SpeciesRichness.Woody.Correlation$ID_Si)%>% rename(PCoA1=PCOA1WOODY, PCoA2=PCOA2WOODY, ID_Si="Resilience.SpeciesRichness.Woody.Correlation$ID_Si") #Eigenvalues for each axis eigWOODY<-eigenvals(test.pcoaWOODY) eigWOODYpor<-(eigWOODY/sum(eigWOODY))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANWOODY <- envfit(test.pcoaWOODY, WOODY, na.rm = TRUE) efitVEGANWOODY.vector<-data.frame(efitVEGANWOODY$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANWOODY.factor<-data.frame(efitVEGANWOODY$factors$centroids)#returns the centroids of the factor #2.2.1.2. Herbaceous ecosystems #Get data Resilience.SpeciesRichness.Herbaceous.Correlation<-read.table("ResilienceSR-herbaceous.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) HERBACEOUS<-Resilience.SpeciesRichness.Herbaceous.Correlation%>% ungroup()%>% select(Biome,Response.category,Disturbance.type,LogDisturbance.days,Community.origin, Study.design,Study.type,Resilience.category,LogTemporal.scale) HERBACEOUS%>%names #Reduce names for simplicity HERBACEOUS<-HERBACEOUS %>% rename(RESP=Response.category, DIST=Disturbance.type, DIST.DAYS=LogDisturbance.days, ORI=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST=case_match(DIST, "Direct biomass removal"~"Biomass removal", "Dry conditions"~"Dry", "Wet conditions"~"Wet")) #characters need to be transformed to factors in order to use Gower distance function HERBACEOUS<-HERBACEOUS%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixHERBACEOUS <- gowdis(HERBACEOUS) #PCoA analysis test.pcoaHERBACEOUS = capscale(test.dist.matrixHERBACEOUS ~ 1) SCORESHERBACEOUS<-scores(test.pcoaHERBACEOUS) #Scores for each site PCOA1HERBACEOUS<-SCORESHERBACEOUS$sites[,1] #Scores for PCoA1 PCOA2HERBACEOUS<-SCORESHERBACEOUS$sites[,2] #Scores for PCoA2 #build a table to plot later PCOAHERBACEOUS<-tibble(PCOA1HERBACEOUS,PCOA2HERBACEOUS,Resilience.SpeciesRichness.Herbaceous.Correlation$ID_Si)%>% rename(PCoA1=PCOA1HERBACEOUS, PCoA2=PCOA2HERBACEOUS, ID_Si="Resilience.SpeciesRichness.Herbaceous.Correlation$ID_Si") #Eigenvalues for each axis eigHERBACEOUS<-eigenvals(test.pcoaHERBACEOUS) eigHERBACEOUSpor<-(eigHERBACEOUS/sum(eigHERBACEOUS))*100 #envfit() fits environmental vectors or factors onto and ordination. efitVEGANHERBACEOUS <- envfit(test.pcoaHERBACEOUS, HERBACEOUS, na.rm = TRUE) efitVEGANHERBACEOUS.vector<-data.frame(efitVEGANHERBACEOUS$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANHERBACEOUS.factor<-data.frame(efitVEGANHERBACEOUS$factors$centroids)#returns the centroids of the factor #2.2.2. Analysis per ecosystem type and Response category #2.2.2.1. Resilience of biomass in woody ecosystems #Get data Resilience.SpeciesRichness.Biomass_W.Correlation<-read.table("ResilienceSR-BiomassW.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) Biomass_W<-Resilience.SpeciesRichness.Biomass_W.Correlation%>% ungroup()%>% select(Biome,Response.component,LogDisturbance.days, Study.design,Study.type,LogTemporal.scale) Biomass_W%>%names #Reduce names for simplicity Biomass_W<-Biomass_W %>% rename(COMP=Response.component, DIST.DAYS=LogDisturbance.days, SD=Study.design, ST=Study.type, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function Biomass_W<-Biomass_W%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixBiomass_W <- gowdis(Biomass_W) #PCoA analysis test.pcoaBiomass_W = capscale(test.dist.matrixBiomass_W ~ 1) SCORESBiomass_W<-scores(test.pcoaBiomass_W) #Scores for each site PCOA1Biomass_W<-SCORESBiomass_W$sites[,1] #Scores for PCoA1 PCOA2Biomass_W<-SCORESBiomass_W$sites[,2] #Scores for PCoA2 #build a table to plot later PCOABiomass_W<-tibble(PCOA1Biomass_W,PCOA2Biomass_W,Resilience.SpeciesRichness.Biomass_W.Correlation$ID_Si)%>% rename(PCoA1=PCOA1Biomass_W, PCoA2=PCOA2Biomass_W, ID_Si="Resilience.SpeciesRichness.Biomass_W.Correlation$ID_Si") #Eigenvalues for each axis eigBiomass_W<-eigenvals(test.pcoaBiomass_W) eigBiomass_Wpor<-(eigBiomass_W/sum(eigBiomass_W))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANBiomass_W <- envfit(test.pcoaBiomass_W, Biomass_W, na.rm = T) efitVEGANBiomass_W.vector<-data.frame(efitVEGANBiomass_W$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANBiomass_W.factor<-data.frame(efitVEGANBiomass_W$factors$centroids)#returns the centroids of the factor #2.2.2.2. Resilience of floristic composition in woody ecosystems #Get data Resilience.SpeciesRichness.Composition_W.Correlation<-read.table("ResilienceSR-CompositionW.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) COMPOSITION_W<-Resilience.SpeciesRichness.Composition_W.Correlation%>% ungroup()%>% select(Biome,Response.component,LogDisturbance.days,Study.design,Study.type,Resilience.category, LogTemporal.scale) COMPOSITION_W%>%names #Reduce names for simplicity COMPOSITION_W<-COMPOSITION_W %>% rename(COMP=Response.component, DIST.DAYS=LogDisturbance.days, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function COMPOSITION_W<-COMPOSITION_W%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixCOMPOSITION_W <- gowdis(COMPOSITION_W) #PCoA analysis test.pcoaCOMPOSITION_W = capscale(test.dist.matrixCOMPOSITION_W ~ 1) SCORESCOMPOSITION_W<-scores(test.pcoaCOMPOSITION_W) #Scores for each site PCOA1COMPOSITION_W<-SCORESCOMPOSITION_W$sites[,1] #Scores for PCoA1 PCOA2COMPOSITION_W<-SCORESCOMPOSITION_W$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOACOMPOSITION_W<-tibble(PCOA1COMPOSITION_W,PCOA2COMPOSITION_W,Resilience.SpeciesRichness.Composition_W.Correlation$ID_Si)%>% rename(PCoA1=PCOA1COMPOSITION_W, PCoA2=PCOA2COMPOSITION_W, ID_Si="Resilience.SpeciesRichness.Composition_W.Correlation$ID_Si") #Eigenvalues for each axis eigCOMPOSITION_W<-eigenvals(test.pcoaCOMPOSITION_W) eigCOMPOSITION_Wpor<-(eigCOMPOSITION_W/sum(eigCOMPOSITION_W))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANCOMPOSITION_W <- envfit(test.pcoaCOMPOSITION_W, COMPOSITION_W, na.rm = TRUE) efitVEGANCOMPOSITION_W.vector<-data.frame(efitVEGANCOMPOSITION_W$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANCOMPOSITION_W.factor<-data.frame(efitVEGANCOMPOSITION_W$factors$centroids)#returns the centroids of the factor #2.2.2.3. Resilience of biomass in herbaceous ecosystems #Get data Resilience.SpeciesRichness.Biomass_H.Correlation<-read.table("ResilienceSR-BiomassH.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) Biomass_H<-Resilience.SpeciesRichness.Biomass_H.Correlation%>% ungroup()%>% select(Biome,Disturbance.type,LogDisturbance.days,Community.origin,Study.design, Study.type,Resilience.category,LogTemporal.scale) Biomass_H%>%names #Reduce names for simplicity Biomass_H<-Biomass_H %>% rename(DIST=Disturbance.type, DIST.DAYS=LogDisturbance.days, ORI=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST=case_match(DIST, "Direct biomass removal"~"Biomass removal", "Dry conditions"~"Dry", "Wet conditions"~"Wet")) #characters need to be transformed to factors in order to use Gower distance function Biomass_H<-Biomass_H%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixBiomass_H <- gowdis(Biomass_H) #PCoA analysis test.pcoaBiomass_H = capscale(test.dist.matrixBiomass_H ~ 1) SCORESBiomass_H<-scores(test.pcoaBiomass_H) #Scores for each site PCOA1Biomass_H<-SCORESBiomass_H$sites[,1] #Scores for PCoA1 PCOA2Biomass_H<-SCORESBiomass_H$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOABiomass_H<-tibble(PCOA1Biomass_H,PCOA2Biomass_H,Resilience.SpeciesRichness.Biomass_H.Correlation$ID_Si)%>% rename(PCoA1=PCOA1Biomass_H, PCoA2=PCOA2Biomass_H, ID_Si="Resilience.SpeciesRichness.Biomass_H.Correlation$ID_Si") #Eigenvalues for each axis eigBiomass_H<-eigenvals(test.pcoaBiomass_H) eigBiomass_Hpor<-(eigBiomass_H/sum(eigBiomass_H))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANBiomass_H <- envfit(test.pcoaBiomass_H, Biomass_H, na.rm = TRUE) efitVEGANBiomass_H.vector<-data.frame(efitVEGANBiomass_H$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANBiomass_H.factor<-data.frame(efitVEGANBiomass_H$factors$centroids)#returns the centroids of the factor ####End of this analysis#### ############################ #2.3. Analysis for the complete data set of Resilience and Indices of Functional Diversity (Results from this analysis are presented in Appendix 5) #Get data Resilience.FunctionalDiversity.Correlation<-read.table("ResilienceFDindices.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) COMPLETE.FD<-Resilience.FunctionalDiversity.Correlation %>% select(Ecosystem.type,Response.category,Biome, Response.component,Disturbance.type, LogDisturbance.days,LogTemporal.scale,Community.origin,Study.design,Study.type, Resilience.category) COMPLETE.FD %>%names #Reduce names for simplicity COMPLETE.FD<-COMPLETE.FD %>% rename(ET=Ecosystem.type, RESP=Response.category, COMP=Response.component, DIST.TYPE=Disturbance.type, DIST.DAYS=LogDisturbance.days, TEMP.SCALE=LogTemporal.scale, ORIG=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST.TYPE=case_match(DIST.TYPE, "Wet conditions"~"Wet", "Dry conditions"~ "Dry", "Direct biomass removal"~"Biomass removal")) #characters need to be transformed to factors in order to use Gower distance function COMPLETE.FD<-COMPLETE.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixCOMPLETE.FD <- gowdis(COMPLETE.FD) #PCoA analysis test.pcoaCOMPLETE.FD = capscale(test.dist.matrixCOMPLETE.FD ~ 1) SCORESCOMPLETE.FD<-scores(test.pcoaCOMPLETE.FD) #Scores for each site PCOA1complete.FD<-SCORESCOMPLETE.FD$sites[,1] #Scores in PCoA1 PCOA2complete.FD<-SCORESCOMPLETE.FD$sites[,2] #Scores in PCoA2 #build a table to plot later PCOAcomplete.FD<-tibble(PCOA1complete.FD,PCOA2complete.FD,Resilience.FunctionalDiversity.Correlation$ID_Si)%>% rename(PCoA1=PCOA1complete.FD, PCoA2=PCOA2complete.FD, ID_Si="Resilience.FunctionalDiversity.Correlation$ID_Si") #Eigenvalues for each axis eigCOMPLETE.FD<-eigenvals(test.pcoaCOMPLETE.FD) eigCOMPLETE.FDpor<-(eigCOMPLETE.FD/sum(eigCOMPLETE.FD))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANcomplete.FD <- envfit(test.pcoaCOMPLETE.FD, COMPLETE.FD, na.rm = TRUE) efitVEGANcomplete.FD.vector<-data.frame(efitVEGANcomplete.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANcomplete.FD.factor<-data.frame(efitVEGANcomplete.FD$factors$centroids) #returns the centroids of the factor #2.3.1. Analysis per ecosystem type (Including both biomass and floristic composition data) #2.3.1.1. Woody ecosystems #Get data Resilience.FunctionalDiversity.Woody.Correlation<-read.table("ResilienceFDindices-woody.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) WOODY.FD<-Resilience.FunctionalDiversity.Woody.Correlation%>% ungroup()%>% select(Response.category,Biome,Response.component,LogDisturbance.days,Study.design,Study.type, Resilience.category,LogTemporal.scale) WOODY.FD%>%names #Reduce names for simplicity WOODY.FD<-WOODY.FD %>% rename(RESP=Response.category, COMP=Response.component, DIST.DAYS=LogDisturbance.days, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function WOODY.FD<-WOODY.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixWOODY.FD <- gowdis(WOODY.FD) #PCoA analysis test.pcoaWOODY.FD = capscale(test.dist.matrixWOODY.FD ~ 1) SCORESWOODY.FD<-scores(test.pcoaWOODY.FD) #Scores for each site PCOA1WOODY.FD<-SCORESWOODY.FD$sites[,1] #Scores for PCoA1 PCOA2WOODY.FD<-SCORESWOODY.FD$sites[,2] #Scores for PCoA2 #build a table to plot later PCOAWOODY.FD<-tibble(PCOA1WOODY.FD,PCOA2WOODY.FD,Resilience.FunctionalDiversity.Woody.Correlation$ID_Si)%>% rename(PCoA1=PCOA1WOODY.FD, PCoA2=PCOA2WOODY.FD, ID_Si="Resilience.FunctionalDiversity.Woody.Correlation$ID_Si") #Eigenvalues for each axis eigWOODY.FD<-eigenvals(test.pcoaWOODY.FD) eigWOODY.FDpor<-(eigWOODY.FD/sum(eigWOODY.FD))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANWOODY.FD <- envfit(test.pcoaWOODY.FD, WOODY.FD, na.rm = TRUE) efitVEGANWOODY.FD.vector<-data.frame(efitVEGANWOODY.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANWOODY.FD.factor<-data.frame(efitVEGANWOODY.FD$factors$centroids)#returns the centroids of the factor #2.3.1.2. Herbaceous ecosystems #Get data Resilience.FunctionalDiversity.Herbaceous.Correlation<-read.table("ResilienceFDindices-herbaceous.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) HERBACEOUS.FD<-Resilience.FunctionalDiversity.Herbaceous.Correlation%>% ungroup()%>% select(Biome,Response.category,Disturbance.type,LogDisturbance.days,Community.origin, Study.design,Study.type,Resilience.category,LogTemporal.scale) HERBACEOUS.FD%>%names #Reduce names for simplicity HERBACEOUS.FD<-HERBACEOUS.FD %>% rename(RESP=Response.category, DIST=Disturbance.type, DIST.DAYS=LogDisturbance.days, ORI=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST=case_match(DIST, "Direct biomass removal"~"Biomass removal", "Dry conditions"~"Dry", "Wet conditions"~"Wet")) #characters need to be transformed to factors in order to use Gower distance function HERBACEOUS.FD<-HERBACEOUS.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixHERBACEOUS.FD <- gowdis(HERBACEOUS.FD) #PCoA analysis test.pcoaHERBACEOUS.FD = capscale(test.dist.matrixHERBACEOUS.FD ~ 1) SCORESHERBACEOUS.FD<-scores(test.pcoaHERBACEOUS.FD) #Scores for each site PCOA1HERBACEOUS.FD<-SCORESHERBACEOUS.FD$sites[,1] #Scores for PCoA1 PCOA2HERBACEOUS.FD<-SCORESHERBACEOUS.FD$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOAHERBACEOUS.FD<-tibble(PCOA1HERBACEOUS.FD,PCOA2HERBACEOUS.FD,Resilience.FunctionalDiversity.Herbaceous.Correlation$ID_Si)%>% rename(PCoA1=PCOA1HERBACEOUS.FD, PCoA2=PCOA2HERBACEOUS.FD, ID_Si="Resilience.FunctionalDiversity.Herbaceous.Correlation$ID_Si") #Eigenvalues for each axis eigHERBACEOUS.FD<-eigenvals(test.pcoaHERBACEOUS.FD) eigHERBACEOUS.FDpor<-(eigHERBACEOUS.FD/sum(eigHERBACEOUS.FD))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANHERBACEOUS.FD <- envfit(test.pcoaHERBACEOUS.FD, HERBACEOUS.FD, na.rm = TRUE) efitVEGANHERBACEOUS.FD.vector<-data.frame(efitVEGANHERBACEOUS.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANHERBACEOUS.FD.factor<-data.frame(efitVEGANHERBACEOUS.FD$factors$centroids)#returns the centroids of the factor #2.3.2. Analysis per ecosystem type and response category #2.3.2.1. Resilience of biomass in woody ecosystems #Get data Resilience.FunctionalDiversity.Biomass_W.Correlation<-read.table("ResilienceFDindices-BiomassW.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) Biomass_W.FD<-Resilience.FunctionalDiversity.Biomass_W.Correlation%>% ungroup()%>% select(Biome,Response.component,LogTemporal.scale, LogPlot.size.m2) Biomass_W.FD%>%names #Reduce names for simplicity Biomass_W.FD<-Biomass_W.FD %>% rename(COMP=Response.component, TEMP.SCALE=LogTemporal.scale, PLOT.Size=LogPlot.size.m2)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function Biomass_W.FD<-Biomass_W.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixBiomass_W.FD <- gowdis(Biomass_W.FD) #PCoA analysis test.pcoaBiomass_W.FD = capscale(test.dist.matrixBiomass_W.FD ~ 1) SCORESBiomass_W.FD<-scores(test.pcoaBiomass_W.FD) #Scores for each site PCOA1Biomass_W.FD<-SCORESBiomass_W.FD$sites[,1] #Scores for PCoA1 PCOA2Biomass_W.FD<-SCORESBiomass_W.FD$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOABiomass_W.FD<-tibble(PCOA1Biomass_W.FD,PCOA2Biomass_W.FD,Resilience.FunctionalDiversity.Biomass_W.Correlation$ID_Si)%>% rename(PCoA1=PCOA1Biomass_W.FD, PCoA2=PCOA2Biomass_W.FD, ID_Si="Resilience.FunctionalDiversity.Biomass_W.Correlation$ID_Si") #Eigenvalues for each axis eigBiomass_W.FD<-eigenvals(test.pcoaBiomass_W.FD) eigBiomass_W.FDpor<-(eigBiomass_W/sum(eigBiomass_W.FD))*100 #envfit() fits environmental vectors or factors into and ordination. efitVEGANBiomass_W.FD <- envfit(test.pcoaBiomass_W.FD, Biomass_W.FD, na.rm = T) efitVEGANBiomass_W.FD.vector<-data.frame(efitVEGANBiomass_W.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANBiomass_W.FD.factor<-data.frame(efitVEGANBiomass_W.FD$factors$centroids)#returns the centroids of the factor #2.3.2.2. Resilience of floristic composition in woody ecosystems #Get data Resilience.FunctionalDiversity.Comp_W.Correlation<-read.table("ResilienceFDindices-CompositionW.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) COMPOSITION_W.FD<-Resilience.FunctionalDiversity.Comp_W.Correlation%>% ungroup()%>% select(Biome,Response.component,LogTemporal.scale, LogPlot.size.m2) COMPOSITION_W.FD%>%names #Reduce names for simplicity COMPOSITION_W.FD<-COMPOSITION_W.FD %>% rename(COMP=Response.component, TEMP.SCALE=LogTemporal.scale, PLOT.SIZE=LogPlot.size.m2)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu")) #characters need to be transformed to factors in order to use Gower distance function COMPOSITION_W.FD<-COMPOSITION_W.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixCOMPOSITION_W.FD <- gowdis(COMPOSITION_W.FD) #PCoA analysis test.pcoaCOMPOSITION_W.FD = capscale(test.dist.matrixCOMPOSITION_W.FD ~ 1) SCORESCOMPOSITION_W.FD<-scores(test.pcoaCOMPOSITION_W.FD) #Scores for each site PCOA1COMPOSITION_W.FD<-SCORESCOMPOSITION_W.FD$sites[,1] #Scores for PCoA1 PCOA2COMPOSITION_W.FD<-SCORESCOMPOSITION_W.FD$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOACOMPOSITION_W.FD<-tibble(PCOA1COMPOSITION_W.FD,PCOA2COMPOSITION_W.FD,Resilience.FunctionalDiversity.Comp_W.Correlation$ID_Si)%>% rename(PCoA1=PCOA1COMPOSITION_W.FD, PCoA2=PCOA2COMPOSITION_W.FD, ID_Si="Resilience.FunctionalDiversity.Comp_W.Correlation$ID_Si") #Eigenvalues for each axis eigCOMPOSITION_W.FD<-eigenvals(test.pcoaCOMPOSITION_W.FD) eigCOMPOSITION_W.FDpor<-(eigCOMPOSITION_W.FD/sum(eigCOMPOSITION_W.FD))*100 #envfit() fits environmental vectors or factors onto and ordination. efitVEGANCOMPOSITION_W.FD <- envfit(test.pcoaCOMPOSITION_W.FD, COMPOSITION_W.FD, na.rm = TRUE) efitVEGANCOMPOSITION_W.FD.vector<-data.frame(efitVEGANCOMPOSITION_W.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANCOMPOSITION_W.FD.factor<-data.frame(efitVEGANCOMPOSITION_W.FD$factors$centroids)#returns the centroids of the factor #2.3.2.3. Resilience of biomass in herbaceous ecosystems #Get data Resilience.FunctionalDiversity.Biomass_H.Correlation<-read.table("ResilienceFDindices-BiomassH.csv", header=TRUE, sep=",", dec=".", na.strings = c("", "NA")) Biomass_H.FD<-Resilience.FunctionalDiversity.Biomass_H.Correlation%>% ungroup()%>% select(Biome,Disturbance.type,LogDisturbance.days,Community.origin,Study.design, Study.type,Resilience.category,LogTemporal.scale) Biomass_H.FD%>%names #Reduce names for simplicity Biomass_H.FD<-Biomass_H.FD %>% rename(DIST=Disturbance.type, DIST.DAYS=LogDisturbance.days, ORI=Community.origin, SD=Study.design, ST=Study.type, RES=Resilience.category, TEMP.SCALE=LogTemporal.scale)%>% mutate(Biome=case_match(Biome, "Boreal Forests/Taiga"~"BF", "Mediterranean Forests, Woodlands and Scrub"~"MF", "Temperate Broadleaf and Mixed Forests"~"TBF", "Temperate Coniferous Forests"~"TCF", "Temperate Grasslands, Savannas and Shrublands"~"TG", "Tropical and subtropical grasslands, savannas, and shrublands"~"TrG", "Tropical and Subtropical Moist Broadleaf Forests "~"TrF", "Tundra"~"Tu"))%>% mutate(DIST=case_match(DIST, "Direct biomass removal"~"Biomass removal", "Dry conditions"~"Dry", "Wet conditions"~"Wet")) #characters need to be transformed to factors in order to use Gower distance function Biomass_H.FD<-Biomass_H.FD%>% mutate(across(where(is_character),as_factor)) #get Gower distances before performing PCoA test.dist.matrixBiomass_H.FD <- gowdis(Biomass_H.FD) #PCoA analysis test.pcoaBiomass_H.FD = capscale(test.dist.matrixBiomass_H.FD ~ 1) SCORESBiomass_H.FD<-scores(test.pcoaBiomass_H.FD) #Scores for each site PCOA1Biomass_H.FD<-SCORESBiomass_H.FD$sites[,1] #Scores for PCoA1 PCOA2Biomass_H.FD<-SCORESBiomass_H.FD$sites[,2] #Scores for PCoA2 #Build a table to plot later PCOABiomass_H.FD<-tibble(PCOA1Biomass_H.FD,PCOA2Biomass_H.FD,Resilience.FunctionalDiversity.Biomass_H.Correlation$ID_Si)%>% rename(PCoA1=PCOA1Biomass_H.FD, PCoA2=PCOA2Biomass_H.FD, ID_Si="Resilience.FunctionalDiversity.Biomass_H.Correlation$ID_Si") #Eigenvalues for each axis eigBiomass_H.FD<-eigenvals(test.pcoaBiomass_H.FD) eigBiomass_H.FDpor<-(eigBiomass_H/sum(eigBiomass_H.FD))*100 #envfit() fits environmental vectors or factors onto and ordination. efitVEGANBiomass_H.FD <- envfit(test.pcoaBiomass_H.FD, Biomass_H.FD, na.rm = TRUE) efitVEGANBiomass_H.FD.vector<-data.frame(efitVEGANBiomass_H.FD$vectors$arrows[c(1:2),c(1:2)])# returns the arrow endpoint of a vector efitVEGANBiomass_H.FD.factor<-data.frame(efitVEGANBiomass_H.FD$factors$centroids)#returns the centroids of the factor ####End of this analysis#### ############################ #3. Meta-analysis of the correlation between Resilience and Richness and Functional Diversity #3.1. Description and objectives of the analysis # For each study site, we calculated Pearson’s correlation coefficients between resilience and the # different functional diversity indices of the monitored vegetation plots. To achieve approximately # normal distributions, we transformed all correlation coefficient values to Fisher’s z-scores # together with the respective measure of effect size variance. # We conducted a series of separated random-effect meta-analyses to calculate the grand # mean effect size (i.e. average Fisher’s z-scores) to test the prediction that resilience would be # positively correlated with 1) a greater variation of trait values (using SR, FR, FDis, FDiv and FEve indices), # and 2) the dominance of fast recovering species (using the small-plant dominance and acquisitive dominance axes). # We included the research site as a random factor to account for the lack of independence between measurements # at the same site. Grand mean effect sizes were considered to be statistically significant if the # approximated 95% confidence intervals (CIs) did not include zero. In addition, we calculated the amount of residual heterogeneity (τ2) # and we visually inspected funnel plots between effect sizes and effect size variances. # We expanded each model by adding the axes of moderator variation as predictor variables. We used an omnibus test of moderators (QM) # to test for the significance of the moderator axes in explaining the heterogeneity of effect sizes between independent sites. # Subgroup-average effect sizes were considered statistically significant if the approximated 95% CIs did not include zero. #3.2. Analysis for the complete data set of Resilience and Species richness #Data Resilience.SpeciesRichness.Correlation ##Getting Fischer Z values for all the sites ResilienceSR.es = escalc(measure="ZCOR", ri=PER_COR, ni = Replicates, data = Resilience.SpeciesRichness.Correlation) #Exclude sites with NAs ResilienceSR.es.NAexc<-ResilienceSR.es[!is.na(ResilienceSR.es$yi),] ResilienceSR.es.NAexc<-ResilienceSR.es.NAexc[!is.na(ResilienceSR.es.NAexc$vi),] #Analysis accounting for within-study dependence (multilevel analysis) of Resilience and SR model.multilevel.analysis.SR = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceSR.es.NAexc) model.multilevel.analysis.SR # forest plots forest(model.multilevel.analysis.SR, slab = ResilienceSR.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.SR) #3.2.1. Analysis per ecosystem type and response category #3.2.1.1. Resilience of biomass in herbaceous ecosystems and Species richness #Data Resilience.SpeciesRichness.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.SpeciesRichness.Biomass_H.CorrelationPCOA<-left_join(Resilience.SpeciesRichness.Biomass_H.Correlation,PCOABiomass_H, by="ID_Si") #geting Fischer Z values for all the sites ResilienceSR.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR, ni = Replicates, data = Resilience.SpeciesRichness.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResilienceSR.Biomass_H.es.NAexc<-ResilienceSR.Biomass_H.es[!is.na(ResilienceSR.Biomass_H.es$yi),] ResilienceSR.Biomass_H.es.NAexc<-ResilienceSR.Biomass_H.es.NAexc[!is.na(ResilienceSR.Biomass_H.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and SR model.multilevel.analysis.SR.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceSR.Biomass_H.es.NAexc) model.multilevel.analysis.SR.Biomass_H # forest plots forest(model.multilevel.analysis.SR.Biomass_H, slab = ResilienceSR.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.SR.Biomass_H) #evalauting PCoA axes model.meta.regr.SR.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceSR.Biomass_H.es.NAexc) model.meta.regr.SR.Biomass_H.PCOA1 model.meta.regr.SR.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceSR.Biomass_H.es.NAexc) model.meta.regr.SR.Biomass_H.PCOA2 #3.2.1.2. Resilience of biomass in woody ecosystems and Species richness #Data Resilience.SpeciesRichness.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.SpeciesRichness.Biomass_W.CorrelationPCOA<-left_join(Resilience.SpeciesRichness.Biomass_W.Correlation,PCOABiomass_W, by="ID_Si") #getting Fischer Z values for all the sites ResilienceSR.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR, ni = Replicates, data = Resilience.SpeciesRichness.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResilienceSR.Biomass_W.es.NAexc<-ResilienceSR.Biomass_W.es[!is.na(ResilienceSR.Biomass_W.es$yi),] ResilienceSR.Biomass_W.es.NAexc<-ResilienceSR.Biomass_W.es.NAexc[!is.na(ResilienceSR.Biomass_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and SR model.multilevel.analysis.SR.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceSR.Biomass_W.es.NAexc) model.multilevel.analysis.SR.Biomass_W # forest plots forest(model.multilevel.analysis.SR.Biomass_W, slab = ResilienceSR.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.SR.Biomass_W) #evalauting PCoA axes model.meta.regr.SR.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceSR.Biomass_W.es.NAexc) model.meta.regr.SR.Biomass_W.PCOA1 model.meta.regr.SR.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceSR.Biomass_W.es.NAexc) model.meta.regr.SR.Biomass_W.PCOA2 #3.2.1.3. Resilience of floristic composition in woody ecosystems and Species richness #Data Resilience.SpeciesRichness.Composition_W.Correlation #Introduce the PCOA values in the table Resilience.SpeciesRichness.Composition_W.CorrelationPCOA<-left_join(Resilience.SpeciesRichness.Composition_W.Correlation,PCOACOMPOSITION_W, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceSR.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR, ni = Replicates, data = Resilience.SpeciesRichness.Composition_W.CorrelationPCOA) #Exclude sites with NAs ResilienceSR.Composition_W.es.NAexc<-ResilienceSR.Composition_W.es[!is.na(ResilienceSR.Composition_W.es$yi),] ResilienceSR.Composition_W.es.NAexc<-ResilienceSR.Composition_W.es.NAexc[!is.na(ResilienceSR.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and SR model.multilevel.analysis.SR.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceSR.Composition_W.es.NAexc) model.multilevel.analysis.SR.Composition_W # forest plots forest(model.multilevel.analysis.SR.Composition_W, slab = ResilienceSR.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.SR.Composition_W) #evalauting PCoA axes model.meta.regr.SR.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceSR.Composition_W.es.NAexc) model.meta.regr.SR.Composition_W.PCOA1 model.meta.regr.SR.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceSR.Composition_W.es.NAexc) model.meta.regr.SR.Composition_W.PCOA2 ####End of this analysis#### ############################ ##3.3. Analysis for the complete data set of Resilience and Functional Richness #Data Resilience.FunctionalDiversity.Correlation #Getting Fischer Z values for all the sites ResilienceFR.es = escalc(measure="ZCOR", ri=PER_COR_FRic, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResilienceFR.es.NAexc<-ResilienceFR.es[!is.na(ResilienceFR.es$yi),] ResilienceFR.es.NAexc<-ResilienceFR.es.NAexc[!is.na(ResilienceFR.es.NAexc$vi),] #analysis accounting for within-study dependence (multilevel analysis) of Resilience and Funtional Richness model.multilevel.analysis.FR = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFR.es.NAexc) model.multilevel.analysis.FR # forest plots forest(model.multilevel.analysis.FR , slab = ResilienceFR.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FR ) #3.3.1. Analysis per ecosystem type and response category #3.3.1.1. Resilience of biomass in herbaceous ecosystems and Functional Richness #Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #geting Fischer Z values for all the sites ResilienceFR.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_FRic, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResilienceFR.Biomass_H.es.NAexc<-ResilienceFR.Biomass_H.es[!is.na(ResilienceFR.Biomass_H.es$yi),] ResilienceFR.Biomass_H.es.NAexc<-ResilienceFR.Biomass_H.es.NAexc[!is.na(ResilienceFR.Biomass_H.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FR model.multilevel.analysis.FR.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFR.Biomass_H.es.NAexc) model.multilevel.analysis.FR.Biomass_H # forest plots forest(model.multilevel.analysis.FR.Biomass_H, slab = ResilienceFR.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FR.Biomass_H) #evalauting PCoA axes model.meta.regr.FR.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFR.Biomass_H.es.NAexc) model.meta.regr.FR.Biomass_H.PCOA1 model.meta.regr.FR.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFR.Biomass_H.es.NAexc) model.meta.regr.FR.Biomass_H.PCOA2 #3.3.1.2. Resilience of biomass in woody ecosystems and Functional Richness #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #getting Fischer Z values for all the sites ResilienceFR.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_FRic, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFR.Biomass_W.es.NAexc<-ResilienceFR.Biomass_W.es[!is.na(ResilienceFR.Biomass_W.es$yi),] ResilienceFR.Biomass_W.es.NAexc<-ResilienceFR.Biomass_W.es.NAexc[!is.na(ResilienceFR.Biomass_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and SR model.multilevel.analysis.FR.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFR.Biomass_W.es.NAexc) model.multilevel.analysis.FR.Biomass_W # forest plots forest(model.multilevel.analysis.FR.Biomass_W, slab = ResilienceFR.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FR.Biomass_W) #evalauting PCoA axes model.meta.regr.FR.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFR.Biomass_W.es.NAexc) model.meta.regr.FR.Biomass_W.PCOA1 model.meta.regr.FR.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFR.Biomass_W.es.NAexc) model.meta.regr.FR.Biomass_W.PCOA2 #3.3.1.3. Resilience of floristic composition in woody ecosystems and Functional Richness #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #getting Fischer Z values for all the sites ResilienceFR.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_FRic, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFR.Composition_W.es.NAexc<-ResilienceFR.Composition_W.es[!is.na(ResilienceFR.Composition_W.es$yi),] ResilienceFR.Composition_W.es.NAexc<-ResilienceFR.Composition_W.es.NAexc[!is.na(ResilienceFR.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FR model.multilevel.analysis.FR.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFR.Composition_W.es.NAexc) model.multilevel.analysis.FR.Composition_W # forest plots forest(model.multilevel.analysis.FR.Composition_W, slab = ResilienceFR.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FR.Composition_W) #evalauting PCoA axes model.meta.regr.FR.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFR.Composition_W.es.NAexc) model.meta.regr.FR.Composition_W.PCOA1 model.meta.regr.FR.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFR.Composition_W.es.NAexc) model.meta.regr.FR.Composition_W.PCOA2 #3.4. Analysis for the complete data set of Resilience and Functional Divergence #Data Resilience.FunctionalDiversity.Correlation #Getting Fischer Z values for all the sites ResilienceFDiv.es = escalc(measure="ZCOR", ri=PER_COR_FDiv, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResilienceFDiv.es.NAexc<-ResilienceFDiv.es[!is.na(ResilienceFDiv.es$yi),] ResilienceFDiv.es.NAexc<-ResilienceFDiv.es.NAexc[!is.na(ResilienceFDiv.es.NAexc$vi),] #Analysis accounting for within-study dependence (multilevel analysis) of Resilience and Funtional Richness model.multilevel.analysis.FDiv = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDiv.es.NAexc) model.multilevel.analysis.FDiv # forest plots forest(model.multilevel.analysis.FDiv , slab = ResilienceFDiv.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDiv ) #3.4.1. Analysis per ecosystem type and response category #3.4.1.1. Resilience of biomass in herbaceous ecosystems and Functional Divergence #Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFDiv.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_FDiv, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResilienceFDiv.Biomass_H.es.NAexc<-ResilienceFDiv.Biomass_H.es[!is.na(ResilienceFDiv.Biomass_H.es$yi),] ResilienceFDiv.Biomass_H.es.NAexc<-ResilienceFDiv.Biomass_H.es.NAexc[!is.na(ResilienceFDiv.Biomass_H.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FDiv.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDiv.Biomass_H.es.NAexc) model.multilevel.analysis.FDiv.Biomass_H # forest plots forest(model.multilevel.analysis.FDiv.Biomass_H, slab = ResilienceFDiv.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDiv.Biomass_H) #evalauting PCoA axes model.meta.regr.FDiv.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDiv.Biomass_H.es.NAexc) model.meta.regr.FDiv.Biomass_H.PCOA1 model.meta.regr.FDiv.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDiv.Biomass_H.es.NAexc) model.meta.regr.FDiv.Biomass_H.PCOA2 #3.4.1.2. Resilience of biomass in woody ecosystems and Functional Divergence #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFDiv.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_FDiv, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFDiv.Biomass_W.es.NAexc<-ResilienceFDiv.Biomass_W.es[!is.na(ResilienceFDiv.Biomass_W.es$yi),] ResilienceFDiv.Biomass_W.es.NAexc<-ResilienceFDiv.Biomass_W.es.NAexc[!is.na(ResilienceFDiv.Biomass_W.es.NAexc$vi),] # Analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FDiv.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDiv.Biomass_W.es.NAexc) model.multilevel.analysis.FDiv.Biomass_W # forest plots forest(model.multilevel.analysis.FDiv.Biomass_W, slab = ResilienceFDiv.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDiv.Biomass_W) #evalauting PCoA axes model.meta.regr.FDiv.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDiv.Biomass_W.es.NAexc) model.meta.regr.FDiv.Biomass_W.PCOA1 model.meta.regr.FDiv.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDiv.Biomass_W.es.NAexc) model.meta.regr.FDiv.Biomass_W.PCOA2 #3.4.1.3. Resilience of floristic composition in woody ecosystems and Functional Divergence #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFDiv.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_FDiv, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFDiv.Composition_W.es.NAexc<-ResilienceFDiv.Composition_W.es[!is.na(ResilienceFDiv.Composition_W.es$yi),] ResilienceFDiv.Composition_W.es.NAexc<-ResilienceFDiv.Composition_W.es.NAexc[!is.na(ResilienceFDiv.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FDiv.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDiv.Composition_W.es.NAexc) model.multilevel.analysis.FDiv.Composition_W # forest plots forest(model.multilevel.analysis.FDiv.Composition_W, slab = ResilienceFDiv.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDiv.Composition_W) #evalauting PCoA axes model.meta.regr.FDiv.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDiv.Composition_W.es.NAexc) model.meta.regr.FDiv.Composition_W.PCOA1 model.meta.regr.FDiv.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDiv.Composition_W.es.NAexc) model.meta.regr.FDiv.Composition_W.PCOA2 #3.5. Analysis for the complete data set of Resilience and Functional Dispersion #Data Resilience.FunctionalDiversity.Correlation ##Getting Fischer Z values for all the sites ResilienceFDis.es = escalc(measure="ZCOR", ri=PER_COR_FDis, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResilienceFDis.es.NAexc<-ResilienceFDis.es[!is.na(ResilienceFDis.es$yi),] ResilienceFDis.es.NAexc<-ResilienceFDis.es.NAexc[!is.na(ResilienceFDis.es.NAexc$vi),] #analysis accounting for within-study dependence (multilevel analysis) of Resilience and Funtional Richness model.multilevel.analysis.FDis = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDis.es.NAexc) model.multilevel.analysis.FDis # forest plots forest(model.multilevel.analysis.FDis , slab = ResilienceFDis.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDis ) #3.5.1. Analysis per ecosystem type and response category #3.5.1.1. Resilience of biomass in herbaceous ecosystems and Functional Dispersion ##Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFDis.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_FDis, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResilienceFDis.Biomass_H.es.NAexc<-ResilienceFDis.Biomass_H.es[!is.na(ResilienceFDis.Biomass_H.es$yi),] ResilienceFDis.Biomass_H.es.NAexc<-ResilienceFDis.Biomass_H.es.NAexc[!is.na(ResilienceFDis.Biomass_H.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FDis.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDis.Biomass_H.es.NAexc) model.multilevel.analysis.FDis.Biomass_H # forest plots forest(model.multilevel.analysis.FDis.Biomass_H, slab = ResilienceFDis.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDis.Biomass_H) #evalauting PCoA axes model.meta.regr.FDis.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDis.Biomass_H.es.NAexc) model.meta.regr.FDis.Biomass_H.PCOA1 model.meta.regr.FDis.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDis.Biomass_H.es.NAexc) model.meta.regr.FDis.Biomass_H.PCOA2 #3.5.1.2. Resilience of biomass in woody ecosystems and Functional Dispersion #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFDis.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_FDis, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFDis.Biomass_W.es.NAexc<-ResilienceFDis.Biomass_W.es[!is.na(ResilienceFDis.Biomass_W.es$yi),] ResilienceFDis.Biomass_W.es.NAexc<-ResilienceFDis.Biomass_W.es.NAexc[!is.na(ResilienceFDis.Biomass_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FDis.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDis.Biomass_W.es.NAexc) model.multilevel.analysis.FDis.Biomass_W # forest plots forest(model.multilevel.analysis.FDis.Biomass_W, slab = ResilienceFDis.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDis.Biomass_W) #evalauting PCoA axes model.meta.regr.FDis.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDis.Biomass_W.es.NAexc) model.meta.regr.FDis.Biomass_W.PCOA1 model.meta.regr.FDis.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDis.Biomass_W.es.NAexc) model.meta.regr.FDis.Biomass_W.PCOA2 #3.5.1.3. Resilience of floristic composition in woody ecosystems and Functional Dispersion #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #getting Fischer Z values for all the sites ResilienceFDis.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_FDis, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFDis.Composition_W.es.NAexc<-ResilienceFDis.Composition_W.es[!is.na(ResilienceFDis.Composition_W.es$yi),] ResilienceFDis.Composition_W.es.NAexc<-ResilienceFDis.Composition_W.es.NAexc[!is.na(ResilienceFDis.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDis model.multilevel.analysis.FDis.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFDis.Composition_W.es.NAexc) model.multilevel.analysis.FDis.Composition_W # forest plots forest(model.multilevel.analysis.FDis.Composition_W, slab = ResilienceFDis.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FDis.Composition_W) #evalauting PCoA axes model.meta.regr.FDis.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFDis.Composition_W.es.NAexc) model.meta.regr.FDis.Composition_W.PCOA1 model.meta.regr.FDis.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFDis.Composition_W.es.NAexc) model.meta.regr.FDis.Composition_W.PCOA2 #3.6. Analysis for the complete data set of Resilience and Functional Evenness #Data Resilience.FunctionalDiversity.Correlation ##Getting Fischer Z values for all the sites ResilienceFEve.es = escalc(measure="ZCOR", ri=PER_COR_FEve, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResilienceFEve.es.NAexc<-ResilienceFEve.es[!is.na(ResilienceFEve.es$yi),] ResilienceFEve.es.NAexc<-ResilienceFEve.es.NAexc[!is.na(ResilienceFEve.es.NAexc$vi),] #analysis accounting for within-study dependence (multilevel analysis) of Resilience and Funtional Richness model.multilevel.analysis.FEve = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFEve.es.NAexc) model.multilevel.analysis.FEve # forest plots forest(model.multilevel.analysis.FEve , slab = ResilienceFEve.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FEve ) #3.6.1. Analysis per ecosystem type and response category #3.6.1.1. Resilience of biomass in herbaceous ecosystems and Functional Evenness #Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFEve.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_FEve, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResilienceFEve.Biomass_H.es.NAexc<-ResilienceFEve.Biomass_H.es[!is.na(ResilienceFEve.Biomass_H.es$yi),] ResilienceFEve.Biomass_H.es.NAexc<-ResilienceFEve.Biomass_H.es.NAexc[!is.na(ResilienceFEve.Biomass_H.es.NAexc$vi),] # Analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FEve.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFEve.Biomass_H.es.NAexc) model.multilevel.analysis.FEve.Biomass_H # forest plots forest(model.multilevel.analysis.FEve.Biomass_H, slab = ResilienceFEve.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FEve.Biomass_H) #evalauting PCoA axes model.meta.regr.FEve.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFEve.Biomass_H.es.NAexc) model.meta.regr.FEve.Biomass_H.PCOA1 model.meta.regr.FEve.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFEve.Biomass_H.es.NAexc) model.meta.regr.FEve.Biomass_H.PCOA2 #3.6.1.2. Resilience of biomass in woody ecosystems and Functional Evenness #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFEve.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_FEve, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFEve.Biomass_W.es.NAexc<-ResilienceFEve.Biomass_W.es[!is.na(ResilienceFEve.Biomass_W.es$yi),] ResilienceFEve.Biomass_W.es.NAexc<-ResilienceFEve.Biomass_W.es.NAexc[!is.na(ResilienceFEve.Biomass_W.es.NAexc$vi),] # Analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.FEve.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFEve.Biomass_W.es.NAexc) model.multilevel.analysis.FEve.Biomass_W # forest plots forest(model.multilevel.analysis.FEve.Biomass_W, slab = ResilienceFEve.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FEve.Biomass_W) #evalauting PCoA axes model.meta.regr.FEve.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFEve.Biomass_W.es.NAexc) model.meta.regr.FEve.Biomass_W.PCOA1 model.meta.regr.FEve.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFEve.Biomass_W.es.NAexc) model.meta.regr.FEve.Biomass_W.PCOA2 #3.6.1.3. Resilience of floristic composition in woody ecosystems and Functional Evenness #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResilienceFEve.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_FEve, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResilienceFEve.Composition_W.es.NAexc<-ResilienceFEve.Composition_W.es[!is.na(ResilienceFEve.Composition_W.es$yi),] ResilienceFEve.Composition_W.es.NAexc<-ResilienceFEve.Composition_W.es.NAexc[!is.na(ResilienceFEve.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDis model.multilevel.analysis.FEve.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFEve.Composition_W.es.NAexc) model.multilevel.analysis.FEve.Composition_W # forest plots forest(model.multilevel.analysis.FEve.Composition_W, slab = ResilienceFEve.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.FEve.Composition_W) #evalauting PCoA axes model.meta.regr.FEve.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResilienceFEve.Composition_W.es.NAexc) model.meta.regr.FEve.Composition_W.PCOA1 model.meta.regr.FEve.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResilienceFEve.Composition_W.es.NAexc) model.meta.regr.FEve.Composition_W.PCOA2 #3.7. Analysis for the complete data set of Resilience and Small plant Dominance #Data Resilience.FunctionalDiversity.Correlation ##Getting Fischer Z values for all the sites ResiliencePC1.es = escalc(measure="ZCOR", ri=PER_COR_PC1, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResiliencePC1.es.NAexc<-ResiliencePC1.es[!is.na(ResiliencePC1.es$yi),] ResilienceFPC1.es.NAexc<-ResiliencePC1.es.NAexc[!is.na(ResiliencePC1.es.NAexc$vi),] #analysis accounting for within-study dependence (multilevel analysis) of Resilience and Funtional Richness model.multilevel.analysis.PC1 = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResilienceFPC1.es.NAexc) model.multilevel.analysis.PC1 # forest plots forest(model.multilevel.analysis.PC1 , slab = ResiliencePC1.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC1 ) #3.7.1. Analysis per ecosystem type and response category #3.7.1.1. Resilience of biomass in herbaceous ecosystems and Small plant Dominance #Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResiliencePC1.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_PC1, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResiliencePC1.Biomass_H.es.NAexc<-ResiliencePC1.Biomass_H.es[!is.na(ResiliencePC1.Biomass_H.es$yi),] ResiliencePC1.Biomass_H.es.NAexc<-ResiliencePC1.Biomass_H.es.NAexc[!is.na(ResiliencePC1.Biomass_H.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and PC1 model.multilevel.analysis.PC1.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC1.Biomass_H.es.NAexc) model.multilevel.analysis.PC1.Biomass_H # forest plots forest(model.multilevel.analysis.PC1.Biomass_H, slab = ResiliencePC1.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC1.Biomass_H) #evalauting PCoA axes model.meta.regr.PC1.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC1.Biomass_H.es.NAexc) model.meta.regr.PC1.Biomass_H.PCOA1 model.meta.regr.PC1.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC1.Biomass_H.es.NAexc) model.meta.regr.PC1.Biomass_H.PCOA2 #3.7.1.2. Resilience of biomass in woody ecosystems and Small plant Dominance #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResiliencePC1.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_PC1, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResiliencePC1.Biomass_W.es.NAexc<-ResiliencePC1.Biomass_W.es[!is.na(ResiliencePC1.Biomass_W.es$yi),] ResiliencePC1.Biomass_W.es.NAexc<-ResiliencePC1.Biomass_W.es.NAexc[!is.na(ResiliencePC1.Biomass_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.PC1.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC1.Biomass_W.es.NAexc) model.multilevel.analysis.PC1.Biomass_W # forest plots forest(model.multilevel.analysis.PC1.Biomass_W, slab = ResiliencePC1.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC1.Biomass_W) #evalauting PCoA axes model.meta.regr.PC1.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC1.Biomass_W.es.NAexc) model.meta.regr.PC1.Biomass_W.PCOA1 model.meta.regr.PC1.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC1.Biomass_W.es.NAexc) model.meta.regr.PC1.Biomass_W.PCOA2 #3.7.1.3. Resilience of floristic composition in woody ecosystems and Small plant Dominance #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #getting Fischer Z values for all the sites ResiliencePC1.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_PC1, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResiliencePC1.Composition_W.es.NAexc<-ResiliencePC1.Composition_W.es[!is.na(ResiliencePC1.Composition_W.es$yi),] ResiliencePC1.Composition_W.es.NAexc<-ResiliencePC1.Composition_W.es.NAexc[!is.na(ResiliencePC1.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and PC1 model.multilevel.analysis.PC1.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC1.Composition_W.es.NAexc) model.multilevel.analysis.PC1.Composition_W # forest plots forest(model.multilevel.analysis.PC1.Composition_W, slab = ResiliencePC1.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC1.Composition_W) #evalauting PCoA axes model.meta.regr.PC1.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC1.Composition_W.es.NAexc) model.meta.regr.PC1.Composition_W.PCOA1 model.meta.regr.PC1.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC1.Composition_W.es.NAexc) model.meta.regr.PC1.Composition_W.PCOA2 #3.8. Analysis for the complete data set of Resilience and Acquisitive plant Dominance #Data Resilience.FunctionalDiversity.Correlation ##Getting Fischer Z values for all the sites ResiliencePC2.es = escalc(measure="ZCOR", ri=PER_COR_PC2, ni = Replicates, data = Resilience.FunctionalDiversity.Correlation) #Exclude sites with NAs ResiliencePC2.es.NAexc<-ResiliencePC2.es[!is.na(ResiliencePC2.es$yi),] ResiliencePC2.es.NAexc<-ResiliencePC2.es.NAexc[!is.na(ResiliencePC2.es.NAexc$vi),] #analysis accounting for within-study dependence (multilevel analysis) of Resilience and Acquisitive Dominance model.multilevel.analysis.PC2 = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC2.es.NAexc) model.multilevel.analysis.PC2 # forest plots forest(model.multilevel.analysis.PC2 , slab = ResiliencePC2.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC2 ) #3.8.1. Analysis per ecosystem type and response category #3.8.1.1. Resilience of biomass in herbaceous ecosystems and Acquisitive plant Dominance #Data Resilience.FunctionalDiversity.Biomass_H.Correlation #Include the PCoA axes calculated in the previous section Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_H.Correlation,PCOABiomass_H.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResiliencePC2.Biomass_H.es = escalc(measure="ZCOR", ri=PER_COR_PC2, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_H.CorrelationPCOA) #Exclude sites with NAs ResiliencePC2.Biomass_H.es.NAexc<-ResiliencePC2.Biomass_H.es[!is.na(ResiliencePC2.Biomass_H.es$yi),] ResiliencePC2.Biomass_H.es.NAexc<-ResiliencePC2.Biomass_H.es.NAexc[!is.na(ResiliencePC2.Biomass_H.es.NAexc$vi),] #Analysis accounting for within-study dependence (multilevel analysis) of Resilience and PC1 model.multilevel.analysis.PC2.Biomass_H = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC2.Biomass_H.es.NAexc) model.multilevel.analysis.PC2.Biomass_H # forest plots forest(model.multilevel.analysis.PC2.Biomass_H, slab = ResiliencePC2.Biomass_H.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC2.Biomass_H) #evaluating PCoA axes model.meta.regr.PC2.Biomass_H.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC2.Biomass_H.es.NAexc) model.meta.regr.PC2.Biomass_H.PCOA1 model.meta.regr.PC2.Biomass_H.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC2.Biomass_H.es.NAexc) model.meta.regr.PC2.Biomass_H.PCOA2 #3.8.1.2. Resilience of biomass in woody ecosystems and Acquisitive plant Dominance #Data Resilience.FunctionalDiversity.Biomass_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Biomass_W.Correlation,PCOABiomass_W.FD, by="ID_Si") #Getting Fischer Z values for all the sites ResiliencePC2.Biomass_W.es = escalc(measure="ZCOR", ri=PER_COR_PC2, ni = Replicates, data = Resilience.FunctionalDiversity.Biomass_W.CorrelationPCOA) #Exclude sites with NAs ResiliencePC2.Biomass_W.es.NAexc<-ResiliencePC2.Biomass_W.es[!is.na(ResiliencePC2.Biomass_W.es$yi),] ResiliencePC2.Biomass_W.es.NAexc<-ResiliencePC2.Biomass_W.es.NAexc[!is.na(ResiliencePC2.Biomass_W.es.NAexc$vi),] # Analysis accounting for within-study dependence (multilevel analysis) of Resilience and FDiv model.multilevel.analysis.PC2.Biomass_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC2.Biomass_W.es.NAexc) model.multilevel.analysis.PC2.Biomass_W # forest plots forest(model.multilevel.analysis.PC2.Biomass_W, slab = ResiliencePC2.Biomass_W.es.NAexc$Site, cex = 0.5) # funnel plot funnel(model.multilevel.analysis.PC2.Biomass_W) #evalauting PCoA axes model.meta.regr.PC2.Biomass_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC2.Biomass_W.es.NAexc) model.meta.regr.PC2.Biomass_W.PCOA1 model.meta.regr.PC2.Biomass_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC2.Biomass_W.es.NAexc) model.meta.regr.PC2.Biomass_W.PCOA2 #3.8.1.3. Resilience of floristic composition in woody ecosystems and Acquisitive plant Dominance #Data Resilience.FunctionalDiversity.Comp_W.Correlation #Introduce the PCOA values in the table Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA<-left_join(Resilience.FunctionalDiversity.Comp_W.Correlation,PCOACOMPOSITION_W.FD, by="ID_Si") #getting Fischer Z values for all the sites ResiliencePC2.Composition_W.es = escalc(measure="ZCOR", ri=PER_COR_PC2, ni = Replicates, data = Resilience.FunctionalDiversity.Comp_W.CorrelationPCOA) #Exclude sites with NAs ResiliencePC2.Composition_W.es.NAexc<-ResiliencePC2.Composition_W.es[!is.na(ResiliencePC2.Composition_W.es$yi),] ResiliencePC2.Composition_W.es.NAexc<-ResiliencePC2.Composition_W.es.NAexc[!is.na(ResiliencePC2.Composition_W.es.NAexc$vi),] # analysis accounting for within-study dependence (multilevel analysis) of Resilience and PC1 model.multilevel.analysis.PC2.Composition_W = rma.mv(yi = yi, V = vi, random = ~ 1|Site, data = ResiliencePC2.Composition_W.es.NAexc) model.multilevel.analysis.PC2.Composition_W # forest plots forest(model.multilevel.analysis.PC2.Composition_W, slab = ResiliencePC2.Composition_W.es.NAexc$Site, cex = 0.5) # funnel plo funnel(model.multilevel.analysis.PC2.Composition_W) #evalauting PCoA axes model.meta.regr.PC2.Composition_W.PCOA1 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA1, data = ResiliencePC2.Composition_W.es.NAexc) model.meta.regr.PC2.Composition_W.PCOA1 model.meta.regr.PC2.Composition_W.PCOA2 = rma.mv(yi = yi, V= vi, random = ~ 1|Site, mods = ~ PCoA2, data = ResiliencePC2.Composition_W.es.NAexc) model.meta.regr.PC2.Composition_W.PCOA2 ####End of this analysis#### ############################ #4. Figures #4.1. Figure 1 - MAP mp <- NULL mapWorld <- borders("world", colour="gray50", fill="gray77") mp <- ggplot() + mapWorld mp <- mp+ geom_point(aes(x=ResilienceSR.es.NAexc$Longitud, y=ResilienceSR.es.NAexc$Latitud), color="black",size = 3) mp <-mp + theme_bw() mp <-mp +theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), axis.title = element_blank()) mp #4.2. Figure 2 - Grand mean effect sizes of the different meta-analyses #graphs for the complete data set #models model.multilevel.analysis.SR model.multilevel.analysis.FR model.multilevel.analysis.FDis model.multilevel.analysis.FDiv model.multilevel.analysis.FEve model.multilevel.analysis.PC1 model.multilevel.analysis.PC2 #data to plot data.to.plot.COMPLETE = data.frame( runningID = 7:1, effect_size = c("Resilience ~ SR", "Resilience ~ FRic", "Resilience ~ FDis", "Resilience ~ FEve", "Resilience ~ FDiv", "Resilience ~ small-plant dominance", "Resilience ~ acquisitive dominance"), yi = c(model.multilevel.analysis.SR$b, model.multilevel.analysis.FR$b, model.multilevel.analysis.FDis$b, model.multilevel.analysis.FEve$b, model.multilevel.analysis.FDiv$b, model.multilevel.analysis.PC1$b, model.multilevel.analysis.PC2$b), ci_lb = c(model.multilevel.analysis.SR$ci.lb, model.multilevel.analysis.FR$ci.lb, model.multilevel.analysis.FDis$ci.lb, model.multilevel.analysis.FEve$ci.lb, model.multilevel.analysis.FDiv$ci.lb, model.multilevel.analysis.PC1$ci.lb, model.multilevel.analysis.PC2$ci.lb), ci_ub = c(model.multilevel.analysis.SR$ci.ub, model.multilevel.analysis.FR$ci.ub, model.multilevel.analysis.FDis$ci.ub, model.multilevel.analysis.FEve$ci.ub, model.multilevel.analysis.FDiv$ci.ub, model.multilevel.analysis.PC1$ci.ub, model.multilevel.analysis.PC2$ci.ub), sample_size= c(model.multilevel.analysis.SR$k, model.multilevel.analysis.FR$k, model.multilevel.analysis.FDis$k, model.multilevel.analysis.FEve$k, model.multilevel.analysis.FDiv$k, model.multilevel.analysis.PC1$k, model.multilevel.analysis.PC2$k), color = c("black", "black", "black", "black", "black", "black", "black")) #setting the labels labs.complete<-c(paste(data.to.plot.COMPLETE$effect_size[1], data.to.plot.COMPLETE$sample_size[1],sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[2], data.to.plot.COMPLETE$sample_size[2],sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[3], data.to.plot.COMPLETE$sample_size[3], sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[4], data.to.plot.COMPLETE$sample_size[4], sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[5], data.to.plot.COMPLETE$sample_size[5], sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[6], data.to.plot.COMPLETE$sample_size[6], sep = "\nn="), paste(data.to.plot.COMPLETE$effect_size[7], data.to.plot.COMPLETE$sample_size[7], sep = "\nn=")) # plotting with ggplot2 ggplot.COMPLETE<-ggplot(data = data.to.plot.COMPLETE) + geom_hline(yintercept = 0, linetype = "dotted") + geom_pointrange(aes(x = runningID, y = yi, ymin = ci_lb, ymax = ci_ub), color = data.to.plot.COMPLETE$color) + scale_x_continuous(breaks = data.to.plot.COMPLETE$runningID, labels = labs.complete) + coord_flip() + ylab("Fishers correlation - Grand mean") + xlab("") + labs(title = " ")+ theme_bw(base_size = 14) + theme(panel.grid = element_blank(), axis.ticks = element_blank()) #4.3. Figure 3 - Grand mean effect sizes of the different meta-analyses grouping by Ecosystem type and Response variable #models model.multilevel.analysis.SR.Biomass_H model.multilevel.analysis.SR.Biomass_W model.multilevel.analysis.SR.Composition_W model.multilevel.analysis.FR.Biomass_H model.multilevel.analysis.FR.Biomass_W model.multilevel.analysis.FR.Composition_W model.multilevel.analysis.FDis.Biomass_H model.multilevel.analysis.FDis.Biomass_W model.multilevel.analysis.FDis.Composition_W model.multilevel.analysis.FDiv.Biomass_H model.multilevel.analysis.FDiv.Biomass_W model.multilevel.analysis.FDiv.Composition_W model.multilevel.analysis.FEve.Biomass_H model.multilevel.analysis.FEve.Biomass_W model.multilevel.analysis.FEve.Composition_W model.multilevel.analysis.PC1.Biomass_H model.multilevel.analysis.PC1.Biomass_W model.multilevel.analysis.PC1.Composition_W model.multilevel.analysis.PC2.Biomass_H model.multilevel.analysis.PC2.Biomass_W model.multilevel.analysis.PC2.Composition_W #4.3.1. Figure 3 - Resilience of biomass in herbaceous ecosystems data.to.ForestPlot.Biomass_H = data.frame( runningID = 7:1, effect_size = c("Resilience ~ SR", "Resilience ~ FRic", "Resilience ~ FDis", "Resilience ~ FDiv", "Resilience ~ FEve", "Resilience ~ small-plant dominance", "Resilience ~ acquisitive dominance"), yi = c(model.multilevel.analysis.SR.Biomass_H$b, model.multilevel.analysis.FR.Biomass_H$b, model.multilevel.analysis.FDis.Biomass_H$b, model.multilevel.analysis.FDiv.Biomass_H$b, model.multilevel.analysis.FEve.Biomass_H$b, model.multilevel.analysis.PC1.Biomass_H$b, model.multilevel.analysis.PC2.Biomass_H$b), ci_lb = c(model.multilevel.analysis.SR.Biomass_H$ci.lb, model.multilevel.analysis.FR.Biomass_H$ci.lb, model.multilevel.analysis.FDis.Biomass_H$ci.lb, model.multilevel.analysis.FDiv.Biomass_H$ci.lb, model.multilevel.analysis.FEve.Biomass_H$ci.lb, model.multilevel.analysis.PC1.Biomass_H$ci.lb, model.multilevel.analysis.PC2.Biomass_H$ci.lb), ci_ub = c(model.multilevel.analysis.SR.Biomass_H$ci.ub, model.multilevel.analysis.FR.Biomass_H$ci.ub, model.multilevel.analysis.FDis.Biomass_H$ci.ub, model.multilevel.analysis.FDiv.Biomass_H$ci.ub, model.multilevel.analysis.FEve.Biomass_H$ci.ub, model.multilevel.analysis.PC1.Biomass_H$ci.ub, model.multilevel.analysis.PC2.Biomass_H$ci.ub), sample_size= c(model.multilevel.analysis.SR.Biomass_H$k, model.multilevel.analysis.FR.Biomass_H$k, model.multilevel.analysis.FDis.Biomass_H$k, model.multilevel.analysis.FDiv.Biomass_H$k, model.multilevel.analysis.FEve.Biomass_H$k, model.multilevel.analysis.PC1.Biomass_H$k, model.multilevel.analysis.PC2.Biomass_H$k), color = c("cyan4", "cyan4", "cyan4", "cyan4", "cyan4", "cyan4", "cyan4")) labsForest.Biomass_H<-c(paste(data.to.ForestPlot.Biomass_H$effect_size[1], data.to.ForestPlot.Biomass_H$sample_size[1],sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[2], data.to.ForestPlot.Biomass_H$sample_size[2],sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[3], data.to.ForestPlot.Biomass_H$sample_size[3], sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[4], data.to.ForestPlot.Biomass_H$sample_size[4], sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[5], data.to.ForestPlot.Biomass_H$sample_size[5], sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[6], data.to.ForestPlot.Biomass_H$sample_size[6], sep = "\nn="), paste(data.to.ForestPlot.Biomass_H$effect_size[7], data.to.ForestPlot.Biomass_H$sample_size[7], sep = "\nn=")) ggplot.Biomass_H<-ggplot(data = data.to.ForestPlot.Biomass_H) + geom_hline(yintercept = 0, linetype = "dotted") + geom_pointrange(aes(x = runningID, y = yi, ymin = ci_lb, ymax = ci_ub), color = data.to.ForestPlot.Biomass_H$color) + scale_x_continuous(breaks = data.to.ForestPlot.Biomass_H$runningID, labels = labsForest.Biomass_H) + coord_flip() + ylab("Fishers correlation - Grand mean") + xlab("") + labs(title = "Biomass_H")+ theme_bw(base_size = 14) + theme(panel.grid = element_blank(), axis.ticks = element_blank()) #4.3.2. Figure 3 - Resilience of biomass in woody ecosystems data.to.ForestPlot.Biomass_W = data.frame( runningID = 7:1, effect_size = c("Resilience ~ SR", "Resilience ~ FRic", "Resilience ~ FDis", "Resilience ~ FDiv", "Resilience ~ FEve", "Resilience ~ small-plant dominance", "Resilience ~ acquisitive dominance"), yi = c(model.multilevel.analysis.SR.Biomass_W$b, model.multilevel.analysis.FR.Biomass_W$b, model.multilevel.analysis.FDis.Biomass_W$b, model.multilevel.analysis.FDiv.Biomass_W$b, model.multilevel.analysis.FEve.Biomass_W$b, model.multilevel.analysis.PC1.Biomass_W$b, model.multilevel.analysis.PC2.Biomass_W$b), ci_lb = c(model.multilevel.analysis.SR.Biomass_W$ci.lb, model.multilevel.analysis.FR.Biomass_W$ci.lb, model.multilevel.analysis.FDis.Biomass_W$ci.lb, model.multilevel.analysis.FDiv.Biomass_W$ci.lb, model.multilevel.analysis.FEve.Biomass_W$ci.lb, model.multilevel.analysis.PC1.Biomass_W$ci.lb, model.multilevel.analysis.PC2.Biomass_W$ci.lb), ci_ub = c(model.multilevel.analysis.SR.Biomass_W$ci.ub, model.multilevel.analysis.FR.Biomass_W$ci.ub, model.multilevel.analysis.FDis.Biomass_W$ci.ub, model.multilevel.analysis.FDiv.Biomass_W$ci.ub, model.multilevel.analysis.FEve.Biomass_W$ci.ub, model.multilevel.analysis.PC1.Biomass_W$ci.ub, model.multilevel.analysis.PC2.Biomass_W$ci.ub), sample_size= c(model.multilevel.analysis.SR.Biomass_W$k, model.multilevel.analysis.FR.Biomass_W$k, model.multilevel.analysis.FDis.Biomass_W$k, model.multilevel.analysis.FDiv.Biomass_W$k, model.multilevel.analysis.FEve.Biomass_W$k, model.multilevel.analysis.PC1.Biomass_W$k, model.multilevel.analysis.PC2.Biomass_W$k), color = c("sienna2", "sienna2", "sienna2", "sienna2", "sienna2", "sienna2", "sienna2")) labsForest.Biomass_W<-c(paste("", data.to.ForestPlot.Biomass_W$sample_size[1],sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[2],sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[3], sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[4], sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[5], sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[6], sep = "\nn="), paste("", data.to.ForestPlot.Biomass_W$sample_size[7], sep = "\nn=")) ggplot.Biomass_W<-ggplot(data = data.to.ForestPlot.Biomass_W) + geom_hline(yintercept = 0, linetype = "dotted") + geom_pointrange(aes(x = runningID, y = yi, ymin = ci_lb, ymax = ci_ub), color = data.to.ForestPlot.Biomass_W$color) + scale_x_continuous(breaks = data.to.ForestPlot.Biomass_W$runningID, labels = labsForest.Biomass_W) + coord_flip() + ylab("Fishers correlation - Grand mean") + xlab("") + labs(title = "Biomass_W")+ theme_bw(base_size = 14) + theme(panel.grid = element_blank(), axis.ticks = element_blank()) #4.3.3. Figure 3 - Resilience of floristic composition in woody ecosystems data.to.ForestPlot.Composition_W = data.frame( runningID = 7:1, effect_size = c("Resilience ~ SR", "Resilience ~ FRic", "Resilience ~ FDis", "Resilience ~ FDiv", "Resilience ~ FEve", "Resilience ~ small-plant dominance", "Resilience ~ acquisitive dominance"), yi = c(model.multilevel.analysis.SR.Composition_W$b, model.multilevel.analysis.FR.Composition_W$b, model.multilevel.analysis.FDis.Composition_W$b, model.multilevel.analysis.FDiv.Composition_W$b, model.multilevel.analysis.FEve.Composition_W$b, model.multilevel.analysis.PC1.Composition_W$b, model.multilevel.analysis.PC2.Composition_W$b), ci_lb = c(model.multilevel.analysis.SR.Composition_W$ci.lb, model.multilevel.analysis.FR.Composition_W$ci.lb, model.multilevel.analysis.FDis.Composition_W$ci.lb, model.multilevel.analysis.FDiv.Composition_W$ci.lb, model.multilevel.analysis.FEve.Composition_W$ci.lb, model.multilevel.analysis.PC1.Composition_W$ci.lb, model.multilevel.analysis.PC2.Composition_W$ci.lb), ci_ub = c(model.multilevel.analysis.SR.Composition_W$ci.ub, model.multilevel.analysis.FR.Composition_W$ci.ub, model.multilevel.analysis.FDis.Composition_W$ci.ub, model.multilevel.analysis.FDiv.Composition_W$ci.ub, model.multilevel.analysis.FEve.Composition_W$ci.ub, model.multilevel.analysis.PC1.Composition_W$ci.ub, model.multilevel.analysis.PC2.Composition_W$ci.ub), sample_size= c(model.multilevel.analysis.SR.Composition_W$k, model.multilevel.analysis.FR.Composition_W$k, model.multilevel.analysis.FDis.Composition_W$k, model.multilevel.analysis.FDiv.Composition_W$k, model.multilevel.analysis.FEve.Composition_W$k, model.multilevel.analysis.PC1.Composition_W$k, model.multilevel.analysis.PC2.Composition_W$k), color = c("sienna2", "sienna2", "sienna2", "sienna2", "sienna2", "sienna2", "sienna2")) labsForest.Composition_W<-c(paste("", data.to.ForestPlot.Composition_W$sample_size[1],sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[2],sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[3], sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[4], sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[5], sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[6], sep = "\nn="), paste("", data.to.ForestPlot.Composition_W$sample_size[7], sep = "\nn=")) ggplot.Composition_W<-ggplot(data = data.to.ForestPlot.Composition_W) + geom_hline(yintercept = 0, linetype = "dotted") + geom_pointrange(aes(x = runningID, y = yi, ymin = ci_lb, ymax = ci_ub), color = data.to.ForestPlot.Composition_W$color) + scale_x_continuous(breaks = data.to.ForestPlot.Composition_W$runningID, labels = labsForest.Composition_W) + coord_flip() + ylab("Fishers correlation - Grand mean") + xlab("") + labs(title = "Composition_W")+ theme_bw(base_size = 14) + theme(panel.grid = element_blank(), axis.ticks = element_blank()) #Multi panel plot ggplot.Biomass_H + ggplot.Biomass_W + ggplot.Composition_W #4.4. Figure 4 - Relationship between the axes of moderator variation and Zscores Herbaceous ecosystems #models model.meta.regr.FDis.Biomass_H.PCOA2 model.meta.regr.PC1.Biomass_H.PCOA1 #4.4.1. Figure 4 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and FDis) #slope to plot slope.FDis.Biomass_H.PCOA2 = model.meta.regr.FDis.Biomass_H.PCOA2$b[2,1] #intercept to plot intercept.FDis.Biomass_H.PCOA2 = model.meta.regr.FDis.Biomass_H.PCOA2$b[1,1] data.to.plot.FDis.Biomass_H.PCOA2 = data.frame( ResilienceFDis.Biomass_H.es.NAexc, ci_lb = predict(model.meta.regr.FDis.Biomass_H.PCOA2, newmods = ResilienceFDis.Biomass_H.es.NAexc$PCoA2, newdata = ResilienceFDis.Biomass_H.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.FDis.Biomass_H.PCOA2, newmods = ResilienceFDis.Biomass_H.es.NAexc$PCoA2,newdata = ResilienceFDis.Biomass_H.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.FDis.Biomass_H.PCOA2<-ggplot() + geom_line(data = data.to.plot.FDis.Biomass_H.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.FDis.Biomass_H.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.FDis.Biomass_H.PCOA2, intercept = intercept.FDis.Biomass_H.PCOA2) + geom_point(data = data.to.plot.FDis.Biomass_H.PCOA2, aes(x = PCoA2, y = yi), colour="cyan4", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.4.2. Figure 4 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and Acquisitive plant dominance) #slope to plot slope.PC1.Biomass_H.PCOA1 = model.meta.regr.PC1.Biomass_H.PCOA1$b[2,1] #intercept to plot intercept.PC1.Biomass_H.PCOA1 = model.meta.regr.PC1.Biomass_H.PCOA1$b[1,1] data.to.plot.PC1.Biomass_H.PCOA1 = data.frame( ResiliencePC1.Biomass_H.es.NAexc, ci_lb = predict(model.meta.regr.PC1.Biomass_H.PCOA1, newmods = ResiliencePC1.Biomass_H.es.NAexc$PCoA1, newdata = ResiliencePC1.Biomass_H.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.PC1.Biomass_H.PCOA1, newmods = ResiliencePC1.Biomass_H.es.NAexc$PCoA1,newdata = ResiliencePC1.Biomass_H.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.PC1.Biomass_H.PCOA1<-ggplot() + geom_line(data = data.to.plot.PC1.Biomass_H.PCOA1, aes(x = PCoA1, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.PC1.Biomass_H.PCOA1, aes(x = PCoA1, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.PC1.Biomass_H.PCOA1, intercept = intercept.PC1.Biomass_H.PCOA1) + geom_point(data = data.to.plot.PC1.Biomass_H.PCOA1, aes(x = PCoA1, y = yi), color="cyan4", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA1)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.5. Figure 5 - Relationship between the axes of moderator variation and Zscores Woody ecosystems #models model.meta.regr.FDiv.Biomass_W.PCOA2 model.meta.regr.PC2.Biomass_W.PCOA2 model.meta.regr.FR.Composition_W.PCOA2 model.meta.regr.FDiv.Composition_W.PCOA2 model.meta.regr.PC2.Composition_W.PCOA2 #4.5.1. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and FDiv) #slope to plot slope.FDiv.Biomass_W.PCOA2 = model.meta.regr.FDiv.Biomass_W.PCOA2$b[2,1] #intercept to plot intercept.FDiv.Biomass_W.PCOA2 = model.meta.regr.FDiv.Biomass_W.PCOA2$b[1,1] data.to.plot.FDiv.Biomass_W.PCOA2 = data.frame( ResilienceFDiv.Biomass_W.es.NAexc, ci_lb = predict(model.meta.regr.FDiv.Biomass_W.PCOA2, newmods = ResilienceFDiv.Biomass_W.es.NAexc$PCoA2, newdata = ResilienceFDiv.Biomass_W.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.FDiv.Biomass_W.PCOA2, newmods = ResilienceFDiv.Biomass_W.es.NAexc$PCoA2,newdata = ResilienceFDiv.Biomass_W.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.FDiv.Biomass_W.PCOA2<-ggplot() + geom_line(data = data.to.plot.FDiv.Biomass_W.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.FDiv.Biomass_W.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.FDiv.Biomass_W.PCOA2, intercept = intercept.FDiv.Biomass_W.PCOA2) + geom_point(data = data.to.plot.FDiv.Biomass_W.PCOA2, aes(x = PCoA2, y = yi), color="sienna2", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.5.2. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of biomass and Acquisitive plant dominance) #slope to plot slope.PC2.Biomass_W.PCOA2 = model.meta.regr.PC2.Biomass_W.PCOA2$b[2,1] #intercept to plot intercept.PC2.Biomass_W.PCOA2 = model.meta.regr.PC2.Biomass_W.PCOA2$b[1,1] data.to.plot.PC2.Biomass_W.PCOA2 = data.frame( ResiliencePC2.Biomass_W.es.NAexc, ci_lb = predict(model.meta.regr.PC2.Biomass_W.PCOA2, newmods = ResiliencePC2.Biomass_W.es.NAexc$PCoA2, newdata = ResiliencePC2.Biomass_W.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.PC2.Biomass_W.PCOA2, newmods = ResiliencePC2.Biomass_W.es.NAexc$PCoA2,newdata = ResiliencePC2.Biomass_W.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.PC2.Biomass_W.PCOA2<-ggplot() + geom_line(data = data.to.plot.PC2.Biomass_W.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.PC2.Biomass_W.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.PC2.Biomass_W.PCOA2, intercept = intercept.PC2.Biomass_W.PCOA2) + geom_point(data = data.to.plot.PC2.Biomass_W.PCOA2, aes(x = PCoA2, y = yi), color="sienna2", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.5.3. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and FRic) #slope to plot slope.FR.Composition_W.PCOA2 = model.meta.regr.FR.Composition_W.PCOA2$b[2,1] #intercept to plot intercept.FR.Composition_W.PCOA2 = model.meta.regr.FR.Composition_W.PCOA2$b[1,1] data.to.plot.FR.Composition_W.PCOA2 = data.frame( ResilienceFR.Composition_W.es.NAexc, ci_lb = predict(model.meta.regr.FR.Composition_W.PCOA2, newmods = ResilienceFR.Composition_W.es.NAexc$PCoA2, newdata = ResilienceFR.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.FR.Composition_W.PCOA2, newmods = ResilienceFR.Composition_W.es.NAexc$PCoA2,newdata = ResilienceFR.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.FR.Composition_W.PCOA2<-ggplot() + geom_line(data = data.to.plot.FR.Composition_W.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.FR.Composition_W.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.FR.Composition_W.PCOA2, intercept = intercept.FR.Composition_W.PCOA2) + geom_point(data = data.to.plot.FR.Composition_W.PCOA2, aes(x = PCoA2, y = yi), color="sienna2", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.5.4. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and FDiv) #slope to plot slope.FDiv.Composition_W.PCOA2 = model.meta.regr.FDiv.Composition_W.PCOA2$b[2,1] #intercept to plot intercept.FDiv.Composition_W.PCOA2 = model.meta.regr.FDiv.Composition_W.PCOA2$b[1,1] data.to.plot.FDiv.Composition_W.PCOA2 = data.frame( ResilienceFDiv.Composition_W.es.NAexc, ci_lb = predict(model.meta.regr.FDiv.Composition_W.PCOA2, newmods = ResilienceFDiv.Composition_W.es.NAexc$PCoA2, newdata = ResilienceFDiv.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.FDiv.Composition_W.PCOA2, newmods = ResilienceFDiv.Composition_W.es.NAexc$PCoA2,newdata = ResilienceFDiv.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.FDiv.Composition_W.PCOA2<-ggplot() + geom_line(data = data.to.plot.FDiv.Composition_W.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.FDiv.Composition_W.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.FDiv.Composition_W.PCOA2, intercept = intercept.FDiv.Composition_W.PCOA2) + geom_point(data = data.to.plot.FDiv.Composition_W.PCOA2, aes(x = PCoA2, y = yi), color="sienna2", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.5.5. Figure 5 - Regression between Fisher Z score and Moderators variation (correlation between resilience of composition and Acquisitive plant dominance) #slope to plot slope.PC2.Composition_W.PCOA2 = model.meta.regr.PC2.Composition_W.PCOA2$b[2,1] #intercept to plot intercept.PC2.Composition_W.PCOA2 = model.meta.regr.PC2.Composition_W.PCOA2$b[1,1] data.to.plot.PC2.Composition_W.PCOA2 = data.frame( ResiliencePC2.Composition_W.es.NAexc, ci_lb = predict(model.meta.regr.PC2.Composition_W.PCOA2, newmods = ResiliencePC2.Composition_W.es.NAexc$PCoA2, newdata = ResiliencePC2.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.lb, ci_ub = predict(model.meta.regr.PC2.Composition_W.PCOA2, newmods = ResiliencePC2.Composition_W.es.NAexc$PCoA2,newdata = ResiliencePC2.Composition_W.es.NAexc, interval="confidence", level = 0.95)$ci.ub) #plotting with ggplot2 plot.PC2.Composition_W.PCOA2<-ggplot() + geom_line(data = data.to.plot.PC2.Composition_W.PCOA2, aes(x = PCoA2, y= ci_lb), linetype = "dotted") + geom_line(data = data.to.plot.PC2.Composition_W.PCOA2, aes(x = PCoA2, y= ci_ub), linetype = "dotted") + geom_abline(slope = slope.PC2.Composition_W.PCOA2, intercept = intercept.PC2.Composition_W.PCOA2) + geom_point(data = data.to.plot.PC2.Composition_W.PCOA2, aes(x = PCoA2, y = yi), color="sienna2", alpha = 0.3, cex=4) + xlab("Moderator variation (PCoA2)") + ylab("Fisher`s Z-score") + theme_bw() + theme(legend.position = "none") #4.6. Appendix 4 - PCoA Data set on resilience and Species Richness #4.6.1. Appendix 4 - Complete dataset data.to.plot.complete<-rbind(efitVEGANcomplete.vector,efitVEGANcomplete.factor) labs<-c(glue("PCoA1({round(eigCOMPLETEpor[1], digits=2)}%)"), glue("PCoA2({round(eigCOMPLETEpor[2], digits=2)}%)")) ggplotpcoa2<-ggplot(data=PCOAcomplete, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4,aes(color=Resilience.SpeciesRichness.Correlation$Ecosystem.type))+ scale_color_manual(values = c("Herbaceous" = "cyan4","Woody"="sienna2"))+ geom_segment(data = data.to.plot.complete, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.complete, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.complete)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x=labs[1],y=labs[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.6.2. Appendix 4 - Woody ecosystems data.to.plot.WOODY<-rbind(efitVEGANWOODY.vector,efitVEGANWOODY.factor) labsWoody<-c(glue("PCoA1({round(eigWOODYpor[1], digits=2)}%)"), glue("PCoA2({round(eigWOODYpor[2], digits=2)}%)")) ggplotpcoaWOODY2<-ggplot(data=PCOAWOODY, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="sienna2")+ geom_segment(data = data.to.plot.WOODY, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.WOODY, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.WOODY)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsWoody[1],y=labsWoody[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.6.3. Appendix 4 - Herbaceous ecosystems data.to.plot.HERBACEOUS<-rbind(efitVEGANHERBACEOUS.vector,efitVEGANHERBACEOUS.factor) labsnowoody<-c(glue("PCoA1({round(eigHERBACEOUSpor[1], digits=2)}%)"), glue("PCoA2({round(eigHERBACEOUSpor[2], digits=2)}%)")) ggplotpcoaHERBACEOUS2<-ggplot(data=PCOAHERBACEOUS, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="cyan4")+ geom_segment(data = data.to.plot.HERBACEOUS, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.HERBACEOUS, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.HERBACEOUS)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsnowoody[1],y=labsnowoody[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.6.4. Appendix 4 - Resilience of biomass in woody ecosystems data.to.plot.Biomass_W<-rbind(efitVEGANBiomass_W.vector,efitVEGANBiomass_W.factor) labs.Biomass_W<-c(glue("PCoA1({round(eigBiomass_Wpor[1], digits=2)}%)"), glue("PCoA2({round(eigBiomass_Wpor[2], digits=2)}%)")) ggplotpcoaBiomass_W2<-ggplot(data=PCOABiomass_W, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4,color="sienna2")+ geom_segment(data = data.to.plot.Biomass_W, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.Biomass_W, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.Biomass_W)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labs.Biomass_W[1],y=labs.Biomass_W[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.6.5. Appendix 4 - Resilience of floristic composition in woody ecosystems data.to.plot.COMPOSITION_W<-rbind(efitVEGANCOMPOSITION_W.vector,efitVEGANCOMPOSITION_W.factor) labs.COMPOSITION_W<-c(glue("PCoA1({round(eigCOMPOSITION_Wpor[1], digits=2)}%)"), glue("PCoA2({round(eigCOMPOSITION_Wpor[2], digits=2)}%)")) ggplotpcoaCOMPOSITION_W2<-ggplot(data=PCOACOMPOSITION_W, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="sienna2")+ geom_segment(data = data.to.plot.COMPOSITION_W, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.COMPOSITION_W, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.COMPOSITION_W)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labs.COMPOSITION_W[1],y=labs.COMPOSITION_W[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.6.6. Appendix 4 - Resilience of biomass herbaceous ecosystems data.to.plot.Biomass_H<-rbind(efitVEGANBiomass_H.vector,efitVEGANBiomass_H.factor) labsBiomass_H<-c(glue("PCoA1({round(eigBiomass_Hpor[1], digits=2)}%)"), glue("PCoA2({round(eigBiomass_Hpor[2], digits=2)}%)")) ggplotpcoaBiomass_H2<-ggplot(data=PCOABiomass_H, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="cyan4")+ geom_segment(data = data.to.plot.Biomass_H, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.Biomass_H, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.Biomass_H)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsBiomass_H[1],y=labsBiomass_H[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7. Appendix 5 - PCoA Data set on resilience and Functional Diverstity indices #4.7.1. Appendix 5 - Complete dataset data.to.plot.complete.FD<-rbind(efitVEGANcomplete.FD.vector,efitVEGANcomplete.FD.factor) labs.FD<-c(glue("PCoA1({round(eigCOMPLETE.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigCOMPLETE.FDpor[2], digits=2)}%)")) ggplotpcoa.FD2<-ggplot(data=PCOAcomplete.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4,aes(color=Resilience.FunctionalDiversity.Correlation$Ecosystem.type))+ scale_color_manual(values = c("Herbaceous" = "cyan4","Woody"="sienna2"))+ geom_segment(data = data.to.plot.complete.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.complete.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.complete.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x=labs.FD[1],y=labs.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7.2. Appendix 5 - Woody ecosystems data.to.plot.WOODY.FD<-rbind(efitVEGANWOODY.FD.vector,efitVEGANWOODY.FD.factor) labsWoody.FD<-c(glue("PCoA1({round(eigWOODY.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigWOODY.FDpor[2], digits=2)}%)")) ggplotpcoaWOODY.FD2<-ggplot(data=PCOAWOODY.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="sienna2")+ geom_segment(data = data.to.plot.WOODY.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.WOODY.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.WOODY.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsWoody.FD[1],y=labsWoody.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7.3. Appendix 5 - Herbaceous ecosystems data.to.plot.HERBACEOUS.FD<-rbind(efitVEGANHERBACEOUS.FD.vector,efitVEGANHERBACEOUS.FD.factor) labsherbaceous.FD<-c(glue("PCoA1({round(eigHERBACEOUS.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigHERBACEOUS.FDpor[2], digits=2)}%)")) ggplotpcoaHERBACEOUS.FD2<-ggplot(data=PCOAHERBACEOUS.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="cyan4")+ geom_segment(data = data.to.plot.HERBACEOUS.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.HERBACEOUS.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.HERBACEOUS.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsherbaceous.FD[1],y=labsherbaceous.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7.4. Appendix 5 - Resilience of biomass in woody ecosystems data.to.plot.Biomass_W.FD<-rbind(efitVEGANBiomass_W.FD.vector,efitVEGANBiomass_W.FD.factor) labs.Biomass_W.FD<-c(glue("PCoA1({round(eigBiomass_W.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigBiomass_W.FDpor[2], digits=2)}%)")) ggplotpcoaBiomass_W.FD2<-ggplot(data=PCOABiomass_W.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4,color="sienna2")+ geom_segment(data = data.to.plot.Biomass_W.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.Biomass_W.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.Biomass_W.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labs.Biomass_W.FD[1],y=labs.Biomass_W.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7.5. Appendix 5 - Resilience of floristic compoitionin woody ecosystems data.to.plot.COMPOSITION_W.FD<-rbind(efitVEGANCOMPOSITION_W.FD.vector,efitVEGANCOMPOSITION_W.FD.factor) labs.COMPOSITION_W.FD<-c(glue("PCoA1({round(eigCOMPOSITION_W.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigCOMPOSITION_W.FDpor[2], digits=2)}%)")) ggplotpcoaCOMPOSITION_W.FD2<-ggplot(data=PCOACOMPOSITION_W.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="sienna2")+ geom_segment(data = data.to.plot.COMPOSITION_W.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.COMPOSITION_W.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.COMPOSITION_W.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labs.COMPOSITION_W.FD[1],y=labs.COMPOSITION_W.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") #4.7.6. Appendix 5 - Resilience of biomass herbaceous ecosystems data.to.plot.Biomass_H.FD<-rbind(efitVEGANBiomass_H.FD.vector,efitVEGANBiomass_H.FD.factor) labsBiomass_H.FD<-c(glue("PCoA1({round(eigBiomass_H.FDpor[1], digits=2)}%)"), glue("PCoA2({round(eigBiomass_H.FDpor[2], digits=2)}%)")) ggplotpcoaBiomass_H.FD2<-ggplot(data=PCOABiomass_H.FD, aes(x=PCoA1, y=PCoA2))+ geom_point(cex=4, color="cyan4")+ geom_segment(data = data.to.plot.Biomass_H.FD, aes(x = 0, xend=MDS1, y=0, yend=MDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", linewidth=0.5) + ggrepel::geom_text_repel(data = data.to.plot.Biomass_H.FD, aes(x=MDS1, y=MDS2, label = rownames(data.to.plot.Biomass_H.FD)), cex = 4, direction = "both", segment.size = 0.25, max.overlaps = 18)+ labs(x= labsBiomass_H.FD[1],y=labsBiomass_H.FD[2])+ labs(title = "Ordination with moderators")+ theme_bw()+ theme(legend.position = "non") ####End of figures#### ######################