knitr::opts_chunk$set(echo = FALSE, cache = TRUE)

Introduction

Methods

Results

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)

Discussion

Conclusion

References



ajrominger/hawaiiMETE documentation built on May 21, 2019, 1:42 p.m.