# This is an "R" script documenting work associated with the following # publication: # Fraser, L.H., Pither, J., Jentsch, A., Sternberg, M., Zobel, M., and the # HerbDivNet # global network. 2015. Worldwide evidence of a unimodal # relationship between productivity and plant species richness. Science. # P.I.: Lauchlan Fraser (Thompson Rivers University) (lfraser@tru.ca) # Script author: Jason Fridley (Syracuse University) # (fridley@syr.edu) # --------------------------------------------------------- # THIS SCRIPT CONDUCTS ANALYSES ASSOCIATED WITH THE # HIERARCHICAL BAYESIAN ANALYSIS OF 95TH PERCENTILE # OUTPUT: FIGURE S3 # --------------------------------------------------------- # --------------------------------------------------------- # REQUIRED FILES: # fraser_plotdata.csv # --------------------------------------------------------- # --------------------------------------------------------- # REQUIRED PACKAGES: library(R2jags) # --------------------------------------------------------- # --------------------------------------------------------- # read in data, and do initial data cleaning # --------------------------------------------------------- alldata <- read.csv("fraser_plotdata.csv", header = T, strip.white = T) # first 10 columns are summary data alldata.summary <- alldata[,c("grid","original.name","community.type","country","pi","plot","biomass","litter","tot.bio","sr")] # column headings are self-explanatory; but note that "pi" is what we equate with "site" throughout the paper # exclude grids that had not litter measurements conducted in all 64 quadrats complete.grid.numbers <- setdiff(unique(alldata.summary$grid),names(table(alldata.summary[is.na(alldata.summary$litter),"grid"])[table(alldata.summary[is.na(alldata.summary$litter),"grid"]) == 64])) num.complete.grids <- length(complete.grid.numbers) # now we have 151 grids with complete or mostly complete measurements # (10 have between 1 and 4 quadrats missing litter) # now subset data to only those grids alldata.litter <- alldata.summary[!is.na(match(alldata.summary$grid,as.numeric(complete.grid.numbers))),] # there are some quadrats with zero species, and two (7023, 7030) # in grid 110 that have valid litter biomass (zero) but no live species; # these remain as VALID data points # put "NA" in all other quadrats na.rownames <- row.names(alldata.litter[alldata.litter$sr == 0,][setdiff(row.names(alldata.litter[alldata.litter$sr == 0,]),c("7023","7030")),]) # remaining quadrats with neither biomass nor species can be "NA" alldata.litter[na.rownames,"biomass"] <- NA alldata.litter[na.rownames,"litter"] <- NA alldata.litter[na.rownames,"tot.bio"] <- NA alldata.litter[na.rownames,"sr"] <- NA # create a log10 total biomass (+1) column alldata.litter$log10.tot.bio <- log10(alldata.litter$tot.bio+1) # --------------------------------------------------------- # DONE data input / cleaning # --------------------------------------------------------- # --------------------------------------------------------- #Bayesian models for total biomass (live + litter), Fridley, 3-17-15 # --------------------------------------------------------- y = alldata.litter$sr x = alldata.litter$log10.tot.bio groups = as.numeric(alldata.litter$pi) N.groups = max(groups) grid = alldata.litter$grid N.grids = max(grid) plot(x,y,col=pi) #8 missing values in x, 2 in y; omit dat = na.exclude(data.frame(x,y,groups,grid)) y = dat$y; x = dat$x; groups = dat$groups; grids = dat$grid #limit modeling to where data are abundant, between 50 and 1500 g/m2 #can skip this step, but may not converge true = x>log10(50)&x