#' Plot mean-variance relationship
#'
#' Creates a scatter plot of the time series log(variance) vs.
#' log(mean). Shows a model fit to the mean-variance data and an
#' extrapolation to the size of the metapopulation. The linear version
#' of this model is referred to as Taylor's power law.
#'
#' @param x A matrix or dataframe of abundance or biomass data. The
#' columns should represent different subpopulations or species. The
#' rows should represent the values through time.
#' @param show A vector of character objects indicating which
#' mean-variance models to show.
#' @param col Colour for the mean-variance model fit. A vector of
#' length 3 with the three values corresponding to \code{linear},
#' \code{quadratic},
#' \code{robust}.
#' @param lty Line type for the mean-variance model fit. A vector of
#' length 3 with the three values corresponding to \code{linear},
#' \code{quadratic}, \code{robust}.
#' @param pch_sa Point type for the extrapolated
#' "single-asset" portfolio. A vector of length 3 with the three
#' values corresponding to \code{linear}, \code{quadratic},
#' \code{robust}.
#' @param ci Add a confidence interval around the model fit? Only
#' appears for the linear fit option.
#' @param pch_subpops Point type for the subpopulations.
#' @param pch_port Point type for the portfolio.
#' @param add_z Logical. Add Taylor's power law z value (based on a
#' linear model)?
#' @param xlab Label for x-axis.
#' @param ylab Label for y-axis.
#' @param ... Other values to pass to \code{plot}.
#' @references
#' Anderson, S.C., A.B. Cooper, N.K. Dulvy. 2013. Ecological prophets:
#' Quantifying metapopulation portfolio effects. Methods in Ecology and
#' Evolution. In Press.
#'
#' Doak, D., D. Bigger, E. Harding, M. Marvier, R. O'Malley, and D.
#' Thomson. 1998. The Statistical Inevitability of Stability-Diversity
#' Relationships in Community Ecology. Amer. Nat. 151:264-276.
#'
#' Tilman, D., C. Lehman, and C. Bristow. 1998. Diversity-Stability
#' Relationships: Statistical Inevitability or Ecological Consequence?
#' Amer. Nat. 151:277-282.
#'
#' Tilman, D. 1999. The Ecological Consequences of Changes in
#' Biodiversity: A Search for General Principles. Ecology
#' 80:1455-1474.
#'
#' Taylor, L. 1961. Aggregation, Variance and the Mean. Nature
#' 189:732-735. doi: 10.1038/189732a0.
#'
#' Taylor, L., I. Woiwod, and J. Perry. 1978. The Density-Dependence
#' of Spatial Behaviour and the Rarity of Randomness. J. Anim. Ecol.
#' 47:383-406.
#' @export
#' @examples
#' data(pinkbr)
#' par(mfrow = c(1,3))
#' plot_mv(pinkbr[,-1], show = "linear")
#' mtext("Linear")
#' plot_mv(pinkbr[,-1], show = "quadratic", add_z = FALSE)
#' mtext("Quadratic")
#' plot_mv(pinkbr[,-1], show = "robust", add_z = FALSE)
#' mtext("Robust linear")
#' @import robustbase
# @importMethodsFrom robustbase predict.lmrob
# @importClassesFrom robustbase lmrob
plot_mv <- function(x, show = c("linear", "quadratic", "robust"), col
= c("#D95F02", "#1B9E77", "#E7298A"), lty = c(1, 1, 1),
pch_sa = c(1, 5, 6), ci = FALSE, pch_subpops = 21,
pch_port = 4, add_z = TRUE, xlab = "log(mean)", ylab =
"log(variance)", ...) {
## get mean and variance of portfolio and assets:
x.long <- reshape::melt.data.frame(x, id.vars = NULL)
overall.d <- apply(x[,-1], 1, sum)
mv <- plyr::ddply(x.long, "variable", function(x) {data.frame(m =
mean(x$value, na.rm = TRUE), v = var(x$value, na.rm = TRUE))})
overall.mean <- mean(overall.d, na.rm = TRUE)
portfolio.var <- var(overall.d, na.rm = TRUE)
m.t <- lm(log(v) ~ log(m), data = mv)
overall.variance <- exp(as.numeric(predict(m.t, newdata =
data.frame(m = overall.mean))))
m.t.quad <- nls(log(v) ~ B0 + B1 * log(m) + B2 * I(log(m) ^ 2), data
= mv, start = list(B0 = 0, B1 = 2, B2 = 0), lower = list(B0 =
-1e9, B1 = 0, B2 = 0), algorithm = "port")
overall.variance.quad <- exp(as.numeric(predict(m.t.quad, newdata =
data.frame(m = overall.mean))))
d1p <- seq(min(mv$m), max(mv$m), length.out = 200)
d2p <- seq(max(mv$m), overall.mean, length.out = 200)
## set up plot:
with(mv, plot((m), (v), xlim = c((min(m)), (overall.mean)*1.15),
ylim = c(min((v)), max(overall.variance, portfolio.var,
overall.variance.quad)*1.25), pch = pch_subpops, log = "xy", xlab =
xlab, ylab =ylab, bg = "#00000020", col = "#00000070", ...))
points((overall.mean), (portfolio.var), col = "black", pch = pch_port, lwd
= 1.2, cex = 1.6)
## now add fits and extrapolations as requested:
if("linear" %in% show) {
m.t.p <- predict(m.t, newdata = data.frame(m = seq(min(mv$m),
max(mv$m), length.out = 2)))
m.t.p2 <- predict(m.t, newdata = data.frame(m = seq(max(mv$m),
overall.mean, length.out = 2)))
p1 <- predict(m.t, newdata = data.frame(m = d1p), se = TRUE)
p2 <- predict(m.t, newdata = data.frame(m = d2p), se = TRUE)
points((overall.mean), (overall.variance), col = col[1], pch = pch_sa[1],
lwd = 1.5, cex = 1.4)
segments((min(mv$m)), exp(m.t.p[1]), (max(mv$m)), exp(m.t.p[2]),
col = col[1], lty = lty[1])
segments((max(mv$m)), exp(m.t.p2[1]), (overall.mean),
exp(m.t.p2[2]), lty = 2, col = col[1])
if(ci) {
polygon(c(d1p, rev(d1p)), c(exp(p1$fit + 1.96*p1$se.fit),
exp(rev(p1$fit - 1.96*p1$se.fit))), border = FALSE, col =
"#00000020")
polygon(c(d2p, rev(d2p)), c(exp(p2$fit + 1.96*p2$se.fit),
exp(rev(p2$fit - 1.96*p2$se.fit))), border = FALSE, col =
"#00000010")
}
}
if("quadratic" %in% show){
p.quad.1 <- predict(m.t.quad, newdata = data.frame(m = d1p), se = FALSE)
p.quad.2 <- predict(m.t.quad, newdata = data.frame(m = d2p), se = FALSE)
points((overall.mean), (overall.variance.quad), col =col[2], pch =
pch_sa[2], lwd = 1.5, cex = 1.4)
lines(d1p, exp(as.numeric(p.quad.1)), col = col[2], lwd = lty[2])
lines(d2p, exp(as.numeric(p.quad.2)), col = col[2], lty = 2, lwd = 1.5)
}
if("robust" %in% show) {
m.t.rob <- lmrob(log(v) ~ log(m), data = mv)
overall.variance.rob <- exp(as.numeric(predict(m.t.rob, newdata =
data.frame(m = overall.mean))))
p.rob.1 <- predict(m.t.rob, newdata = data.frame(m = d1p), se = FALSE)
p.rob.2 <- predict(m.t.rob, newdata = data.frame(m = d2p), se = FALSE)
points((overall.mean), (overall.variance.rob), col = col[3], pch =
pch_sa[3], lwd = 1.5, cex = 1.1)
lines(d1p, exp(p.rob.1), col = col[3], lty = lty[3])
lines(d2p, exp(p.rob.2), col = col[3], lty = 2)
}
if(add_z) {
mtext(paste("z =", formatC(round(coef(m.t)[2], 2), digits = 1,
format = "f")), side = 1, adj = 0.9, line = -1.2)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.