# plot.R - ggplot2-based plot methods for FLCore classes
# ggplotFL/R/plot.R
# Copyright 2012-2018 FLR Team. Distributed under the GPL 2
# Maintainer: Iago Mosqueira (EC JRC) <iago.mosqueira@ec.europa.eu
globalVariables("density")
# plot(FLQuant) {{{
#' ggplot versions of FLR class plot() methods
#'
#' New basic plot methods for some FLR classes are defined in ggplotFL.
#'
#' The coertion to *data.frame* that is carried out in the plot methods sets
#' the argument `date=TRUE`. This generates a new column of class `POSIXct` for
#' the first day of the first month of each season. If the `season` dimension of
#' the object being plotted is of length greater than one, `date` will be used
#' as variable on the x axis of the plot. Otherwise, it will be `year`. Keep this
#' in mind when adding extra elements to the plot (see examples below).
#'
#' A similar mechanism is used for the *y* axis, depending on the length of the
#' `iter` dimension. For objects with no *iters*, a single line is plotted for
#' each *FLQuant*, and the *y* axis is mapped to the `data` column of the
#' *data.frame*. For objects with iterations, i.e. with length greater than 1 on
#' the `iter` dimension, the default plots show the quantiles of the distribution
#' and the *y* axis is mapped to the middle quantile, by default `50%`. See the
#' examples below on how to refer to these variables when adding elements to the
#' plot.
#'
#' @param x FLR object to plot
#' @param y FLR object to plot
#' @param na.rm Should NAs be deleted in quantile calculations?, defaults to TRUE.
#' @param probs Quantiles to be plotted if object has iters, defaults to c(0.10, 0.33, 0.50, 0.66, 0.90).
#' @param ... Other arguments to be passed to the corresponding ggplot call.
#' @param metrics FlQuants computed from complex objects (e.g. FLStock)
#' @param iter Individual iterations to show as worm plots over the quantiles.
#' @param worm Individual iterations to show as worm plots over the quantiles.
#'
#' @aliases plot,FLQuant,missing-method
#' @seealso \code{\link{ISOdate}} \code{\link{ggplot}}
#' @docType methods
#' @rdname plot
#' @examples
#'
#' # Plot a single FLQuant
#' data(ple4)
#' plot(catch.n(ple4))
#'
#' # Plot an FLQuant with iters, shows quantiles
#' flq <- rnorm(100, catch(ple4), 60000)
#' plot(flq)
#'
#' # Specify quantiles, default is c(0.10, 0.33, 0.50, 0.66, 0.90)
#' plot(flq, probs=c(0.05, 0.25, 0.50, 0.75, 0.95))
#'
#' # Adding extra elements to an FLQuant plot, with seasons
#' flq <- FLQuant(runif(200), dim=c(1,15,1,4))
#' plot(flq) + geom_point(aes(x=date, y=data, colour=season))
#'
#' # or without them
#' flq <- FLQuant(runif(200), dim=c(1,15))
#' plot(flq) + geom_point(aes(x=year, y=data))
#'
#' # For an object with iter
#' flq <- rlnorm(100, flq, 0.4)
#' plot(flq) + geom_point(aes(x=year, y=data))
#'
#' # To plot(FLQuant) as in previous versions of ggplotFL
#' plot(rnorm(300, catch(ple4), catch(ple4)/2), probs=c(0.10, 0.5, 0.90)) +
#' geom_flquantiles(probs=c(0.01), linetype=3, colour="red", alpha=0.1) +
#' geom_flquantiles(probs=c(0.99), linetype=3, colour="red", alpha=0.1)
setMethod("plot", signature(x="FLQuant", y="missing"),
function(x, probs=c(0.05, 0.25, 0.50, 0.75, 0.95), na.rm=FALSE, ...) {
# GET base plot from plot(FLQuants)
p <- plot(FLQuants(x), probs=probs, na.rm=na.rm, ...)
# PARSE dimensions > 1 for new facets
ldi <- c(names(x)[-c(2,3,4,6)][dim(x)[-c(2,3,4,6)] > 1])
# PLOT(FLQuants(x)) & reset facets
if(length(ldi) == 0) {
p <- p + theme(strip.text.y = element_blank())
}
else if(length(ldi) == 1) {
p <- p + facet_grid(as.formula(paste0(ldi, "~.")), scales="free",
labeller=label_both)
}
else if (length(ldi) > 1) {
p <- p + facet_grid(as.formula(paste0(ldi[1], "~", paste(ldi[-1],
collapse= "+"))), scales="free", labeller=label_both)
}
return(p)
})
# }}}
# plot(FLQuant, FLQuant) {{{
#' @rdname plot
#' @examples
#' # plot(FLQuant, FLQuant, ...) to place in one facet
#' plot(catch(ple4), landings(ple4))
#' # Add legend by hand
#' plot(rnorm(200, landings(ple4), 8000), discards(ple4)) +
#' scale_colour_discrete(name="Yield (t)", labels=c("Landings", "Discards")) +
#' theme(legend.position="bottom")
setMethod("plot", signature(x="FLQuant", y="FLQuant"),
function(x, y, ..., probs=c(0.05, 0.25, 0.50, 0.75, 0.95), na.rm=FALSE, iter=NULL) {
# ASSEMBLE FLQuants
fqs <- FLQuants(c(list(x, y), list(...)))
# PLOT as
ggplot(fqs, aes(x=year, y=data, fill=qname, colour=qname)) +
geom_flquantiles(alpha=0.3, probs=probs, na.rm=na.rm) +
theme(legend.position="none") +
xlab("") + ylab("")
}
)
# }}}
# plot(FLQuants) {{{
#' @aliases plot,FLQuants,missing-method
#' @rdname plot
#' @examples
#' # Plot an FLQuants created from ple4 FLStock
#' data(ple4)
#' plot(FLQuants(SSB=ssb(ple4), rec=rec(ple4)))
#' plot(FLQuants(SSB=ssb(ple4), rec=rec(ple4)), probs = NULL)
setMethod("plot", signature(x="FLQuants", y="missing"),
function(x, probs=c(0.05, 0.25, 0.50, 0.75, 0.95), na.rm=FALSE, worm=iter,
iter=NULL) {
# SWITCH off ribbons
if(is.null(probs)) {
probs <- rep(0, 3)
if(missing(worm))
worm <- TRUE
} else if(length(probs) == 1) {
probs <- rep(probs, 3)
}
# CHECK probs length is odd
if(is.integer(length(probs)/2))
stop("quantile probs can only be a vector of odd length")
# FIND center of probs
idx <- ceiling(length(probs)/2)
# GET max dimensions
mds <- apply(do.call(rbind, lapply(x, dim)), 2, max)
# USE year or date for x axis
xvar <- sym(ifelse(mds[4] == 1, "year", "date"))
if(mds[6] == 1) {
# NO ITERS? PLOT central line by unit
p <- if(mds[3] == 1) {
ggplot(x, aes(x=!!xvar, y=data)) +
geom_line(na.rm=na.rm)
} else {
ggplot(x, aes(x=!!xvar, y=data, fill=unit, colour=unit)) +
geom_line(na.rm=na.rm)
}
} else {
# worm=TRUE? PLOT central ribbon and line by unit
if(isTRUE(worm)) {
p <- if(mds[3] == 1) {
ggplot(x, aes(x=!!xvar, y=data, fill=flpalette_colours(1))) +
geom_line(aes(group=iter), alpha=0.2, linewidth=1, colour="#adadad") +
geom_flquantiles(alpha=0.5,
probs=probs[seq(idx - 1, idx + 1)], na.rm=na.rm)
} else {
ggplot(x, aes(x=!!xvar, y=data, fill=unit, colour=unit)) +
geom_line(aes(group=iter), alpha=0.1, linewidth=1) +
geom_flquantiles(alpha=0.5,
probs=probs[seq(idx - 1, idx + 1)], na.rm=na.rm)
}
} else {
p <- if(mds[3] == 1) {
ggplot(x, aes(x=!!xvar, y=data, fill=flpalette_colours(1))) +
geom_flquantiles(alpha=0.5,
probs=probs[seq(idx - 1, idx + 1)], na.rm=na.rm)
} else {
ggplot(x, aes(x=!!xvar, y=data, fill=unit, colour=unit)) +
geom_flquantiles(alpha=0.5,
probs=probs[seq(idx - 1, idx + 1)], na.rm=na.rm)
}
}
}
# PLOT other ribbons, if any
if(length(probs) > 3 & mds[6] > 1) {
geoms <- lapply(seq((length(probs)-3)/2), function(x) {
geom_flquantiles(probs=probs[seq(idx - x - 1, idx + x + 1)],
alpha=0.3)
})
p <- p + geoms
}
# SHOW NAs in x axis, only if no iters
if(mds[6] == 1 & na.rm == FALSE) {
if(any(is.na(p$data)))
p <- p + geom_point(aes(y=0), cex=0.6, colour='darkgrey',
data=subset(p$data, is.na(data)))
}
# PLOT iter worms
if(is.numeric(worm) | is.character(worm)) {
idata <- subset(p$data, iter %in% worm)
# SELECT by position if numbers
if(is.numeric(worm)) {
# FIND iter dimnames for given positions
ma <- data.table(p$data)[, .(iter=unique(iter)[worm]), by=qname]
# SUBSET for combinations
idata <- ma[data.table(p$data), , nomatch=0L, on=c("qname", "iter")]
# RENAME iters so colours (factor) match
idata[, iter:=as.character(rep(rep(seq(length(worm)),
each=uniqueN(idata$year)), uniqueN(idata$qname)))]
}
p <- p + geom_line(data=idata, aes(x=!!xvar, y=data, colour=iter),
na.rm=TRUE)
}
# BUILD facet formula
ldi <- c("qname", names(x[[1]])[-c(2, 3, 4, 6)][mds[-c(2, 3, 4, 6)] > 1])
if(length(ldi) == 1) {
p <- p + facet_grid(as.formula(paste0(ldi, "~.")), scales="free",
labeller=label_flqs(x))
}
else if (length(ldi) > 1) {
p <- p + facet_grid(as.formula(paste0(ldi[1], "~", paste(ldi[-1],
collapse= "+"))), scales="free", labeller=label_flqs(x))
}
# ASSEMBLE plot
p <- p + xlab("") + ylab("") + theme(legend.position="none")
return(p)
})
# }}}
# plot(FLQuants, FLPar) {{{
#' @aliases plot,FLQuants,FLPar-method
#' @rdname plot
#' @examples
#' # plot for FLQuants, FLPar
#' data(ple4)
#' rps <- FLPar(F=0.14, Catch=1.29e5, Rec=9.38e5, SSB=1.8e5)
#' fqs <- metrics(ple4)
#' plot(fqs, rps)
#' # Works also if reptsa are given for some panels
#' rps <- FLPar(F=0.14, Catch=1.29e5, SSB=1.8e5)
#' plot(fqs, rps)
# TODO REFORMAT
setMethod("plot", signature(x="FLQuants", y="FLPar"),
function(x, y, ...) {
p <- plot(x)
# GET name variable mapped to y axis
dat <- quo_name(p$mapping$y)
# CREATE df with right name
rpa <- data.frame(dat=c(y), qname=dimnames(y)$params, stringsAsFactors=FALSE)
colnames(rpa)[1] <- dat
# FIX mixmatch between refpts and FLStock slots naming
if('yield' %in% rpa$qname)
rpa$qname[rpa$qname == 'yield'] <- 'catch'
p <- p + geom_hline(data=rpa, aes(yintercept=data), colour="blue", linetype=2)
return(p)
}
) # }}}
# plot (FLQuants, FLPars) {{{
#' @aliases plot,FLQuants,FLPars-method
#' @rdname plot
#' @examples
#' # plot for FLQuants, FLPars
#' data(ple4)
#' rps <- FLPars(F=FLPar(Fmsy=0.14, Fpa=0.35), SSB=FLPar(SBmsy=1.8e5, SBlim=1.1e5))
#' fqs <- metrics(ple4, list(SSB=ssb, F=fbar))
#' plot(fqs, rps) + ylim(c(0, NA))
setMethod("plot", signature(x="FLQuants", y="FLPars"),
function(x, y, ...) {
p <- plot(x)
# GET name variable mapped to y axis
dnm <- quo_name(p$mapping$y)
# CREATE df with right names
dat <- do.call(rbind, c(mapply(function(i, j)
cbind(as.data.frame(i, drop=TRUE), qname=j), y, names(y),
SIMPLIFY=FALSE), make.row.names = FALSE))
colnames(dat)[2] <- dnm
# DEBUG factors to character
dat[,"params"] <- as.character(dat[,"params"])
# FIX common mixmatch between refpts and FLStock slots naming
if('yield' %in% dat$qname)
dat$qname[dat$qname == 'yield'] <- 'catch'
# MERGE refpts with same value
counts <- table(dat$data)
dups <- counts[counts > 1]
if(length(dups) > 0) {
idx <- lapply(names(dups), function(x) which(dat$data == x))
dat[unlist(lapply(idx, '[', 1)), "params"] <-
unlist(lapply(idx, function(x) paste(as.character(dat[x, "params"]),
collapse=" - ")))
# DROP merged params
dat <- dat[-unlist(lapply(idx, '[', -1)),]
}
# SET y nudge up
lim <- do.call(rbind, mapply(function(i, j) data.frame(max=max(i),
min=min(i), qname=as.character(j)), x, names(x), SIMPLIFY=FALSE))
dat <- merge(dat, lim)
p <- p + geom_hline(data=dat, aes(yintercept=data), linetype=2,
linewidth=0.25) +
geom_text(data=dat, aes(y=data + ((max-min) * 0.05), label=params),
x=dims(x[[1]])$minyear - 1, size=3, hjust="inward")
return(p)
}
) # }}}
# plot(FLQuantPoint) {{{
#' @aliases plot,FLQuantPoint,missing-method
#' @rdname plot
#' @examples
#' # plot for FLQuantPoint
#' fqp <- FLQuantPoint(rlnorm(300, log(catch(ple4)), 0.20))
#' plot(fqp)
setMethod("plot", signature(x="FLQuantPoint", y="missing"),
function(x, mean=TRUE, median=TRUE) {
# BASE plot
p <- ggplot(x, aes(x=date))
# ELEMENTS to add
if(!all(is.na(p$data$median)) & mean)
p <- p + geom_line(aes(y=median), linetype=1)
if(!all(is.na(p$data$mean)) & median)
p <- p + geom_line(aes(y=mean), linetype=2)
if(!all(is.na(p$data[, c("lowq", "uppq")])))
p <- p + geom_ribbon(aes(ymin=lowq, ymax=uppq), colour="gray",
alpha=0.20, linetype=3)
# PARSE dimensions > 1 for new facets
ldi <- c(names(x)[-c(2,3,4,6)][dim(x)[-c(2,3,4,6)] > 1])
# PLOT(FLQuants(x)) & reset facets
if(length(ldi) == 0) {
p <- p + theme(strip.text.y = element_blank())
}
else if(length(ldi) == 1) {
p <- p + facet_grid(as.formula(paste0(ldi, "~.")), scales="free",
labeller=label_both)
}
else if (length(ldi) > 1) {
p <- p + facet_grid(as.formula(paste0(ldi[1], "~", paste(ldi[-1],
collapse= "+"))), scales="free", labeller=label_both)
}
# ASSEMBLE plot
p <- p + xlab("") + ylab("")
return(p)
}
) # }}}
# plot(FLQuantPoint, FLQuant) {{{
#' @aliases plot,FLQuantPoint,FLQuant-method
#' @rdname plot
#' @examples
#' # plot for FLQuantPoint, FLQuant
#' plot(fqp, rlnorm(3, log(catch(ple4)), 0.20))
setMethod("plot", signature(x="FLQuantPoint", y="FLQuant"),
function(x, y, na.rm=FALSE, ...) {
if(dim(y)[6] > 1)
plot(x, divide(y), ...)
else
plot(x, FLQuants(y=y), ...)
}
) # }}}
# plot(FLQuantPoint, FLQuants) {{{
#' @aliases plot,FLQuantPoint,FLQuants-method
#' @rdname plot
#' @examples
#' # plot for FLQuantPoint, FLQuants
#' fqp <- FLQuantPoint(rlnorm(300, log(catch(ple4)), 0.20))
#' fqs <- divide(rlnorm(3, log(catch(ple4)), 0.20))
#' plot(fqp, fqs)
setMethod("plot", signature(x="FLQuantPoint", y="FLQuants"),
function(x, y, na.rm=FALSE, mean=TRUE, median=TRUE, ...) {
# PLOT FLQuantPoint
p <- plot(x, mean=mean, median=median)
# EXTRACT FLQuants
dat <- as.data.frame(y, drop=TRUE, date=TRUE)
# ADD FLQuants as lines
p <- p + geom_flquantiles(data=dat, aes(x=date, y=data, colour=qname))
# PARSE dimensions > 1 for new facets
ldi <- c(names(x)[-c(2,3,4,6)][dim(x)[-c(2,3,4,6)] > 1])
# PLOT(FLQuants(x)) & reset facets
if(length(ldi) == 0) {
p <- p + theme(strip.text.y = element_blank())
}
else if(length(ldi) == 1) {
p <- p + facet_grid(as.formula(paste0(ldi, "~.")), scales="free",
labeller=label_both)
}
else if (length(ldi) > 1) {
p <- p + facet_grid(as.formula(paste0(ldi[1], "~", paste(ldi[-1],
collapse= "+"))), scales="free", labeller=label_both)
}
# ASSEMBLE plot
p <- p + xlab("") + ylab("") + theme(legend.position="none")
return(p)
})
# }}}
# plot(FLPar, missing) {{{
#' @rdname plot
#' @examples
#' par <- FLPar(alpha=rnorm(200, 0.6, 0.2), beta=rlnorm(200, 0.8, 0.3))
#' plot(par)
setMethod("plot", signature(x="FLPar", y="missing"),
function(x, names=NULL) {
# EXTRACT units
ups <- units(x)
ggplot(as.data.frame(x), aes(x=data)) +
geom_density(alpha=.2, aes(fill=params)) +
geom_histogram(aes(y=after_stat(density)), bins=20,
colour="black", fill=NA) +
facet_wrap(~params, scales="free", labeller=
format_label_flqs(ups, dimnames(x)$params)) +
xlab("") + ylab("") +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
scale_x_continuous(n.breaks = 7) +
guides(fill="none")
}
)
# }}}
# plot(FLStock) {{{
#' @aliases plot,FLStock,missing-method
#' @rdname plot
#' @examples
#' # plot of an FLStock
#' data(ple4)
#' plot(ple4)
setMethod("plot", signature(x="FLStock", y="missing"),
function(x, metrics=list(Rec=rec, SSB=ssb, Catch=catch, F=fbar), na.rm=TRUE,
...) {
# PLOT unfitted
if(all(is.na(harvest(x))) & missing(metrics))
metrics <- list(Catch=catch, Landings=landings, Discards=discards)
metrics <- metrics(x, metrics=metrics)
# HACK for F units
if("F" %in% names(metrics))
units(metrics$F) <- paste0(range(x, c("minfbar", "maxfbar")),
collapse="-")
# ADAPT for 2-sex model
if("SSB" %in% names(metrics)) {
if(all(dimnames(metrics$SSB)$unit %in% c("F", "M"))) {
# FIND spawning season, if it exists
#if(dim(metrics$SSB)[4] > 1) {
# metrics$SSB[is.na(metrics$SSB)] <- 0
#}
# DROP M ssb if missing
metrics$SSB <- metrics$SSB[,,'F'] + metrics$SSB[,,'M']
# SUM rec across units
if("Rec" %in% names(metrics))
metrics$Rec <- unitSums(metrics$Rec)
}
}
# ADAPT for seasonal recruitment
# if("Rec" %in% names(metrics)) {
# if(dim(metrics$Rec)[4] > 1) {
# metrics$Rec[metrics$Rec == 0] <- NA
# }
# }
p <- plot(metrics, na.rm=na.rm, ...) + ylim(c(0, NA))
# ADD legend if 2 sexes
if("SSB" %in% names(metrics))
if(all(dimnames(metrics$SSB)$unit %in% c("F", "M"))) {
return(p +
theme(legend.position="bottom", legend.key=element_blank()) +
labs(color="Sex") +
scale_color_manual(name="",
labels=c("Both", "F", "M"),
values=flpalette_colours(3))
)
}
return(p)
}
) # }}}
# plot(FLStock, FLStock) {{{
#' @aliases plot,FLStock,FLStock-method
#' @rdname plot
setMethod("plot", signature(x="FLStock", y="FLStock"),
function(x, y, metrics=list(Rec=rec, SSB=ssb, Catch=catch, F=fbar),
probs=c(0.10, 0.33, 0.50, 0.66, 0.90), na.rm=TRUE, iter=NULL, ...) {
args <- list(...)
sts <- do.call("FLStocks", c(list(x, y), args))
names(sts) <- unlist(lapply(sts, name))
p <- plot(sts, probs=probs, na.rm=na.rm, iter=iter, metrics=metrics)
return(p)
}
) # }}}
# plot(FLStock, FLPar) {{{
#' @aliases plot,FLStock,FLPar-method
#' @rdname plot
#' @examples
#' # plot for FLStock, FLPar
#' data(ple4)
#' rps <- FLPar(F=0.14, Catch=1.29e5, Rec=9.38e5, SSB=1.8e5)
#' plot(ple4, rps)
setMethod("plot", signature(x="FLStock", y="FLPar"),
function(x, y, metrics=list(Rec=rec, SSB=ssb, Catch=catch, F=fbar), ...) {
p <- plot(metrics(x, metrics=metrics), y, ...)
return(p)
}
) # }}}
# plot(FLStocks) {{{
#' @aliases plot,FLStocks,missing-method
#' @rdname plot
#' @param metrics function returning an FLQuants for each FLStock
#' @param probs Quantiles to calculate along the iter dimension. A vector of length 5, for the lower outer, lower inner, central, upper inner and upper outer quantiles. Defaults to the 66 and 80 percent quantiles, plus median line.
#' @param alpha alpha values for the quantile ribbons, defaults to 0.10 and 0.40.
#' @examples
#' # plot for FLStocks
#' data(ple4)
#' pls <- FLStocks(runA=ple4, runB=qapply(ple4, function(x) x*1.10))
#' plot(pls)
#' # geom_flpar can be used draw refpts lines and labels
#' plot(pls, metrics=list(SSB=ssb, F=fbar)) +
#' facet_grid(qname~stock, scales='free') +
#' geom_flpar(data=FLPars(SSB=FLPar(Blim=300000, Bpa=230000),
#' F=FLPar(FMSY=0.21)), x=c(1960), stock='runA', fill=alpha('white', 0.4))
setMethod("plot", signature(x="FLStocks", y="missing"),
function(x, metrics=list(Rec=rec, SSB=ssb, Catch=catch, F=fbar),
probs=c(0.10, 0.33, 0.50, 0.66, 0.90), alpha=c(0.10, 0.40), worm=iter,
iter=NULL, ...) {
# CHECK names not repeated
dup <- duplicated(names(x))
if(any(dup)) {
names(x)[dup] <- paste(names(x)[dup], LETTERS[seq(sum(dup))], sep='_')
# warning('Duplicated names in object, changed to differentiate')
}
# EXTRACT slots by stock
fqs <- lapply(x, "metrics", metrics)
# HACK for F units
if("F" %in% names(fqs[[1]]))
fqs <- lapply(fqs, function(fq) {
units(fq$F) <- paste0(range(x[[1]],
c("minfbar", "maxfbar")), collapse="-")
return(fq)
})
# GET labels
labeller <- label_flqs(fqs[[1]])
# ASSEMBLE data
data <- lapply(fqs, as.data.frame, date=TRUE, drop=FALSE)
# SET stock names
stk <- rep.int(names(fqs), unlist(lapply(data, nrow)))
# RBIND dfs
data <- do.call(rbind, data)
rownames(data) <- NULL
# USE year or date for x axis
xvar <- sym(ifelse(length(unique(data$season)) == 1, "year", "date"))
# ADD stock names
data <- transform(data, stock=factor(stk, levels=names(x)))
# PLOT using geom_flquantiles
p <- ggplot(data, aes(x=!!xvar, y=data, fill=stock, colour=stock)) +
facet_grid(qname~., labeller=labeller, scales="free_y") +
# outer quantile
geom_flquantiles(probs=probs[c(1, 5)], alpha=alpha[1],
colour="white") +
# inner quantile
geom_flquantiles(probs=probs[c(2, 4)], alpha=alpha[2],
colour="white") +
# median
geom_flquantiles(probs=probs[3], alpha=1) +
xlab("") + ylab("") +
# SET limits to include 0
expand_limits(y=0) +
# SET legend with no title
theme(legend.title = element_blank()) +
# and only with lines and no title
guides(fill = "none")
# PLOT iter worms
if(is.numeric(worm)) {
idata <- p$data[p$data$iter %in% worm,]
p <- p + geom_line(data=idata, aes(x=!!xvar, y=data, colour=iter))
}
return(p)
}
) # }}}
# plot(FLStocks) {{{
#' @aliases plot,FLStocks,missing-method
#' @rdname plot
#' @param metrics function returning an FLQuants for each FLStock
#' @param probs Quantiles to calculate along the iter dimension. A vector of length 5, for the lower outer, lower inner, central, upper inner and upper outer quantiles. Defaults to the 66 and 80 percent quantiles, plus median line.
#' @param alpha alpha values for the quantile ribbons, defaults to 0.10 and 0.40.
#' @examples
#' # plot for FLStocks
#' data(ple4)
#' pls <- FLStocks(runA=ple4, runB=qapply(ple4, function(x) x*1.10))
#' plot(pls)
#' # geom_flpar can then be used draw refpts lines and labels
#' plot(pls, metrics=list(SSB=ssb, F=fbar)) +
#' facet_grid(qname~stock, scales='free') +
#' geom_flpar(data=FLPars(SSB=FLPar(Blim=300000, Bpa=230000),
#' F=FLPar(FMSY=0.21)), x=c(1960), stock='runA', fill=alpha('white', 0.4))
setMethod("plot", signature(x="FLStocks", y="missing"),
function(x, metrics=list(Rec=function(x) unitSums(rec(x)),
SB=function(x) unitSums(ssb(x)), C=function(x) unitSums(catch(x)),
F=function(x) unitMeans(fbar(x))),
probs=c(0.10, 0.33, 0.50, 0.66, 0.90), alpha=c(0.10, 0.40), worm=iter,
iter=NULL, ...) {
# EXTRACT slots by stock
fqs <- lapply(x, "metrics", metrics=metrics)
p <- plotListFLQuants(fqs, probs=probs, alpha=alpha, worm=worm, iter=iter)
return(p)
}
) # }}}
# plot(FLStocks, FLPar) {{{
# TODO Move to geom_flquantiles
#' @aliases plot,FLStocks,FLPar-method
#' @rdname plot
setMethod("plot", signature(x="FLStocks", y="FLPar"),
function(x, y, na.rm=TRUE,
metrics= function(x, y) FLQuants(SSB=ssb(x)/y[,'ssb',], F=fbar(x)/y[,'harvest',],
Catch=catch(x))) {
# check names not repeated
dup <- duplicated(names(x))
if(any(dup)) {
names(x)[dup] <- paste(names(x)[dup], LETTERS[seq(sum(dup))], sep='_')
warning('Duplicated names in object, changed to differentiate')
}
# extract slots by stock
fqs <- lapply(x, metrics, y)
# get median & 85% quantiles if iters
its <- unlist(lapply(x, function(x) dims(x)$iter))
if(any(its > 1))
{
# quantiles
fqs <- lapply(fqs, function(y) as.data.frame(lapply(y, quantile,
c(0.10, 0.50, 0.90), na.rm=TRUE)))
} else {
fqs <- lapply(fqs, as.data.frame)
fqs <- lapply(fqs, function(x) {x$iter <- "50%"; return(x)})
}
# stock names
stk <- rep.int(names(fqs), unlist(lapply(fqs, nrow)))
# rbind dfs
fqs <- do.call(rbind, fqs)
rownames(fqs) <- NULL
# add stock names
fqs <- transform(fqs, stock=stk)
# compute quantiles
df <- reshape(fqs, timevar="iter", direction="wide",
idvar=c(names(fqs)[1:5], "qname", "stock", "date"))
names(df) <- gsub("data.", "", names(df))
# plot data vs. date + facet on qname +
p <- ggplot(data=na.omit(df), aes_string(x='date', y='`50%`', group='stock')) +
facet_grid(qname~., scales="free") +
# line + xlab + ylab +
geom_line(aes_string(colour='stock'), na.rm=na.rm) + xlab("") + ylab("") +
# limits to include 0 + no legend
expand_limits(y=0) + theme(legend.title = element_blank())
# object w/ iters?
if(any(unlist(lapply(x, function(y) dims(y)$iter)) > 1)) {
p <- p +
# 75% quantile ribbon in red, alpha=0.25
geom_ribbon(aes_string(x='date', ymin = '`10%`', ymax = '`90%`', group='stock',
colour='stock', fill='stock'), alpha = .20, linetype = 0, na.rm=na.rm)
# 90% quantile ribbon in red, aplha=0.10
}
return(p)
}
) # }}}
# plot(FLStock, FLStocks) {{{
#' @aliases plot,FLStock,FLStocks,missing-method
#' @rdname plot
setMethod("plot", signature(x="FLStock", y="FLStocks"),
function(x, y, ...) {
plot(FLStocks(c(x, y)), ...)
}
) # }}}
# plot(FLSR) {{{
#' @aliases plot,FLSR,missing-method
#' @docType methods
#' @rdname plot
#' @examples
#' # plot for FLSR
#' data(nsher)
#' plot(nsher)
setMethod('plot', signature(x='FLSR', y='missing'),
function(x, ...) {
dat <- model.frame(FLQuants(SSB=ssb(x), Rec=rec(x), Residuals=residuals(x),
RecHat=fitted(x)))
uns <- units(x)
unr <- ifelse(uns$rec == 'NA', 'Recruits', as.expression(paste0('Recruits (',
sub('*', 'A', uns$rec, fixed=TRUE), ')')))
uns <- ifelse(uns$ssb == 'NA', 'SSB', as.expression(paste0('SSB (', sub('*',
'%*%', uns$ssb, fixed=TRUE), ')')))
# SSB vs. REC
p1 <- ggplot(data=na.omit(dat), aes_string(x='SSB', y='Rec')) + geom_point() +
geom_smooth(formula= y ~ x, method='loess', span=3) + xlab(uns) + ylab(unr) +
expand_limits(y=0) + expand_limits(x=0)
# model fit line
form <- as.list(model(x))[[3]]
pars <- as(params(x), 'list')
fmo <- function(x) {
c(eval(form, c(list(ssb=FLQuant(x)), pars)))
}
p1 <- p1 + stat_function(fun=fmo, colour='red', linewidth=0.5)
# P2
p2 <- ggplot(data=na.omit(dat), aes_string(x='year', y='Residuals')) +
geom_point() + geom_smooth(formula= y ~ x, method='loess', span=3) + xlab("Year")
# P3
p3 <- ggplot(data=na.omit(data.frame(res1=dat$Residuals[-length(dat$Residuals)],
res2=dat$Residuals[-1])), aes_string(x='res1', y='res2')) + geom_point() +
xlab(expression(Residuals[t])) + ylab(expression(Residuals[t + 1])) +
geom_smooth(formula= y ~ x, method='lm')
# P4
p4 <- ggplot(data=na.omit(dat), aes_string(x='SSB', y='Residuals')) +
geom_point() + geom_smooth(formula= y ~ x, method='loess', span=3)
# P5
p5 <- ggplot(data=na.omit(dat), aes_string(sample = 'Residuals')) +
stat_qq(color="red", alpha=1) +
geom_abline(aes_q(intercept = quote(mean(Residuals)),
slope = quote(sd(Residuals)))) + xlab("Theoretical") + ylab("Sample")
# P6
p6 <- ggplot(data=na.omit(dat), aes_string(x='RecHat', y='Residuals')) +
geom_point() + geom_smooth(formula= y ~ x, method='loess', span=3) +
xlab(expression(hat(Recruits)))
# ASSEMBLE grid and CONVERT to gg
res <- arrangeGrob(p1, p2, p3, p4, p5, p6, ncol=2)
return(ggdraw() + draw_grob(grid::grobTree(res)))
# return(p)
}
) # }}}
# plot(FLSRs) {{{
#' @aliases plot,FLSRs,missing-method
#' @docType methods
#' @rdname plot
#' @param legend_label function to create the legend labels
#' @param facets
#' @examples
#' # plot for FLSRs
#' data(nsher)
#' srs <- FLSRs(sapply(c('segreg', 'bevholt'), function(x) {
#' y <- nsher
#' model(y) <- x
#' return(fmle(y))
#' }))
#' plot(srs, facets=TRUE)
#' plot(srs, legend_label=eqlabel)
#' plot(srs, legend_label=modlabel)
setMethod("plot", signature(x="FLSRs"),
function(x, legend_label=names(x), facets=FALSE, ...) {
uns <- units(x[[1]])
# DIFFERENT data?
# if(all(unlist(lapply(x[-1],
# function(y) isTRUE(all.equal(rec(y), rec(x[[1]]))))))) {
# dat <- cbind(sr=NA, model.frame(FLQuants(ssb=ssb(x[[1]]),
# rec=rec(x[[1]]))))
# } else {
dat <- Reduce(rbind, Map(function(x, i)
cbind(sr=i, model.frame(FLQuants(ssb=ssb(x), rec=rec(x)), drop=TRUE)),
x, names(x)))
# }
# EXTRACT models & pars
mods <- lapply(x, 'model')
pars <- lapply(x, 'params')
ssbs <- seq(0, max(dat$ssb), length=100)
# RESULTS
res <- lapply(names(mods), function(i) {
if(facets)
ssbs <- seq(0, max(dat[dat$sr == i, 'ssb']), length=100)
data.frame(sr=i, ssb=ssbs, rec=c(eval(as.list(mods[[i]])[[3]],
c(list(ssb=ssbs), as(pars[[i]], 'list'))))
)
})
#
if(!is(legend_label, 'function')) {
legend_label <- function(model, params)
return(setNames(nm=names(model)))
}
res <- Reduce('rbind', res)
# GET plot
p <- ggplot(na.omit(res), aes(x=ssb, y=rec, colour=sr)) +
geom_line(aes(group=sr, color=sr)) +
geom_point(data=dat) +
xlab(as.expression(paste0("SSB (", sub('\\*', '%.%', uns$ssb), ")"))) +
ylab(as.expression(paste0("Recruits (", sub('\\*', '%.%', uns$rec), ")"))) +
scale_color_discrete(name="", breaks=names(x),
labels=do.call(legend_label, list(model=mods, param=pars))) +
theme(legend.position="bottom") +
guides(color=guide_legend(nrow=length(mods), byrow=TRUE))
return(p)
}
) # }}}
# plot(FLBiol) {{{
#' @aliases plot,FLBiol,missing-method
#' @docType methods
#' @rdname plot
setMethod("plot", signature(x="FLBiol", y="missing"),
function(x, metrics=list(Rec=function(x) n(x)[1,], B=tsb), ...) {
flqs <- metrics(x, metrics)
p <- plot(flqs) + ylim(c(0,NA))
# TODO ADD SRR
# TODO ADD mat, fec, m, wt by age
return(p)
}
) # }}}
# plot(FLBiols) {{{
#' @aliases plot,FLBiols,missing-method
#' @docType methods
#' @rdname plot
setMethod("plot", signature(x="FLBiols", y="missing"),
function(x, metrics=list(Rec=function(x) n(x)[1,], B=tsb), ...) {
fqs <- lapply(x, function(x) metrics(x, metrics))
# GET labels
labeller <- label_flqs(fqs[[1]])
dfs <- lapply(fqs, as.data.frame, date=TRUE, units=TRUE)
data <- do.call("rbind", c(mapply(`[<-`, dfs, "biol", value=names(fqs),
SIMPLIFY=FALSE), list(make.row.names = FALSE)))
# USE year or date for x axis
xvar <- sym(ifelse(length(unique(data$season)) == 1, "year", "date"))
# PLOT using geom_flquantiles
p <- ggplot(data, aes(x=!!xvar, y=data, fill=.data$biol, colour=.data$biol)) +
facet_grid(qname~., labeller=labeller, scales="free_y") +
geom_flquantiles() + xlab("") + ylab("") +
# SET limits to include 0
expand_limits(y=0) +
# SET legend with no title
theme(legend.title = element_blank()) +
# and only with lines and no title
guides(fill = "none")
return(p)
}
) # }}}
# plot(FLIndexBiomass) {{{
#' @aliases plot,FLIndexBiomass,missing-method
#' @docType methods
#' @rdname plot
setMethod("plot", signature(x="FLIndexBiomass", y="missing"),
function(x, ...) {
flqs <- FLQuants(Index=index(x))
p <- plot(flqs, ...) + geom_smooth(formula=y ~ x, na.rm=TRUE, method="loess")
return(p)
}
) # }}}
# plot(FLIndex) {{{
# TODO ADD index.var/se
#' @aliases plot,FLIndex,missing-method
#' @docType methods
#' @rdname plot
#' @examples
#' # Plot a FLIndex object
#' data(ple4.index)
#' plot(ple4.index)
setMethod("plot", signature(x="FLIndex", y="missing"),
function(x) {
if(!units(index(x)) %in% c("NA", ""))
ylab <- paste0("Abundance (", units(index(x)), ")")
else
ylab <- "Abundance"
ggplot(index(x), aes(x=date, y=data)) +
geom_line() +
facet_wrap(~age, scales="free_y") +
xlab("") + ylab(ylab)
}
) # }}}
# plot(FLIndices) {{{
#' @aliases plot,FLIndices,missing-method
#' @docType methods
#' @rdname plot
#' @examples
#' # Plot a FLIndices object
#' data(ple4.indices)
#' plot(ple4.indices)
#' plot(ple4.indices) +
#' geom_smooth(formula=y ~ x, se=FALSE, method="loess", linewidth=0.2)
setMethod("plot", signature(x="FLIndices", y="missing"),
function(x) {
fqs <- lapply(x, function(x) (index(x) %-% yearMeans(index(x)) %/%
sqrt(yearVars(index(x)))))
aes_(quote(mpg), quote(wt), col = quote(cyl))
# CHOOSE xvar = date if seasons
if(dim(fqs[[1]])[4] > 1)
xvar <- quote(date)
else
xvar <- quote(year)
p <- ggplot(fqs, aes_(x=xvar, y=quote(data), group=quote(qname),
colour=quote(qname), fill=flpalette_colours(quote(qname)))) +
geom_flquantiles(alpha=0.3) +
ylab("Standardized relative abundance") + xlab("") +
theme(legend.title=element_blank())
if(all(unlist(lapply(x, is, "FLIndexBiomass")))) {
return(p)
} else {
return(p + facet_wrap(~age, scales="free_y"))
}
}
) # }}}
# plotListFLQuants {{{
plotListFLQuants <- function(x, probs=c(0.10, 0.33, 0.50, 0.66, 0.90),
alpha=c(0.10, 0.40), worm=iter, iter=NULL, fages=NULL) {
# CHECK names not repeated
dup <- duplicated(names(x))
if(any(dup)) {
names(x)[dup] <- paste(names(x)[dup], LETTERS[seq(sum(dup))], sep='_')
}
# HACK for F units
if("F" %in% names(x[[1]]))
x <- lapply(x, function(i) {
units(i$F) <- paste0(fages, collapse="-")
return(i)
})
# GET labels
labeller <- label_flqs(x[[1]])
# ASSEMBLE data
data <- lapply(x, as.data.frame, date=TRUE, drop=FALSE)
# SET stock names
stk <- rep.int(names(x), unlist(lapply(data, nrow)))
# RBIND dfs
data <- do.call(rbind, data)
rownames(data) <- NULL
# USE year or date for x axis
xvar <- sym(ifelse(length(unique(data$season)) == 1, "year", "date"))
# ADD stock names
data <- transform(data, stock=factor(stk, levels=names(x)))
# PLOT using geom_flquantiles
p <- ggplot(data, aes(x=!!xvar, y=data, fill=stock, colour=stock)) +
facet_grid(qname~., labeller=labeller, scales="free_y") +
xlab("") + ylab("") +
# SET limits to include 0
expand_limits(y=0) +
# SET legend with no title
theme(legend.title = element_blank()) +
# and only with lines and no title
guides(fill = "none")
if(length(probs) == 5) {
# outer quantile
p <- p + geom_flquantiles(probs=probs[c(1, 5)], alpha=alpha[1],
colour="white") +
# inner quantile
geom_flquantiles(probs=probs[c(2, 4)], alpha=alpha[2],
colour="white") +
# median
geom_flquantiles(probs=probs[3], alpha=1)
} else if (length(probs) == 3) {
p <- p + geom_flquantiles(probs=probs[c(1, 3)], alpha=alpha[2],
colour="white") +
geom_flquantiles(probs=probs[2], alpha=1)
} else if (length(probs) == 1) {
p <- p + geom_flquantiles(probs=probs, alpha=1)
} else {
stop("probs can only be of length 1, 3 or 5")
}
# PLOT iter worms
if(is.numeric(worm)) {
idata <- p$data[p$data$iter %in% worm,]
p <- p + geom_line(data=idata, aes(x=!!xvar, y=data, colour=iter))
}
return(p)
}
# }}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.