knitr::opts_chunk$set(echo = FALSE, cache = TRUE)
library(meteR) library(socorro) library(hawaiiMETE) grun <- grun[grun$Stage == 'adult', ] # METE objects and state vars meteByST <- plyr::dlply(grun[grun$trophic %in% c('P', 'H', 'D'), ], c('Site', 'trophic'), function(x) { meteESF(x$SpeciesCode, x$Abundance, x$IND_BIOM^0.75, minE = min(grun$IND_BIOM)^0.75) }) stateVarByST <- lapply(meteByST, function(x) c(x$state.var, nmax = max(tapply(x$data$n, x$data$s, sum)), emax = max(x$data$e))) stateVarByST <- as.data.frame(do.call(rbind, stateVarByST))
layout(matrix(1:15, nrow = 3)) par(mar = rep(0.3, 4), oma = c(3.5, 3.5, 2, 2), mgp = c(2, 0.75, 0)) for(i in grunSite$code) { for(j in c('P', 'H', 'D')) { thisMETE <- meteByST[[paste(i, j, sep = '.')]] plot(sad(thisMETE), ptype = 'rad', log = 'y', add.legend = TRUE, xlim = c(1, max(stateVarByST$S0)), ylim = c(1, max(stateVarByST$nmax)), xaxt = 'n', yaxt = 'n') if(i == 'KA') mtext(j, side = 4, line = 1) if(j == 'P') mtext(i, side = 3, line = 1) if(i == 'VO') logAxis(2, expLab = TRUE) if(j == 'D') axis(1) } } mtext('Rank', side = 1, outer = TRUE, line = 2) mtext('Abundance', side = 2, outer = TRUE, line = 2)
layout(matrix(1:15, nrow = 3)) par(mar = rep(0.3, 4), oma = c(3.5, 3.5, 2, 2), mgp = c(2, 0.75, 0)) for(i in grunSite$code) { for(j in c('P', 'H', 'D')) { thisMETE <- meteByST[[paste(i, j, sep = '.')]] plot(ipd(thisMETE), ptype = 'rad', log = 'y', add.legend = FALSE, xlim = c(1, max(stateVarByST$N0)), ylim = c(1, max(stateVarByST$emax)), xaxt = 'n', yaxt = 'n') if(i == 'KA') mtext(j, side = 4, line = 1) if(j == 'P') mtext(i, side = 3, line = 1) if(i == 'VO') logAxis(2, expLab = TRUE) if(j == 'D') axis(1) } } mtext('Rank', side = 1, outer = TRUE, line = 2) mtext('Metabolic rate', side = 2, outer = TRUE, line = 2)
# get bootstrapped MSE for SAD sadBootByST <- lapply(meteByST, function(x) { b <- bootMSE.meteDist(sad(x), type = 'cumulative', relative = FALSE, log = TRUE) o <- c(mse(sad(x), type = 'cumulative', relative = FALSE, log = TRUE), mean(b), quantile(b, probs = c(0.025, 0.975))) names(o) <- c('obs', 'mean', 'ciLo', 'ciHi') return(o) }) sadBootByST <- as.data.frame(do.call(rbind, sadBootByST)) sadBootByST <- cbind(site = gsub('\\..*', '', rownames(sadBootByST)), troph = gsub('.*\\.', '', rownames(sadBootByST)), sadBootByST) rownames(sadBootByST) <- NULL layout(matrix(1:3, nrow = 1)) par(mar = c(0.1, 3, 0, 0), oma = c(3, 0.5, 0.5, 0.5)) for(i in c('P', 'H', 'D')) { dat <- sadBootByST[sadBootByST$troph == i, ] plot(grunSite$age[match(grunSite$code, dat$site)], dat$mean, ylim = range(sadBootByST[, c('ciLo', 'ciHi')]), panel.first = segments(x0 = grunSite$age[match(grunSite$code, dat$site)], y0 = dat$ciLo, y1 = dat$ciHi), pch = 21, bg = 'white', log = 'x') }
sadZByST <- lapply(meteByST, function(x) { z <- logLikZ(sad(x))$z o <- c(logLik(sad(x)), z) names(o) <- c('obs', 'z') return(o) }) sadZByST <- as.data.frame(do.call(rbind, sadZByST)) sadZByST <- cbind(site = gsub('\\..*', '', rownames(sadZByST)), troph = gsub('.*\\.', '', rownames(sadZByST)), sadZByST) rownames(sadZByST) <- NULL layout(matrix(1:3, nrow = 1)) par(mar = c(3, 0.5, 0, 0), oma = c(0.5, 3, 0.5, 0.5), mgp = c(2, 0.75, 0)) for(i in c('P', 'H', 'D')) { plot(grunSite$age[match(grunSite$code, sadZByST$site[sadZByST$troph == i])], sadZByST$z[sadZByST$troph == i], log = 'x', xaxt = 'n', yaxt = 'n', ylim = range(sadZByST$z)) logAxis(1, expLab = TRUE) if(i == 'P') axis(2) } mtext('Age', side = 1, outer = TRUE, line = 0) mtext('log Likelihood z^2', side = 2, outer = TRUE, line = 2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.