# 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 Pither (University of British Columbia, Okanagan campus) # (jason.pither@ubc.ca) # --------------------------------------------------------- # THIS SCRIPT CONDUCTS ANALYSES ASSOCIATED WITH THE # ANALYSIS OF MAXIMUM QUADRAT RICHNESS VERSUS BIOMASS # OUTPUT: FIGURE S4 # --------------------------------------------------------- # --------------------------------------------------------- # REQUIRED FILES: # fraser_plotdata.csv # --------------------------------------------------------- # --------------------------------------------------------- # REQUIRED PACKAGES: library(quantreg) # --------------------------------------------------------- # --------------------------------------------------------- # 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 # remaining columns are pres/abs data alldata.prab <- alldata[,(dim(alldata.summary)[[2]]+1):length(alldata)] # convert NA values that should be zeros: alldata.prab[is.na(alldata.prab)] <- 0 # 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 # --------------------------------------------------------- numgrids <- length(unique(alldata.litter$grid)) # number of unique grids gridnumbers <- as.character(unique(alldata.litter$grid)) # vector of grid identifiers numspec <- length(alldata.prab) # number of unique species across all grids numgridcells <- 64 # number in 8 x 8 grid # create species richness array, dimensions 8 x 8 x number of complete grids (151) # this literally stacks each grid on one-another, and places richness values in cells; sr.1m.array <- array(dim = c(sqrt(numgridcells),sqrt(numgridcells),numgrids),0) # give dimnames as described in document, and with appropriate grid number identifiers dimnames(sr.1m.array) <- list(1:8,c("A","B","C","D","E","F","G","H"),as.character(gridnumbers)); for (i in 1:numgrids){ sr.1m.array[,,gridnumbers[i]] <- matrix(nrow = sqrt(numgridcells),ncol = sqrt(numgridcells),alldata.litter[alldata.litter$grid == gridnumbers[i],"sr"]) print(i); } # DO SAME FOR TOTAL BIOMASS # create biomass array, dimensions 8 x 8 x number of complete grids (141) tot.biomass.array <- array(dim = c(sqrt(numgridcells),sqrt(numgridcells),numgrids),0); dimnames(tot.biomass.array) <- list(1:8,c("A","B","C","D","E","F","G","H"),as.character(gridnumbers)); for (i in 1:numgrids){ tot.biomass.array[,,gridnumbers[i]] <- matrix(nrow = sqrt(numgridcells),ncol = sqrt(numgridcells),alldata.litter[alldata.litter$grid == gridnumbers[i],"tot.bio"]) print(i); } # ------------------- # Find the maximum richness in a grid, then match it with its biomass: # ------------------- max.rich.matching.totbio <- data.frame(max.rich = rep(0,dim(sr.1m.array)[[3]]),totbio = rep(0,dim(sr.1m.array)[[3]])) for (i in 1:dim(sr.1m.array)[[3]]){ max.rich.matching.totbio[i,] <- c( sr.1m.array[,,i][order(c(sr.1m.array[,,i]),decreasing = T)[1]], tot.biomass.array[,,i][order(c(sr.1m.array[,,i]),decreasing = T)[1]] ) } # so for all analyses below, N = 151; max.rich.matching.totbio$log10.totbio <- log10(max.rich.matching.totbio$totbio+1); max.rich.matching.totbio$log10.max.rich <- log10(max.rich.matching.totbio$max.rich); # ------------------- # Least squares regression of maximum quadrat richness vs. associated biomass: # ------------------- log.maxrich.vs.matching.biomass.lm <- lm(log10.max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio) summary(log.maxrich.vs.matching.biomass.lm) # how does eliminating the two bottom right points influence the regression? summary(lm(log10.max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio[-c(86,88),])) # *** It cuts the adjusted R-square in half, but the form of relationship remains the same: confint(lm(log10.max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio[-c(86,88),])) # 2.5 % 97.5 % #(Intercept) -4.4007424 -0.6284784 #log10.totbio 1.7042512 4.7183167 #I(log10.totbio^2) -0.9721211 -0.3756569 confint(lm(log10.max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio)) # 2.5 % 97.5 % #(Intercept) -5.781107 -2.3147721 #log10.totbio 3.129082 5.8618797 #I(log10.totbio^2) -1.205542 -0.6716587 # ------------------- # Quantile regression of maximum quadrat richness vs. associated biomass: # ------------------- temp.quad <- rq(max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio,tau = c(0.8,0.9,0.95)); # summary(temp.quad,se = "boot",R = 1000) # now just get 95 quantile for plotting temp.quad.95 <- rq(max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio,tau = 0.95); # calculate approximate "r-square" temp.quad.95.1 <- rq(max.rich~1,max.rich.matching.totbio,tau = 0.95); rho <- function(u,tau = 0.5)u*(tau-(u<0)); R1 <- 1 - temp.quad.95$rho/temp.quad.95.1$rho R1 #[1] 0.1359595 # and eliminating the lower right points? # first check adjusted R-square: temp.quad.95.1.no.out <- rq(max.rich~1,max.rich.matching.totbio[-c(86,88),],tau = 0.95) temp.quad.95.no.out <- rq(max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio[-c(86,88),],tau = 0.95) 1 - temp.quad.95.no.out$rho/temp.quad.95.1.no.out$rho # [1] 0.128473 ** minimal influence on pseudo Rsquare temp.quad.no.out <- rq(max.rich~log10.totbio+I(log10.totbio^2),max.rich.matching.totbio[-c(86,88),],tau = c(0.8,0.9,0.95)); summary(temp.quad.no.out,se = "boot",R = 1000) # And eliminating the right-hand points does not affect form # ------------------- # Begin Figure S4 # ------------------- # set default par par.default <- par(no.readonly = TRUE) par(par.default) par(mfrow = c(1,1),mar = c(5,5,3,3)) plot(log10.max.rich~log10.totbio,max.rich.matching.totbio,las = 1,cex.lab = 1.4,cex.axis = 1.2,xlab = expression("Total biomass associated with maximum richness" ~ log[10]~(g~m^{-2})),cex.axis = 1,cex.lab = 1.1,ylab = expression("Maximum grid richness " ~ (log[10])),cex = 1.2,pch = 1,col = "darkgrey",lwd = 1.3); lines(max.rich.matching.totbio$log10.totbio[order(max.rich.matching.totbio$log10.totbio)],predict(log.maxrich.vs.matching.biomass.lm,newdata = max.rich.matching.totbio[order(max.rich.matching.totbio$log10.totbio),],type = "response"),col = "lightgrey",lwd = 2) # add quantile 95 regression: lines(max.rich.matching.totbio$log10.totbio[order(max.rich.matching.totbio$log10.totbio)],log10(predict(temp.quad.95,newdata = max.rich.matching.totbio[order(max.rich.matching.totbio$log10.totbio),],type = "response")),lwd = 2); # ------------------- # END Figure S4 # -------------------