R/vector_plotting.R
In ctkremer/priceTools: Tools for analyzing biodiversity-ecosystem function relationships
Defines functions
process.data.cafe
process.data.bef
process.data.price
leap.zig
leap.zig.both
leap.zig.bef
leap.zig.cafe
leap.zig.price
Documented in
leap.zig
leap.zig.bef
leap.zig.both
leap.zig.cafe
leap.zig.price
process.data.bef
process.data.cafe
process.data.price
######################## Vector Plotting ########################
#' Process data from a single treatment prior to making a vector plot
#'
#' These functions take the data generated by pairwise.price() and processes it into the
#' terms needed in CAFE, BEF, or Price component vector plots.
#'
#' @param data Pairwise Price data
#' @param group.vars A vector of grouping variables, if any
#' @param standardize Should ecosystem function values be standardized against baseline? T/F
#'
#' @return A list of three data sets with different levels of aggregation used in subsequent vector plots. Although still exported, this function has become essentially an internal function, called by members of the leap.zig family of functions, and may rarely be useful to call directly.
#'
#' @examples
#'
#' # write one
#' set.seed(36)
#'
#' # Data frame containing multiple communities we want to compare
#' cms<-data.frame(comm.id=sort(rep(seq(1,3),6)),
#' species=rep(LETTERS[seq(1,6)],3),
#' func=rpois(6*3,lambda = 2))
#'
#' #Identify one (or more) grouping columns
#' cms<-group_by(cms,comm.id)
#'
#' # Perform pairwise comparisons of all communities in cms identified by comm.id
#' pp<-pairwise.price(cms,species='species',func='func')
#' pp<-group.columns(pp,gps=c('comm.id'))
#'
#' process.data.cafe(data=pp,group.vars='comm.id')
#' process.data.bef(data=pp,group.vars='comm.id')
#' process.data.price(data=pp,group.vars='comm.id')
#'
#' @export
#' @import dplyr
process.data.cafe<-function(data, group.vars=NULL, standardize=TRUE){
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE","SL","SG","SR","CE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
data$SG <- data$SL + data$SG # net SG change:
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SL <- data$x.func + data$SL
data$SG <- data$SL + data$SG
}
# CAFE component richness function
# base = c(x.rich,x.func)
# SL = c(c.rich,SL)
# SG = c(y.rich,SL+SG)
# CDE = c(y.rich,y.func)
cols <- c(group.vars,'x.func','SL','SG','y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols],id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add richness column:
p2$rich <- ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SL", p2$c.rich,
ifelse(p2$variable == "SG", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA))))
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable,levels=c("x.func","SL","SG","y.func"),
labels=c("baseline","SL","SG","comparison"))
p3b$variable <- factor(p3b$variable,levels=c("x.func","SL","SG","y.func"),
labels=c("baseline","SL","SG","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SG",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SL","SG"),
labels=c("SL vector","SG vector","CDE vector"))
p2$variable <- factor(p2$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p3b$variable <- factor(p3b$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p4$variable <- factor(p4$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' @describeIn process.data.cafe Process data for plotting BEF components.
#' @export
#' @import dplyr
process.data.bef<-function(data, group.vars=NULL, standardize=TRUE){
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE","SL","SG","SR","CE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SR <- data$x.func + data$SR
}
# BEF component
# base = c(x.rich,x.func)
# SR = c(x.rich,x.func + SRE.L + SRE.G)
# CE = c(y.rich,y.func) = c(y.rich, x.func + SRE.L + SRE.G + SIEs + CDE)
cols <- c(group.vars,'x.func','SR','y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols], id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add composite richness column:
p2$rich <- ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SR", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA)))
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>%
summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
# stagger rows:
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable,levels=c("x.func","SR","y.func"),
labels=c("baseline","SR","comparison"))
p3b$variable <- factor(p3b$variable,levels=c("x.func","SR","y.func"),
labels=c("baseline","SR","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SR",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SR"),
labels=c("SR vector","CE vector"))
p2$variable <- factor(p2$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p3b$variable <- factor(p3b$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p4$variable <- factor(p4$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' @describeIn process.data.cafe Process data for plotting 5-part Price components.
#' @export
#' @import dplyr
process.data.price<-function(data,group.vars=NULL,standardize=TRUE){
### standardize
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
data$SIE.L <- data$SRE.L + data$SIE.L # create net change vectors
data$SRE.G <- data$SIE.L + data$SRE.G
data$SIE.G <- data$SRE.G + data$SIE.G
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SRE.L <- data$x.func + data$SRE.L
data$SIE.L <- data$SRE.L + data$SIE.L
data$SRE.G <- data$SIE.L + data$SRE.G
data$SIE.G <- data$SRE.G + data$SIE.G
}
# base = c(x.rich,x.func)
# SRE.L = c(c.rich,0 + SRE.L)
# SIE.L = c(c.rich,0 + SRE.L + SIE.L = SL)
# SRE.G = c(y.rich,0 + SRE.L + SIE.L + SRE.G = SL + SRE.G)
# SIE.G = c(y.rich,0 + SRE.L + SIE.L + SRE.G + SIE.G = SL + SG)
# CDE = c(y.rich,y.func)
cols <- c(group.vars,'x.func','SRE.L','SIE.L','SRE.G','SIE.G',
'y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols], id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add richness column:
p2$rich<-ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SRE.L", p2$c.rich,
ifelse(p2$variable == "SIE.L", p2$c.rich,
ifelse(p2$variable == "SRE.G", p2$y.rich,
ifelse(p2$variable == "SIE.G", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA)
)
)
)
)
)
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
# stagger rows:
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable, levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G",
"y.func"),
labels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G","comparison"))
p3b$variable <- factor(p3b$variable, levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G",
"y.func"),
labels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SIE.G",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G"),
labels=c("SRE.L vector","SIE.L vector","SRE.G vector","SIE.G vector","CDE vector"))
p2$variable <- factor(p2$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p3b$variable <- factor(p3b$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p4$variable <- factor(p4$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' Plot changes in ecosystem function between communities as vector diagram.
#'
#' This wrapper function takes the data generated by \code{\link{pairwise.price()}} and produces vector plots using either BEF, CAFE, or 5-part Price components.
#'
#' @param data Pairwise Price data
#' @param type Specify the type of vector diagram to draw ("cafe","bef","both", or "price")
#' @param group.vars A vector of grouping variables, if any
#' @param standardize Should ecosystem function values be standardized against baseline
#' @param ... Additional graphical options passed to lower-level functions. See \code{\link{leap.zig.bef}}, \code{\link{leap.zig.cafe}}, \code{\link{leap.zig.both}}, \code{\link{leap.zig.price}}
#'
#' @return A ggplot object.
#'
#' @examples
#'
#' # Data setup
#' data(cedarcreek)
#' head(cedarcreek)
#'
#' #Identify one grouping columns
#' cc2<-group_by(cedarcreek,NTrt,NAdd,Plot)
#'
#' # Perform pairwise comparisons of all communities in cms identified by comm.id
#' # (takes ~30 sec)
#' pp<-pairwise.price(cc2,species='Species',func='Biomass')
#'
#' # Organize/format the results, and pull out a subset using NAdd.x=0 as the control/baseline site
#' pp<-group.columns(pp,gps=c('NTrt','NAdd'))
#' pp<-pp[pp$NAdd.x=='0',]
#' dat1<-pp[pp$NAdd %in% c('0 27.2'),]
#' dat1.ctrl<-pp[pp$NAdd %in% c('0 0'),]
#'
#' # Demonstrate vector plotting for comparisons between the control level of N addition (0) and the highest level of N addition (27.2):
#' leap.zig(dat1,type='cafe')
#'
#' # Or sets of vectors grouped by N addition level:
#' leap.zig(pp,type='cafe',group.vars=c('NAdd.y'),raw.points=F,ylim=c(-100,700))
#'
#' # The plots returned are ggplot objects, so additional design specifications can be added on:
#' leap.zig(dat1,type='cafe')+
#' ggtitle('Enrichment \n(0 vs. 27.2)')
#'
#' # Control plot window
#' leap.zig(dat1,type='cafe',xlim=c(3,18),ylim=c(-100,700))
#'
#' # Turn on errorbars associated with vector endpoints
#' leap.zig(dat1,type='cafe',xlim=c(3,18),ylim=c(-100,700),error.bars=T)
#'
#' # Turn on/off standardization of vector magnitude (taken as % change relative to baseline)
#' leap.zig(dat1,type='cafe',standardize=F)
#'
#' # Make a plot without showing individual points:
#' leap.zig(dat1,type='cafe',standardize=F,raw.points=F)
#'
#' # Use other styles of vector arrangements:
#' leap.zig(dat1,type='bef',standardize=F,raw.points=F)
#' leap.zig(dat1,type='price',standardize=F,raw.points=F)
#'
#' # BUSTED!
#' # leap.zig(dat1,type='both')
#'
#' # Turn legend off
#' leap.zig(dat1,type='price',standardize=F,raw.points=F,legend=F)
#'
#' # Combine multiple plots in separate panels
#' # (note: faceting currently doesn't work with leap.zig)
#' s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=T,vectors=T,raw.points = F,legend=F)
#' s2 <- leap.zig(dat1,type='bef', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=T,vectors=T,raw.points = F,legend=F)
#'
#' library(gridExtra)
#' grid.arrange(s1,s2,nrow=1)
#'
#' # Or on top of each other
#' s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=F,vectors=T,raw.points = F,legend=F)
#' leap.zig(dat1.ctrl,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=F,vectors=T,raw.points = F,legend=F,
#' add=T,old.plot=s1)
#'
#' # BUSTED - due to conflicting color scales
#' # s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' # error.bars=T,vectors=T,raw.points = F,legend=F)
#' # leap.zig(dat1,type='bef', xlim=c(3,18),ylim=c(-100,700),
#' # error.bars=T,vectors=T,raw.points = F,legend=F,
#' # add=T,old.plot=s1)
#'
#' @export
leap.zig<-function(data, type="cafe", group.vars=NULL, standardize=TRUE, ...){
switch(type,
cafe={
tmp <- process.data.cafe(data, group.vars, standardize)
leap.zig.cafe(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
bef={
tmp <- process.data.bef(data, group.vars, standardize)
leap.zig.bef(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
both={
tmp <- process.data.cafe(data, group.vars, standardize)
tmp.bef <- process.data.bef(data, group.vars, standardize)
tmp[[1]] <- unique(rbind(tmp[[1]], tmp.bef[[1]]))
tmp[[2]] <- unique(rbind(tmp[[2]], tmp.bef[[2]]))
tmp[[3]] <- unique(rbind(tmp[[3]], tmp.bef[[3]]))
lvls <- c("baseline","SL","SG","SR","comparison","SL vector","SG vector","CDE vector",
"SR vector","CE vector")
tmp[[1]]$variable <- factor(tmp[[1]]$variable, levels=lvls)
leap.zig.both(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
price={
tmp <- process.data.price(data, group.vars, standardize)
leap.zig.price(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
"Error! Invalid plot method in leap.zig()"
)
}
#' Internal functions for generating vector plots of ecosystem function components.
#'
#' This suite of functions is accessed by \code{\link{leap.zig()}} and produces vector plots
#' using BEF, CAFE, and Price components. \code{leap.zig.both} plots both the BEF and CAFE vectors.
#'
#' @param tmp Data to plot
#' @param xlim Plot's x limits
#' @param ylim Plot's y limits
#' @param loc.standarize Are these standardized vectors
#' @param error.bars Plot error bars
#' @param raw.points Plot raw data points at level of community pairs
#' @param vectors Plot averaged vectors
#' @param legend Show legend
#' @param old.plot ggplot object from previous \code{leap.zig()} call
#' @param add Add new plot to object provided in old.plot option
#'
#' @return A ggplot object.
#'
#' @examples
#'
#' # Load data and run pairwise comparisons of communities
#' cc2<-group_by(cedarcreek,NTrt,NAdd,Plot)
#' pp<-pairwise.price(cc2,species='Species',func='Biomass')
#'
#' # Organize/format the results, and pull out a subset using NTrt=1 as the control site
#' pp<-group.columns(pp,gps=c('NTrt','NAdd'))
#' pp<-pp[pp$NTrt.x==1,]
#' dat1<-pp[pp$NAdd %in% c('0 27.2'),]
#'
#' # Process and plot the result
#' tmp <- process.data.bef(dat1, group.vars=NULL, standardize=F)
#' leap.zig.bef(tmp, loc.standardize=F, group.vars=NULL,raw.points = F,legend = F)
#'
#' @export
#' @import ggplot2
leap.zig.both<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE,group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp[[3]]$variable <- factor(as.character(tmp[[3]]$variable),
levels=c("SL vector","SG vector","CDE vector","SR vector","CE vector"))
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot() + geom_hline(yintercept=0, linetype=2)
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
colour='gray') +
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x, y=mean.y),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp[[3]], aes(colour=variable, x=mean.x, y=mean.y),
arrow=arrow(length=unit(0.2,"cm"), ends="first"))
}
cols <- c(alpha('#d95f02'),alpha('#1b9e77'),alpha('#7570b3'),alpha('#7fcdbb'),alpha('#e34a33'))
trcols <- c(alpha('#d95f02',0.1),alpha('#1b9e77',0.1),alpha('#7fcdbb',0.1),alpha('#7570b3',0.1))
if(raw.points==FALSE & vectors==TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA,NA,NA,NA),
linetype=c(1,1,1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend( override.aes=list(
shape=c(19,1,1,1,1,NA,NA,NA,NA,NA),
linetype=c(0,0,0,0,0,1,1,1,1,1),
colour=c('black',cols[1:4],cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim) == 2){
xlim.outer <- c(min(xlim.D[1], xlim[1]), max(xlim.D[2], xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1], ylim[1]), max(ylim.D[2], ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + theme_bw() + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for BEF components
#' @export
#' @import ggplot2
leap.zig.bef<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable), levels=c("SR vector","CE vector"))
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.bef()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#008002'),alpha('#56F256'))
trcols <- c(alpha('#008002',0.1),alpha('#56F256',0.1))
if(raw.points == FALSE & vectors == TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA),linetype=c(1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend( override.aes=list(shape=c(19,1,1,NA,NA),
linetype=c(0,0,0,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim)==2){
xlim.outer <- c(min(xlim.D[1],xlim[1]), max(xlim.D[2],xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1],ylim[1]), max(ylim.D[2],ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}
if(legend == FALSE){
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for CAFE components
#' @export
#' @import ggplot2
leap.zig.cafe<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable),
levels=c("SL vector","SG vector","CDE vector"))
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.cafe()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#FE0000'), alpha('#0025FF'), alpha('#A01AE5'))
trcols <- c(alpha('#FE0000',0.1), alpha('#0025FF',0.1), alpha('#A01AE5',0.1))
if(raw.points==FALSE & vectors==TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA,NA),
linetype=c(1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols))+
guides(colour=guide_legend( override.aes=list(shape=c(19,1,1,1,NA,NA,NA),
linetype=c(0,0,0,0,1,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim)==2){
xlim.outer <- c(min(xlim.D[1],xlim[1]), max(xlim.D[2],xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1],ylim[1]), max(ylim.D[2],ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for Price components
#' @export
#' @import ggplot2
leap.zig.price <- function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable),
levels=c("SRE.L vector","SIE.L vector","SRE.G vector",
"SIE.G vector","CDE vector"))
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.cafe()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#DB0102'), alpha('#FF6500'), alpha('#001ECC'),
alpha('#48C5E9'), alpha('#A01AE5'))
trcols <- c(alpha('#DB0102',0.1), alpha('#FF6500',0.1), alpha('#001ECC',0.1),
alpha('#48C5E9',0.1), alpha('#A01AE5',0.1))
if(raw.points == FALSE & vectors == TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend(override.aes=list(shape=c(NA,NA,NA,NA,NA),
linetype=c(1,1,1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend(
override.aes=list(shape=c(19,1,1,1,1,1,NA,NA,NA,NA,NA),
linetype=c(0,0,0,0,0,0,1,1,1,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim) == 2){
xlim.outer <- c(min(xlim.D[1], xlim[1]), max(xlim.D[2], xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1], ylim[1]), max(ylim.D[2], ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
ctkremer/priceTools documentation built on May 28, 2019, 7:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions process.data.cafe process.data.bef process.data.price leap.zig leap.zig.both leap.zig.bef leap.zig.cafe leap.zig.price
Documented in leap.zig leap.zig.bef leap.zig.both leap.zig.cafe leap.zig.price process.data.bef process.data.cafe process.data.price
######################## Vector Plotting ########################
#' Process data from a single treatment prior to making a vector plot
#'
#' These functions take the data generated by pairwise.price() and processes it into the
#' terms needed in CAFE, BEF, or Price component vector plots.
#'
#' @param data Pairwise Price data
#' @param group.vars A vector of grouping variables, if any
#' @param standardize Should ecosystem function values be standardized against baseline? T/F
#'
#' @return A list of three data sets with different levels of aggregation used in subsequent vector plots. Although still exported, this function has become essentially an internal function, called by members of the leap.zig family of functions, and may rarely be useful to call directly.
#'
#' @examples
#'
#' # write one
#' set.seed(36)
#'
#' # Data frame containing multiple communities we want to compare
#' cms<-data.frame(comm.id=sort(rep(seq(1,3),6)),
#' species=rep(LETTERS[seq(1,6)],3),
#' func=rpois(6*3,lambda = 2))
#'
#' #Identify one (or more) grouping columns
#' cms<-group_by(cms,comm.id)
#'
#' # Perform pairwise comparisons of all communities in cms identified by comm.id
#' pp<-pairwise.price(cms,species='species',func='func')
#' pp<-group.columns(pp,gps=c('comm.id'))
#'
#' process.data.cafe(data=pp,group.vars='comm.id')
#' process.data.bef(data=pp,group.vars='comm.id')
#' process.data.price(data=pp,group.vars='comm.id')
#'
#' @export
#' @import dplyr
process.data.cafe<-function(data, group.vars=NULL, standardize=TRUE){
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE","SL","SG","SR","CE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
data$SG <- data$SL + data$SG # net SG change:
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SL <- data$x.func + data$SL
data$SG <- data$SL + data$SG
}
# CAFE component richness function
# base = c(x.rich,x.func)
# SL = c(c.rich,SL)
# SG = c(y.rich,SL+SG)
# CDE = c(y.rich,y.func)
cols <- c(group.vars,'x.func','SL','SG','y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols],id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add richness column:
p2$rich <- ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SL", p2$c.rich,
ifelse(p2$variable == "SG", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA))))
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable,levels=c("x.func","SL","SG","y.func"),
labels=c("baseline","SL","SG","comparison"))
p3b$variable <- factor(p3b$variable,levels=c("x.func","SL","SG","y.func"),
labels=c("baseline","SL","SG","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SG",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SL","SG"),
labels=c("SL vector","SG vector","CDE vector"))
p2$variable <- factor(p2$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p3b$variable <- factor(p3b$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p4$variable <- factor(p4$variable, levels=c("baseline","SL","SG","comparison",
"SL vector","SG vector","CDE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' @describeIn process.data.cafe Process data for plotting BEF components.
#' @export
#' @import dplyr
process.data.bef<-function(data, group.vars=NULL, standardize=TRUE){
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE","SL","SG","SR","CE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SR <- data$x.func + data$SR
}
# BEF component
# base = c(x.rich,x.func)
# SR = c(x.rich,x.func + SRE.L + SRE.G)
# CE = c(y.rich,y.func) = c(y.rich, x.func + SRE.L + SRE.G + SIEs + CDE)
cols <- c(group.vars,'x.func','SR','y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols], id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add composite richness column:
p2$rich <- ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SR", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA)))
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>%
summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
# stagger rows:
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable,levels=c("x.func","SR","y.func"),
labels=c("baseline","SR","comparison"))
p3b$variable <- factor(p3b$variable,levels=c("x.func","SR","y.func"),
labels=c("baseline","SR","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SR",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SR"),
labels=c("SR vector","CE vector"))
p2$variable <- factor(p2$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p3b$variable <- factor(p3b$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p4$variable <- factor(p4$variable,levels=c("baseline","SR","comparison","SR vector",
"CE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' @describeIn process.data.cafe Process data for plotting 5-part Price components.
#' @export
#' @import dplyr
process.data.price<-function(data,group.vars=NULL,standardize=TRUE){
### standardize
if(standardize == TRUE){
comps <- c("SRE.L","SRE.G","SIE.L","SIE.G","CDE")
data[,comps] <- 100*data[,comps]/data$x.func # X function scaled
data$y.func <- 100*(data$y.func - data$x.func)/data$x.func # Y function scaled
data$x.func <- 0 # X function set as baseline.
data$SIE.L <- data$SRE.L + data$SIE.L # create net change vectors
data$SRE.G <- data$SIE.L + data$SRE.G
data$SIE.G <- data$SRE.G + data$SIE.G
} else{
data$x.func <- data$x.func
data$y.func <- data$y.func
data$SRE.L <- data$x.func + data$SRE.L
data$SIE.L <- data$SRE.L + data$SIE.L
data$SRE.G <- data$SIE.L + data$SRE.G
data$SIE.G <- data$SRE.G + data$SIE.G
}
# base = c(x.rich,x.func)
# SRE.L = c(c.rich,0 + SRE.L)
# SIE.L = c(c.rich,0 + SRE.L + SIE.L = SL)
# SRE.G = c(y.rich,0 + SRE.L + SIE.L + SRE.G = SL + SRE.G)
# SIE.G = c(y.rich,0 + SRE.L + SIE.L + SRE.G + SIE.G = SL + SG)
# CDE = c(y.rich,y.func)
cols <- c(group.vars,'x.func','SRE.L','SIE.L','SRE.G','SIE.G',
'y.func','x.rich','c.rich','y.rich')
p2 <- reshape2::melt(data[,cols], id.vars=c(group.vars,'x.rich','c.rich','y.rich'))
# add richness column:
p2$rich<-ifelse(p2$variable == "x.func", p2$x.rich,
ifelse(p2$variable == "SRE.L", p2$c.rich,
ifelse(p2$variable == "SIE.L", p2$c.rich,
ifelse(p2$variable == "SRE.G", p2$y.rich,
ifelse(p2$variable == "SIE.G", p2$y.rich,
ifelse(p2$variable == "y.func", p2$y.rich, NA)
)
)
)
)
)
# summarize raw points:
if(!is.null(group.vars)){
p3b <- p2 %>% group_by_(.dots=c(group.vars,'variable')) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}else{
p3b <- p2 %>% group_by(variable) %>% summarise(mean.y=mean(value),
y.qt.lw=quantile(value, probs=0.025),
y.qt.up=quantile(value, probs=0.975),
mean.x=mean(rich),
x.qt.lw=quantile(rich, probs=0.025),
x.qt.up=quantile(rich, probs=0.975))
}
# stagger rows:
p4 <- p3b
# Organize factor levels for plotting:
p2$variable <- factor(p2$variable, levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G",
"y.func"),
labels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G","comparison"))
p3b$variable <- factor(p3b$variable, levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G",
"y.func"),
labels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G","comparison"))
# updated code:
p4$variable <- as.character(p4$variable)
p4$variable <- ifelse(p4$variable=="y.func","SIE.G",p4$variable)
p4$variable <- factor(p4$variable,levels=c("x.func","SRE.L","SIE.L","SRE.G","SIE.G"),
labels=c("SRE.L vector","SIE.L vector","SRE.G vector","SIE.G vector","CDE vector"))
p2$variable <- factor(p2$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p3b$variable <- factor(p3b$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p4$variable <- factor(p4$variable, levels=c("baseline","SRE.L","SIE.L","SRE.G","SIE.G",
"comparison","SRE.L vector","SIE.L vector",
"SRE.G vector","SIE.G vector","CDE vector"))
p3b <- p3b[p3b$variable != "baseline",]
return(list(p2, p3b, p4))
}
#' Plot changes in ecosystem function between communities as vector diagram.
#'
#' This wrapper function takes the data generated by \code{\link{pairwise.price()}} and produces vector plots using either BEF, CAFE, or 5-part Price components.
#'
#' @param data Pairwise Price data
#' @param type Specify the type of vector diagram to draw ("cafe","bef","both", or "price")
#' @param group.vars A vector of grouping variables, if any
#' @param standardize Should ecosystem function values be standardized against baseline
#' @param ... Additional graphical options passed to lower-level functions. See \code{\link{leap.zig.bef}}, \code{\link{leap.zig.cafe}}, \code{\link{leap.zig.both}}, \code{\link{leap.zig.price}}
#'
#' @return A ggplot object.
#'
#' @examples
#'
#' # Data setup
#' data(cedarcreek)
#' head(cedarcreek)
#'
#' #Identify one grouping columns
#' cc2<-group_by(cedarcreek,NTrt,NAdd,Plot)
#'
#' # Perform pairwise comparisons of all communities in cms identified by comm.id
#' # (takes ~30 sec)
#' pp<-pairwise.price(cc2,species='Species',func='Biomass')
#'
#' # Organize/format the results, and pull out a subset using NAdd.x=0 as the control/baseline site
#' pp<-group.columns(pp,gps=c('NTrt','NAdd'))
#' pp<-pp[pp$NAdd.x=='0',]
#' dat1<-pp[pp$NAdd %in% c('0 27.2'),]
#' dat1.ctrl<-pp[pp$NAdd %in% c('0 0'),]
#'
#' # Demonstrate vector plotting for comparisons between the control level of N addition (0) and the highest level of N addition (27.2):
#' leap.zig(dat1,type='cafe')
#'
#' # Or sets of vectors grouped by N addition level:
#' leap.zig(pp,type='cafe',group.vars=c('NAdd.y'),raw.points=F,ylim=c(-100,700))
#'
#' # The plots returned are ggplot objects, so additional design specifications can be added on:
#' leap.zig(dat1,type='cafe')+
#' ggtitle('Enrichment \n(0 vs. 27.2)')
#'
#' # Control plot window
#' leap.zig(dat1,type='cafe',xlim=c(3,18),ylim=c(-100,700))
#'
#' # Turn on errorbars associated with vector endpoints
#' leap.zig(dat1,type='cafe',xlim=c(3,18),ylim=c(-100,700),error.bars=T)
#'
#' # Turn on/off standardization of vector magnitude (taken as % change relative to baseline)
#' leap.zig(dat1,type='cafe',standardize=F)
#'
#' # Make a plot without showing individual points:
#' leap.zig(dat1,type='cafe',standardize=F,raw.points=F)
#'
#' # Use other styles of vector arrangements:
#' leap.zig(dat1,type='bef',standardize=F,raw.points=F)
#' leap.zig(dat1,type='price',standardize=F,raw.points=F)
#'
#' # BUSTED!
#' # leap.zig(dat1,type='both')
#'
#' # Turn legend off
#' leap.zig(dat1,type='price',standardize=F,raw.points=F,legend=F)
#'
#' # Combine multiple plots in separate panels
#' # (note: faceting currently doesn't work with leap.zig)
#' s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=T,vectors=T,raw.points = F,legend=F)
#' s2 <- leap.zig(dat1,type='bef', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=T,vectors=T,raw.points = F,legend=F)
#'
#' library(gridExtra)
#' grid.arrange(s1,s2,nrow=1)
#'
#' # Or on top of each other
#' s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=F,vectors=T,raw.points = F,legend=F)
#' leap.zig(dat1.ctrl,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' error.bars=F,vectors=T,raw.points = F,legend=F,
#' add=T,old.plot=s1)
#'
#' # BUSTED - due to conflicting color scales
#' # s1 <- leap.zig(dat1,type='cafe', xlim=c(3,18),ylim=c(-100,700),
#' # error.bars=T,vectors=T,raw.points = F,legend=F)
#' # leap.zig(dat1,type='bef', xlim=c(3,18),ylim=c(-100,700),
#' # error.bars=T,vectors=T,raw.points = F,legend=F,
#' # add=T,old.plot=s1)
#'
#' @export
leap.zig<-function(data, type="cafe", group.vars=NULL, standardize=TRUE, ...){
switch(type,
cafe={
tmp <- process.data.cafe(data, group.vars, standardize)
leap.zig.cafe(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
bef={
tmp <- process.data.bef(data, group.vars, standardize)
leap.zig.bef(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
both={
tmp <- process.data.cafe(data, group.vars, standardize)
tmp.bef <- process.data.bef(data, group.vars, standardize)
tmp[[1]] <- unique(rbind(tmp[[1]], tmp.bef[[1]]))
tmp[[2]] <- unique(rbind(tmp[[2]], tmp.bef[[2]]))
tmp[[3]] <- unique(rbind(tmp[[3]], tmp.bef[[3]]))
lvls <- c("baseline","SL","SG","SR","comparison","SL vector","SG vector","CDE vector",
"SR vector","CE vector")
tmp[[1]]$variable <- factor(tmp[[1]]$variable, levels=lvls)
leap.zig.both(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
price={
tmp <- process.data.price(data, group.vars, standardize)
leap.zig.price(tmp, loc.standardize=standardize, group.vars=group.vars, ...)
},
"Error! Invalid plot method in leap.zig()"
)
}
#' Internal functions for generating vector plots of ecosystem function components.
#'
#' This suite of functions is accessed by \code{\link{leap.zig()}} and produces vector plots
#' using BEF, CAFE, and Price components. \code{leap.zig.both} plots both the BEF and CAFE vectors.
#'
#' @param tmp Data to plot
#' @param xlim Plot's x limits
#' @param ylim Plot's y limits
#' @param loc.standarize Are these standardized vectors
#' @param error.bars Plot error bars
#' @param raw.points Plot raw data points at level of community pairs
#' @param vectors Plot averaged vectors
#' @param legend Show legend
#' @param old.plot ggplot object from previous \code{leap.zig()} call
#' @param add Add new plot to object provided in old.plot option
#'
#' @return A ggplot object.
#'
#' @examples
#'
#' # Load data and run pairwise comparisons of communities
#' cc2<-group_by(cedarcreek,NTrt,NAdd,Plot)
#' pp<-pairwise.price(cc2,species='Species',func='Biomass')
#'
#' # Organize/format the results, and pull out a subset using NTrt=1 as the control site
#' pp<-group.columns(pp,gps=c('NTrt','NAdd'))
#' pp<-pp[pp$NTrt.x==1,]
#' dat1<-pp[pp$NAdd %in% c('0 27.2'),]
#'
#' # Process and plot the result
#' tmp <- process.data.bef(dat1, group.vars=NULL, standardize=F)
#' leap.zig.bef(tmp, loc.standardize=F, group.vars=NULL,raw.points = F,legend = F)
#'
#' @export
#' @import ggplot2
leap.zig.both<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE,group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp[[3]]$variable <- factor(as.character(tmp[[3]]$variable),
levels=c("SL vector","SG vector","CDE vector","SR vector","CE vector"))
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot() + geom_hline(yintercept=0, linetype=2)
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
colour='gray') +
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x, y=mean.y),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp[[3]], aes(colour=variable, x=mean.x, y=mean.y),
arrow=arrow(length=unit(0.2,"cm"), ends="first"))
}
cols <- c(alpha('#d95f02'),alpha('#1b9e77'),alpha('#7570b3'),alpha('#7fcdbb'),alpha('#e34a33'))
trcols <- c(alpha('#d95f02',0.1),alpha('#1b9e77',0.1),alpha('#7fcdbb',0.1),alpha('#7570b3',0.1))
if(raw.points==FALSE & vectors==TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA,NA,NA,NA),
linetype=c(1,1,1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend( override.aes=list(
shape=c(19,1,1,1,1,NA,NA,NA,NA,NA),
linetype=c(0,0,0,0,0,1,1,1,1,1),
colour=c('black',cols[1:4],cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim) == 2){
xlim.outer <- c(min(xlim.D[1], xlim[1]), max(xlim.D[2], xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1], ylim[1]), max(ylim.D[2], ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + theme_bw() + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for BEF components
#' @export
#' @import ggplot2
leap.zig.bef<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable), levels=c("SR vector","CE vector"))
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.bef()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#008002'),alpha('#56F256'))
trcols <- c(alpha('#008002',0.1),alpha('#56F256',0.1))
if(raw.points == FALSE & vectors == TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA),linetype=c(1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend( override.aes=list(shape=c(19,1,1,NA,NA),
linetype=c(0,0,0,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim)==2){
xlim.outer <- c(min(xlim.D[1],xlim[1]), max(xlim.D[2],xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1],ylim[1]), max(ylim.D[2],ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}
if(legend == FALSE){
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for CAFE components
#' @export
#' @import ggplot2
leap.zig.cafe<-function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable),
levels=c("SL vector","SG vector","CDE vector"))
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.cafe()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#FE0000'), alpha('#0025FF'), alpha('#A01AE5'))
trcols <- c(alpha('#FE0000',0.1), alpha('#0025FF',0.1), alpha('#A01AE5',0.1))
if(raw.points==FALSE & vectors==TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend( override.aes=list(shape=c(NA,NA,NA),
linetype=c(1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols))+
guides(colour=guide_legend( override.aes=list(shape=c(19,1,1,1,NA,NA,NA),
linetype=c(0,0,0,0,1,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim)==2){
xlim.outer <- c(min(xlim.D[1],xlim[1]), max(xlim.D[2],xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1],ylim[1]), max(ylim.D[2],ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
#' @describeIn leap.zig.both Plot vectors for Price components
#' @export
#' @import ggplot2
leap.zig.price <- function(tmp, xlim=NA, ylim=NA, loc.standardize=TRUE, error.bars=FALSE,
raw.points=TRUE, vectors=TRUE, group.vars=NULL,
legend=TRUE, old.plot=NA, add=FALSE){
# rename to simplify code:
tmp3<-tmp[[3]]
# Trim out un-needed factor levels
if(raw.points == FALSE & vectors == TRUE){
tmp3$variable <- factor(as.character(tmp3$variable),
levels=c("SRE.L vector","SIE.L vector","SRE.G vector",
"SIE.G vector","CDE vector"))
}
# Plot it:
if(add == TRUE){
lzp <- old.plot
}else{
lzp <- ggplot()
}
# Set up group column?
if(!is.null(group.vars)){
if(length(group.vars)>1){
print("Error in leap.zig.cafe()! Only a single grouping variable implemented")
break()
}
tmp3$gps <- data.frame(tmp3[,group.vars])[,1]
}else{ # with no user-specified grouping variable, create a grouping variable with a single level as dummy column
tmp3$gps <- 1
}
# Add points
if(raw.points){
lzp <- lzp + geom_point(data=tmp[[1]], aes(colour=variable, x=rich, y=value))
}
# Add error bars.
if(error.bars){
lzp <- lzp + geom_errorbarh(data=tmp[[2]], aes(xmin=x.qt.lw, xmax=x.qt.up, x=mean.x, y=mean.y),
height=15, colour='gray')+
geom_errorbar(data=tmp[[2]], aes(ymin=y.qt.lw, ymax=y.qt.up, x=mean.x),
width=1, colour='gray')
}
# Add vectors
if(vectors){
lzp <- lzp + geom_path(data=tmp3, aes(colour=variable, x=mean.x, y=mean.y, group=gps),
arrow=arrow(length=unit(0.2,"cm"), ends="last"))
}
cols <- c(alpha('#DB0102'), alpha('#FF6500'), alpha('#001ECC'),
alpha('#48C5E9'), alpha('#A01AE5'))
trcols <- c(alpha('#DB0102',0.1), alpha('#FF6500',0.1), alpha('#001ECC',0.1),
alpha('#48C5E9',0.1), alpha('#A01AE5',0.1))
if(raw.points == FALSE & vectors == TRUE){
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=cols) +
guides(colour=guide_legend(override.aes=list(shape=c(NA,NA,NA,NA,NA),
linetype=c(1,1,1,1,1))))
}else{
lzp <- lzp + scale_color_manual("Component\n", drop=FALSE, values=c('black',trcols,cols)) +
guides(colour=guide_legend(
override.aes=list(shape=c(19,1,1,1,1,1,NA,NA,NA,NA,NA),
linetype=c(0,0,0,0,0,0,1,1,1,1,1),
colour=c('black',cols,cols))))
}
# Figure out plot ranges, based on data ranges:
xlim.D <- range(tmp[[1]]$rich)*c(0.98,1.02)
ylim.D <- range(tmp[[1]]$value)*c(0.98,1.02)
# If user supplied xlim's
if(length(xlim) == 2){
xlim.outer <- c(min(xlim.D[1], xlim[1]), max(xlim.D[2], xlim[2]))
xlim.inner <- xlim
}else{
xlim.outer <- xlim.inner <- xlim.D
}
# If user supplied ylim's
if(length(ylim) == 2){
ylim.outer <- c(min(ylim.D[1], ylim[1]), max(ylim.D[2], ylim[2]))
ylim.inner <- ylim
}else{
ylim.outer <- ylim.inner <- ylim.D
}
# Adjust inner plotting range
lzp <- lzp + coord_cartesian(xlim=xlim.inner, ylim=ylim.inner)
# Select color and label options:
lzp <- lzp + scale_x_continuous("Species richness", limits=xlim.outer)
if(loc.standardize == TRUE){
lzp <- lzp + scale_y_continuous("% change in EF vs. baseline", limits=ylim.outer) +
geom_hline(yintercept=0, linetype=2)
}else{
lzp <- lzp + scale_y_continuous("Ecosystem function", limits=ylim.outer)
}
if(legend == TRUE){
lzp <- lzp + theme(legend.title=element_blank())
}else{
lzp <- lzp + theme(legend.position="none")
}
return(lzp)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.