Intro

This is a step-by-step walk-through of the analysis of the GermanBIODEPTH site in the paper to illustrate the different analyses weconducted to evaluate the relationship between biodiversity andecosystem multifunctionality.

We begin by loading themultifunc package. If it is notcurrently installed, it can be found on github, and installed using thedevtools library.

library(devtools)install_github("multifunc",username="jebyrnes")

library(multifunc)#for plottinglibrary(ggplot2)library(patchwork)#for datalibrary(tidyr)library(dplyr)library(purrr)library(forcats)#for analysuslibrary(car)

Next, we load in the BIODEPTH data. We use some of the helperfunctions in the package to identify the column numbers of the relevantfunctions for later use.

#Read in data  to run sample analyses on the biodepth datadata(all_biodepth)allVars<-c("biomassY3","root3","N.g.m2","light3","N.Soil","wood3","cotton3")varIdx<-which(names(all_biodepth)%in% allVars)

As we want to work with just the German data, we next subset down thedata to just Germany. We then look at the whether species were seeded ornot, and determine which species were indeed seeded in Germany. If acolumn of species is full of zeroes, then it was not seeded in theseplots. We will need the column indices of relevant species for later usein the overlap analyses.

######## Now, specify what data we're working with - Germany# and the relevant variables we'll be working with#######germany<- all_biodepth%>%filter(location=="Germany")# which variables have > 2/3 of the values not NA?german_vars<-whichVars(germany, allVars)# What are the names of species in this dataset# that have at least some values > 0?species<-relevantSp(germany,26:ncol(germany))spIDX<-which(names(germany)%in% species)#in case we need these

Single Function Approach

First, we will demonstrate the qualitative single function approach.As we’re usingggplot2 for plotting, we’ll usetidyrto pivot the data long for a nice faceted plot of each relationship.

germanyForPlotting<- germany%>%select(Diversity, german_vars)%>%pivot_longer(cols = german_vars,names_to ="variable",values_to ="value")%>%mutate(variable =factor(variable))#make the levels of the functions into something nice for plotslevels(germanyForPlotting$variable)<-c('Aboveground Biomass','Root Biomass','Cotton Decomposition','Soil Nitrogen','Plant Nitrogen')

Nest, as we want to display additional information about the fit ofdiversity to each function, we’ll need to iterate over each function,derive a fit, and then make some labels to be used in the eventualplot.

germanyLabels<- germanyForPlotting%>%group_by(variable)%>%nest()%>%mutate(fits =map(data,~lm(value~ Diversity,data=.x)),stats =map(fits, broom::glance))%>%unnest(stats)%>%ungroup()%>%mutate(labels =paste0("p = ",round(p.value,3),", ",expression(R^2)," = ",round(r.squared,2)),labels =gsub("p = 0 ","p < 0.001 ", labels),Diversity =7,value=c(2000,1200,1.5,6,20))

Finally, we will put it all together into a single plot.

ggplot(aes(x=Diversity,y=value),data=germanyForPlotting)+geom_point(size=3)+facet_wrap(~variable,scales="free")+theme_bw(base_size=15)+stat_smooth(method="lm",colour="black",size=2)+xlab("Species Richness")+ylab("Value of Function")+geom_text(data=germanyLabels,aes(label=labels))+theme(panel.grid =element_blank())

Next, we ask whether the best performing species differs betweenfunctions. To do that, we have to pull out all monoculture plots, thenreshape the data and run a quick qualitative analysis on it to determinethe best performer.

#pull out only those plots planted in monoculturemonoGermany<- germany%>%filter(rowSums(select(., ACHMIL1:VICTET1))==1)%>%#get mono namespivot_longer(ACHMIL1:VICTET1,names_to ="mono")%>%filter(value!=0)%>%select(-value)%>%#pivot by functionpivot_longer(german_vars,names_to ="variable")%>%#get the mean for each monoculturegroup_by(variable, mono)%>%summarize(value =mean(value))#now, who is the best performer?monoGermany%>%group_by(variable)%>%filter(value==max(value))

# # A tibble: 5 × 3# # Groups:   variable [5]#   variable  mono       value#   <chr>     <chr>      <dbl># 1 N.Soil    TRIREP1    4.81 # 2 N.g.m2    TRIPRA1   16.9  # 3 biomassY3 TRIPRA1  818.   # 4 cotton3   TRIREP1    0.854# 5 root3     FESRUB1 1072.

Overlap Approach

One of the functions we are looking at, soil N, first needs to bereflected before we can apply the overlap approach.

germany<- germany%>%mutate(N.Soil =-1*N.Soil+max(N.Soil,na.rm=T))

We can examine any single species and generate coefficient values and‘effects’ - i.e., whether a species has a positive (1), negative (-1),or neutral (0) effect on a function. Let’s examine this for biomassproduction in year 3.

spList<-sAICfun("biomassY3", species, germany)spList

# $pos.sp# [1] "DACGLO1" "FESPRA1" "LATPRA1" "LOTCOR1" "TRIPRA1" "TRIREP1"# # $neg.sp# [1] "CAMPAT1" "PHLPRA1" "RANACR1"# # $neu.sp#  [1] "ACHMIL1" "ALOPRA1" "ANTODO1" "ARRELA1" "BROHOR1" "CENJAC1" "CHRLEU1"#  [8] "CREBIE1" "CYNCRI1" "FESRUB1" "GERPRA1" "HOLLAN1" "KNAARV1" "LEOAUT1"# [15] "LOLPER1" "LYCFLO1" "PIMMAJ1" "PLALAN1" "RUMACE1" "TAROFF1" "VICCRA1"# [22] "VICSEP1"# # $functions# [1] "biomassY3"# # $coefs#     ACHMIL1     ALOPRA1     ANTODO1     ARRELA1     BROHOR1     CAMPAT1 #      0.0000      0.0000      0.0000      0.0000      0.0000   -391.5405 #     CENJAC1     CHRLEU1     CREBIE1     CYNCRI1     DACGLO1     FESPRA1 #      0.0000      0.0000      0.0000      0.0000    104.6404    135.6735 #     FESRUB1     GERPRA1     HOLLAN1     KNAARV1     LATPRA1     LEOAUT1 #      0.0000      0.0000      0.0000      0.0000    248.5566      0.0000 #     LOLPER1     LOTCOR1     LYCFLO1     PHLPRA1     PIMMAJ1     PLALAN1 #      0.0000    375.6694      0.0000   -297.7452      0.0000      0.0000 #     RANACR1     RUMACE1     TAROFF1     TRIPRA1     TRIREP1     VICCRA1 #   -140.2996      0.0000      0.0000    492.6311    325.0173      0.0000 #     VICSEP1 (Intercept) #      0.0000    335.8382 # # $effects# ACHMIL1 ALOPRA1 ANTODO1 ARRELA1 BROHOR1 CAMPAT1 CENJAC1 CHRLEU1 CREBIE1 CYNCRI1 #       0       0       0       0       0      -1       0       0       0       0 # DACGLO1 FESPRA1 FESRUB1 GERPRA1 HOLLAN1 KNAARV1 LATPRA1 LEOAUT1 LOLPER1 LOTCOR1 #       1       1       0       0       0       0       1       0       0       1 # LYCFLO1 PHLPRA1 PIMMAJ1 PLALAN1 RANACR1 RUMACE1 TAROFF1 TRIPRA1 TRIREP1 VICCRA1 #       0      -1       0       0      -1       0       0       1       1       0 # VICSEP1 #       0

Using this as our starting point, we can use the getRedundancyfunction to generate a) an effect matrix for all functions, b) acoefficient matrix for all functions, and c) a standardized coefficientmatrix for all functions.

redund<-getRedundancy(german_vars, species, germany)coefs<-getRedundancy(german_vars, species, germany,output="coef")stdCoefs<-stdEffects(coefs, germany, german_vars, species)#for exampleredund

#           ACHMIL1 ALOPRA1 ANTODO1 ARRELA1 BROHOR1 CAMPAT1 CENJAC1 CHRLEU1# biomassY3       0       0       0       0       0      -1       0       0# root3           1       0      -1       0       1       0       0       0# N.g.m2          1       0       0      -1       0       0       0       0# N.Soil         -1       0      -1       0       1      -1       0       0# cotton3         0       1       1       0      -1       1       0       0#           CREBIE1 CYNCRI1 DACGLO1 FESPRA1 FESRUB1 GERPRA1 HOLLAN1 KNAARV1# biomassY3       0       0       1       1       0       0       0       0# root3           0       0       0       0       1       0      -1       0# N.g.m2         -1       0       0       1       0       0      -1       0# N.Soil          1       0       1      -1       0       0       0       0# cotton3         0       0       0       0       0       0       0       0#           LATPRA1 LEOAUT1 LOLPER1 LOTCOR1 LYCFLO1 PHLPRA1 PIMMAJ1 PLALAN1# biomassY3       1       0       0       1       0      -1       0       0# root3           1       0       1       0       0       0       0      -1# N.g.m2          0       0       1       1       0      -1       0      -1# N.Soil          1       0       0       1       0       1       0       0# cotton3        -1       0      -1       0       0       0       0       1#           RANACR1 RUMACE1 TAROFF1 TRIPRA1 TRIREP1 VICCRA1 VICSEP1# biomassY3      -1       0       0       1       1       0       0# root3          -1       0       0      -1       0       0       0# N.g.m2         -1       0       0       1       1       0       0# N.Soil          1       0       0      -1      -1       0       0# cotton3         0       0       0       1       1       0       0

Using the effect matrix, we then estimate the average number ofspecies affecting all possible combinations of functions - both positiveand negative - and then plot the results.

#plot the num. functions by fraction of the species pool neededposCurve<-divNeeded(redund,type="positive")posCurve$div<-posCurve$div/ncol(redund)pc<-qplot(nfunc, div,data=posCurve,group=nfunc,geom=c("boxplot"))+geom_jitter(size=4,position =position_jitter(height=0.001,width = .04))+ylab("Fraction of Species Pool\nwith Positive Effect\n")+xlab("\nNumber of Functions")+theme_bw(base_size=24)+ylim(c(0,1))negCurve<-divNeeded(redund,type="negative")negCurve$div<-negCurve$div/ncol(redund)nc<-qplot(nfunc, div,data=negCurve,group=nfunc,geom=c("boxplot"))+geom_jitter(size=4,position =position_jitter(height=0.001,width = .04))+ylab("Fraction of Species Pool\nwith Negative Effect\n")+xlab("\nNumber of Functions")+theme_bw(base_size=24)+ylim(c(0,1))#combine these into one plotpc+ nc+plot_annotation(tag_levels ='A')

We could have also looked at overlap indices, such as the Sorensonindex.

allOverlapPos<-map_df(2:nrow(redund),~getOverlapSummary(redund,m = .x,denom ="set"),.id ="m")allOverlapPos

# # A tibble: 4 × 4#   m     meanOverlap sdOverlap     n#   <chr>       <dbl>     <dbl> <dbl># 1 1          0.283      0.206    10# 2 2          0.0658     0.116    10# 3 3          0          0         5# 4 4          0         NA         1

One we’ve examined the consequences of functional overlap, we canthen move on to compare the number and size of these effects. First,what is the relative balance of positive to negative contributions foreach species?

posNeg<-data.frame(Species =colnames(redund),Pos =colSums(filterOverData(redund)),Neg =colSums(filterOverData(redund,type="negative")))#plot itggplot(aes(x=Pos,y= Neg),data=posNeg)+geom_jitter(position =position_jitter(width =0.05,height =0.05),size=5,shape=1)+theme_bw(base_size=18)+xlab("\n# of Positive Contributions per species\n")+ylab("# of Negative Contributions per species\n")+annotate("text",x=0,y=3,label="a)")+stat_smooth(method="glm",colour="black",size=2,method.args =list(family=poisson(link="identity")))+geom_abline(slope =1,intercept =0,lty =2,color ="red")

So, a generally positive relationship. We can see from the figurethat the slope is less than 1:1, so positive effects accumulate faster.Still, more positive effects = more negative effects. What about thesize of those positive versus negative effects?

#now get the standardized effect sizesposNeg<-within(posNeg, {  stdPosMean<-colSums(filterCoefData(stdCoefs))/Pos  stdPosMean[which(is.nan(stdPosMean))]<-0  stdNegMean<-colSums(filterCoefData(stdCoefs,type="negative"))/Neg  stdNegMean[which(is.nan(stdNegMean))]<-0})ggplot(aes(x=stdPosMean,y= stdNegMean),data=posNeg)+#  geom_point(size=3) +geom_jitter(position =position_jitter(width = .02,height = .02),size=5,shape=1)+theme_bw(base_size=18)+xlab("\nAverage Standardized Size of\nPositive Contributions")+ylab("Average Standardized Size of\nNegative Contributions\n")+stat_smooth(method="lm",colour="black",size=2)+annotate("text",x=0,y=0.25,label="b)")+geom_abline(slope=-1,intercept=0,size=1,colour="black",lty=3)

Averaging Approach

Now we’ll try the averaging approach. The first step is to create aset of columns where we have standardized all of the functions ofinterest, and then create an average of those standardizedfunctions.

#add on the new functions along with the averaged multifunctional indexgermany<-cbind(germany,getStdAndMeanFunctions(germany, german_vars))#germany<-cbind(germany, getStdAndMeanFunctions(germany, vars, standardizeZScore))

We can then plot the averaged multifunctionality quite simply

#plot itggplot(aes(x=Diversity,y=meanFunction),data=germany)+geom_point(size=3)+theme_bw(base_size=15)+stat_smooth(method="lm",colour="black",size=2)+xlab("\nSpecies Richness")+ylab("Average Value of Standardized Functions\n")

We may want to look at all of the functions on a standardized scale,and add in the averaged line for comparison as well. We can do this byreshaping the data, and plotting it.

#reshape for plotting everything with ggplot2germanyMeanForPlotting<- germany%>%select(Diversity, biomassY3.std:meanFunction)%>%pivot_longer(cols =c(biomassY3.std:meanFunction),names_to ="variable")%>%mutate(variable =fct_inorder(variable))#nice names for plottinglevels(germanyMeanForPlotting$variable)<-c('Aboveground Biomass','Root Biomass','Cotton Decomposition','Soil Nitrogen','Plant Nitrogen','Mean Multifuncion Index')#plot itggplot(aes(x=Diversity,y=value),data=germanyMeanForPlotting)+geom_point(size=3)+facet_grid(~variable)+theme_bw(base_size=15)+stat_smooth(method="lm",colour="black",size=2)+xlab("\nSpecies Richness")+ylab("Standardized Value of Function\n")

Last, let’s test the statistical fit of the effect of diversity onthe averaged multifunctionality index.

#statistical fitaveFit<-lm(meanFunction~ Diversity,data=germany)Anova(aveFit)

# Anova Table (Type II tests)# # Response: meanFunction#            Sum Sq Df F value    Pr(>F)    # Diversity 0.29502  1   32.63 4.042e-07 ***# Residuals 0.52440 58                      # ---# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aveFit)

# # Call:# lm(formula = meanFunction ~ Diversity, data = germany)# # Residuals:#       Min        1Q    Median        3Q       Max # -0.200533 -0.069321  0.006506  0.062717  0.288942 # # Coefficients:#             Estimate Std. Error t value Pr(>|t|)    # (Intercept) 0.385422   0.017049  22.606  < 2e-16 ***# Diversity   0.015360   0.002689   5.712 4.04e-07 ***# ---# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # Residual standard error: 0.09509 on 58 degrees of freedom# Multiple R-squared:   0.36,   Adjusted R-squared:  0.349 # F-statistic: 32.63 on 1 and 58 DF,  p-value: 4.042e-07

Threshold Approach

First, we need to create a new data set that looks at the number offunctions greater than or equal to a threshold across a wide range ofthresholds.

germanyThresh<-getFuncsMaxed(germany, german_vars,threshmin=0.05,threshmax=0.99,prepend=c("plot","Diversity"),maxN=7)

Next, let’s perform an analysis on just the 0.8 threshold data.

mfuncGermanyLinear08<-glm(funcMaxed~ Diversity,data=subset(germanyThresh, germanyThresh$thresholds=="0.8"),family=quasipoisson(link="identity"))Anova(mfuncGermanyLinear08,test.statistic="F")

# Analysis of Deviance Table (Type II tests)# # Response: funcMaxed# Error estimate based on Pearson residuals # #           Sum Sq Df F value Pr(>F)  # Diversity  6.375  1  6.4928 0.0135 *# Residuals 56.944 58                 # ---# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mfuncGermanyLinear08)

# # Call:# glm(formula = funcMaxed ~ Diversity, family = quasipoisson(link = "identity"), #     data = subset(germanyThresh, germanyThresh$thresholds == #         "0.8"))# # Deviance Residuals: #     Min       1Q   Median       3Q      Max  # -1.4741  -0.9796  -0.8101   0.8364   2.0621  # # Coefficients:#             Estimate Std. Error t value Pr(>|t|)  # (Intercept)  0.27756    0.11302   2.456   0.0171 *# Diversity    0.05056    0.02476   2.042   0.0457 *# ---# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # (Dispersion parameter for quasipoisson family taken to be 0.9817929)# #     Null deviance: 66.909  on 59  degrees of freedom# Residual deviance: 60.534  on 58  degrees of freedom# AIC: NA# # Number of Fisher Scoring iterations: 6

To get a better sense of how robust these results are to thresholdchoice, let’s plot the diversity versus number of functions greater thanor equal to a threshold for four different thresholds.

gcPlot<-subset(germanyThresh, germanyThresh$thresholds%in%qw(0.2,0.4,0.6,0.8))#note, using qw as %in% is a string comparison operatorgcPlot$percent<-paste(100*gcPlot$thresholds,"%",sep="")qplot(Diversity, funcMaxed,data=gcPlot,facets=~percent)+stat_smooth(method="glm",method.args =list(family=quasipoisson(link="identity")),colour="red",lwd=1.2)+ylab(expression("Number of Functions">= Threshold))+xlab("Species Richness")+theme_bw(base_size=14)+geom_text(data=data.frame(percent =unique(gcPlot$percent),lab =paste(letters[1:4],")",sep=""),Diversity=2,funcMaxed=6                       ),mapping=aes(x=Diversity,y=funcMaxed,label=lab))

Given that variation, let’s look at the entire spread ofthresholds.

germanyThresh$percent<-100*germanyThresh$thresholdsggplot(data=germanyThresh,aes(x=Diversity,y=funcMaxed,group=percent))+ylab(expression("Number of Functions">= Threshold))+xlab("Species Richness")+stat_smooth(method="glm",method.args =list(family=quasipoisson(link="identity")),lwd=0.8,fill=NA,aes(color=percent))+theme_bw(base_size=14)+scale_color_gradient(name="Percent of\nMaximum",low="blue",high="red")

Multiple Threshold Approach

To systematically explore the variation in the relationship based onthreshold choice, let’s look at the slope of the relationship and itsconfidence interval over all thresholds. We will then plot thisrelationship.

germanyLinearSlopes<-getCoefTab(funcMaxed~ Diversity,data = germanyThresh,coefVar ="Diversity",family =quasipoisson(link="identity"))####### Plot the values of the diversity slope at# different levels of the threshold######germanSlopes<-ggplot(germanyLinearSlopes,aes(x=thresholds*100,y = estimate,ymax = estimate+1.96* std.error,ymin = estimate-1.96* std.error))+geom_ribbon(fill="grey50")+geom_point()+ylab("Change in Number of Functions per Addition of 1 Species\n")+xlab("\nThreshold (%)")+geom_abline(intercept=0,slope=0,lwd=1,linetype=2)+theme_bw(base_size=14)germanSlopes

We can easily see a number of indices mentioned in the text - Tmin,Tmax, Tmde, etc. Let’s pull them out with our indices function.

germanIDX<-getIndices(germanyLinearSlopes, germanyThresh, funcMaxed~ Diversity)germanIDX

#   Tmin Tmax Tmde Rmde.linear Pmde.linear     Mmin    Mmax    Mmde nFunc# 1 0.11 0.81 0.41   0.1734241    0.554957 5.050544 1.06222 4.53135     5

We can now add annotations to our earlier plot to make it more rich,and show us those crucial threshold values.

germanyLinearSlopes$estimate[which(germanyLinearSlopes$thresholds==germanIDX$Tmde)]

# [1] 0.1734241

germanyThresh$IDX<-0germanyThresh$IDX [which(germanyThresh$thresholds%in%c(germanIDX$Tmin, germanIDX$Tmax, germanIDX$Tmde))]<-1ggplot(data=germanyThresh,aes(x=Diversity,y=funcMaxed,group=percent))+ylab(expression("Number of Functions">= Threshold))+xlab("Species Richness")+geom_smooth(method="glm",method.args =list(family=quasipoisson(link="identity")),fill=NA,aes(color=percent,lwd=IDX))+theme_bw(base_size=14)+scale_color_gradient(name="Percent of\nMaximum",low="blue",high="red")+scale_size(range=c(0.3,5),guide="none")+annotate(geom="text",x=0,y=c(0.2,2,4.6),label=c("Tmax","Tmde","Tmin"))+annotate(geom="text",x=16.7,y=c(germanIDX$Mmin, germanIDX$Mmax, germanIDX$Mmde),label=c("Mmin","Mmax","Mmde"))

Or we can see where they are on the curve of slopes at differentthresholds.

germanSlopes+annotate(geom="text",y=c(-0.01,-0.01,-0.01, germanIDX$Rmde.linear+0.02),x=c(germanIDX$Tmin*100, germanIDX$Tmde*100, germanIDX$Tmax*100, germanIDX$Tmde*100),label=c("Tmin","Tmde","Tmax","Rmde"),color="black")

Comparison of Sites

Furthermore, we can extend this to the more of the BIODEPTH datasetto compare the strength of multifunctionality across sites. We begin byloading in the data and calculating slopes at different thresholds.We’ll subset down to Sheffield, Portugal, and Sweden for simplicity, butthis could be done for all sites quite easily.

##### Now we will look at the entire BIODEPTH dataset## Note the use of dplyr to generate the new data frame using all locations from biodepth.# If you're not using dplyr for your data aggregation, you should be!# It will save you a lot of time and headaches in the future.#####Read in data  to run sample analyses on the biodepth datadata(all_biodepth)sub_biodepth<-subset(all_biodepth, all_biodepth$location%in%c("Sheffield","Portugal","Sweden"))sub_biodepth$location<-factor(sub_biodepth$location,levels=c("Portugal","Sweden","Sheffield"))allVars<-qw(biomassY3, root3, N.g.m2,  light3, N.Soil, wood3, cotton3)varIdx<-which(names(sub_biodepth)%in% allVars)#re-normalize so that everything is on the same#sign-scale (e.g. the maximum level of a function is#the "best" function)sub_biodepth<- sub_biodepth%>%group_by(location)%>%mutate(    light3<-ifelse(sum(is.na(light3)==0),-1*light3+max(light3,na.rm=T),NA),    N.Soil<-ifelse(sum(is.na(N.Soil)==0),-1*N.Soil+max(N.Soil,na.rm=T),NA)  )%>%ungroup()#get thresholdsbdThreshes<-sub_biodepth%>%group_by(location)%>%nest()%>%mutate(fmaxed =map(data, getFuncsMaxed,vars=allVars,maxN=8))%>%ungroup()%>%unnest(fmaxed)####look at slopes#note, maxIT argument is due to some fits needing more iterations to convergebdLinearSlopes<-getCoefTab(funcMaxed~ Diversity,data=bdThreshes,groupVar=c("location","thresholds"),coefVar="Diversity",family=quasipoisson(link="identity"),control=list(maxit=800))

# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values# Error : no valid set of coefficients has been found: please supply starting values

Great, now we can look at the indices of these three new sites. We’lladd Germany in so we can compare them against the German values.

indexTable<-lapply(levels(bdLinearSlopes$location),function(x){  slopedata<-subset(bdLinearSlopes, bdLinearSlopes$location==x)  threshdata<-subset(bdThreshes, bdThreshes$location==x)  ret<-getIndices(slopedata, threshdata, funcMaxed~ Diversity)  ret<-cbind(location=x, ret)  ret})indexTable<-do.call(rbind, indexTable)indexTable<-rbind(data.frame(location="Germany", germanIDX), indexTable)indexTable

#     location Tmin Tmax Tmde Rmde.linear Pmde.linear     Mmin     Mmax     Mmde# 1    Germany 0.11 0.81 0.41   0.1734241   0.5549570 5.050544 1.062220 4.531350# 13  Portugal   NA 0.29 0.19   0.1715245   0.4002239       NA 4.466414 5.705561# 11    Sweden 0.05 0.47 0.11   0.1486650   0.2973299 5.943906 1.998405 5.725218# 12 Sheffield 0.06   NA 0.55   0.1198584   0.3595751 3.709067       NA 3.212588#    nFunc# 1      5# 13     6# 11     6# 12     4

Finally, let’s compare a visualization of the relevant curves atdifferent thresholds to the relationship between threshold and slope ofcurve to give the table of numbers some visual meaning.

library(gridExtra)bdCurvesAll<-qplot(Diversity, funcMaxed,data=bdThreshes,group=thresholds,alpha=I(0))+facet_wrap(~location)+#, scales="free") +scale_color_gradient(low="blue",high="red",name="Proportion of\nMaximum",guide=FALSE)+stat_smooth(method="glm",lwd=0.8,fill=NA,method.args =list(family=gaussian(link="identity"),control=list(maxit=200)),aes(color=thresholds))+ylab("\nNumber of Functions ≥ Threshold\n\n")+xlab("Species Richness")+theme_bw(base_size=15)+geom_text(data=data.frame(location =levels(bdThreshes$location),lab =paste(letters[1:3],")",sep=""),thresholds=c(NA,NA,NA)),x=1,y=7,mapping=aes(label=lab))#Plot it!slopePlot<-ggplot(bdLinearSlopes,aes(x=thresholds*100,y=estimate,ymin = estimate-2*std.error,ymax = estimate+2*std.error))+geom_ribbon(fill="grey50")+geom_point()+ylab("Change in Number of Functions\nper Addition of 1 Species\n")+xlab("Threshold (%)")+facet_wrap(~location)+#, scale="free") +geom_abline(intercept=0,slope=0,lwd=0.6,linetype=2)+theme_bw(base_size=15)+geom_text(data=data.frame(location =levels(bdThreshes$location),lab =paste(letters[4:6],")",sep=""),thresholds=c(NA,NA,NA),std.error=c(NA,NA,NA),estimate =c(NA,NA,NA)                            ),x=0.05,y=0.27,mapping=aes(label=lab))###Plot them in a single figuregrid.arrange(bdCurvesAll,             slopePlot)