Fix some warnings and notes from PEcAn.MA package.

docker/depends/pecan.depends.R

@@ -45,6 +45,7 @@ wanted <- c(
 'geonames',
 'getPass',
 'ggmap',
+'ggmcmc',
 'ggplot2',
 'ggrepel',
 'glue',

modules/meta.analysis/DESCRIPTION

@@ -25,21 +25,25 @@ Description: The Predictive Ecosystem Carbon Analyzer (PEcAn) is a scientific
     streamline the interaction between data and models, and to improve the
     efficacy of scientific investigation. The PEcAn.MA package contains
     the functions used in the Bayesian meta-analysis of trait data.
-Depends:
+Imports:
+    coda (>= 0.18),
     XML,
     lattice,
-    MASS,
     PEcAn.utils,
-    PEcAn.DB
-Imports:
-    coda (>= 0.18),
+    PEcAn.DB,
     PEcAn.logger,
+    MASS,
     PEcAn.settings,
     rjags
 Suggests:
+    ggmcmc,
+    ggplot2,
+    knitr,
+    rmarkdown,
     testthat (>= 1.0.2)
 SystemRequirements: JAGS
 License: BSD_3_clause + file LICENSE
+VignetteBuilder: knitr
 Copyright: Authors
 LazyData: FALSE
 Encoding: UTF-8

modules/meta.analysis/R/approx.posterior.R

@@ -34,7 +34,7 @@ approx.posterior <- function(trait.mcmc, priors, trait.data = NULL, outdir = NUL
   posteriors <- priors
   do.plot <- !is.null(outdir)
   if (do.plot == TRUE) {
-    pdf(file.path(outdir, paste("posteriors", filename.flag, ".pdf", sep = "")))
+    grDevices::pdf(file.path(outdir, paste("posteriors", filename.flag, ".pdf", sep = "")))
   }
 
   ## loop over traits
@@ -55,21 +55,21 @@ approx.posterior <- function(trait.mcmc, priors, trait.data = NULL, outdir = NUL
     zerobound <- c("exp", "gamma", "lnorm", "weibull")
     if (pdist %in% "beta") {
       m   <- mean(dat)
-      v   <- var(dat)
+      v   <- stats::var(dat)
       k   <- (1 - m)/m
       a   <- (k / ((1 + k) ^ 2 * v) - 1) / (1 + k)
       b   <- a * k
       fit <- try(suppressWarnings(MASS::fitdistr(dat, "beta", list(shape1 = a, shape2 = b))), silent = TRUE)
 
       if (do.plot) {
         x <- seq(0, 1, length = 1000)
-        plot(density(dat), col = 2, lwd = 2, main = trait)
+        plot(stats::density(dat), col = 2, lwd = 2, main = trait)
         if (!is.null(trait.data)) {
-          rug(trait.data[[trait]]$Y, lwd = 2, col = "purple")
+          graphics::rug(trait.data[[trait]]$Y, lwd = 2, col = "purple")
         }
-        lines(x, dbeta(x, fit$estimate[1], fit$estimate[2]), lwd = 2, type = "l")
-        lines(x, dbeta(x, pparm[1], pparm[2]), lwd = 3, type = "l", col = 3)
-        legend("topleft", 
+        graphics::lines(x, dbeta(x, fit$estimate[1], fit$estimate[2]), lwd = 2, type = "l")
+        graphics::lines(x, dbeta(x, pparm[1], pparm[2]), lwd = 3, type = "l", col = 3)
+        graphics::legend("topleft", 
                legend = c("data", "prior", "post", "approx"), 
                col = c("purple", 3, 2, 1), lwd = 2)
       }
@@ -111,15 +111,15 @@ approx.posterior <- function(trait.mcmc, priors, trait.data = NULL, outdir = NUL
       ## default: NORMAL
       posteriors[trait, "distn"] <- "norm"
       posteriors[trait, "parama"] <- mean(dat)
-      posteriors[trait, "paramb"] <- sd(dat)
+      posteriors[trait, "paramb"] <- stats::sd(dat)
       if (do.plot) {
         .dens_plot(posteriors, priors, ptrait, dat, trait, trait.data)
       }
     }
   }  ## end trait loop
 
   if (do.plot) {
-    dev.off()
+    grDevices::dev.off()
   }
 
   return(posteriors)
@@ -159,13 +159,13 @@ approx.posterior <- function(trait.mcmc, priors, trait.data = NULL, outdir = NUL
     rng <- range(trait.data[[trait]]$Y)
   }
 
-  plot(density(dat), col = 2, lwd = 2, main = trait, xlim = rng)
+  plot(stats::density(dat), col = 2, lwd = 2, main = trait, xlim = rng)
   if (!is.null(trait.data)) {
-    rug(trait.data[[trait]]$Y, lwd = 2, col = "purple")
+    graphics::rug(trait.data[[trait]]$Y, lwd = 2, col = "purple")
   }
-  lines(x, f(x), lwd = 2, type = "l")
-  lines(x, fp(x), lwd = 3, type = "l", col = 3)
-  legend("topleft", 
+  graphics::lines(x, f(x), lwd = 2, type = "l")
+  graphics::lines(x, fp(x), lwd = 3, type = "l", col = 3)
+  graphics::legend("topleft", 
          legend = c("data", "prior", "post", "approx"), 
          col = c("purple", 3, 2, 1), lwd = 2)
 }

modules/meta.analysis/R/meta.analysis.R

@@ -27,6 +27,9 @@
 ##' @param prior.distns `data.frame` of prior distributions generated
 ##'   by call to [PEcAn.DB::query.priors()]
 ##' @param taupriors priors on variance parameters, can be scaled as needed with data mean
+##' @param data data frame generated by jagify function
+##'   with indexed values for greenhouse, treatment, and site (ghs, trt, site)
+##'   as well as Y, SE, and n for each observation or summary statistic.
 ##' @param j.iter number of MCMC samples
 ##' @param outdir output directory
 ##' @param random use random effects, `FALSE` by default
@@ -124,10 +127,10 @@ pecan.ma <- function(trait.data, prior.distns,
                        "\nmean:", signif(mean(data$Y, na.rm = TRUE), 3),
                        "\nn:", length(data$Y)))
       writeLines("stem plot of data points")
-      writeLines(paste(stem(data$Y)))
+      writeLines(paste(graphics::stem(data$Y)))
       if (any(!is.na(data$obs.prec)) && all(!is.infinite(data$obs.prec))) {
         writeLines("stem plot of obs.prec:")
-        writeLines(paste(stem(data[["obs.prec"]] ^ 2)))
+        writeLines(paste(graphics::stem(data[["obs.prec"]] ^ 2)))
       } else {
         writeLines(paste("no estimates of SD for", trait.name))
       }
@@ -152,7 +155,7 @@ pecan.ma <- function(trait.data, prior.distns,
       print(summary(jags.out))
     }
 
-    jags.out.trunc <- window(jags.out, start = j.iter / 2)
+    jags.out.trunc <- stats::window(jags.out, start = j.iter / 2)
 
     mcmc.object[[trait.name]] <- jags.out.trunc
   }

modules/meta.analysis/R/meta.analysis.summary.R

@@ -27,7 +27,6 @@
 ##' }
 ##' @author David LeBauer, Shawn Serbin
 pecan.ma.summary <- function(mcmc.object, pft, outdir, threshold = 1.2, gg = FALSE) {
-
   fail <- rep(FALSE, length(mcmc.object))
   names(fail) <- names(mcmc.object)
   not.converged <- data.frame()
@@ -36,13 +35,18 @@ pecan.ma.summary <- function(mcmc.object, pft, outdir, threshold = 1.2, gg = FAL
   for (trait in names(mcmc.object)) {
 
     if (gg) {
-      gg <- require(ggmcmc)
+      if (!requireNamespace("ggmcmc", quietly = TRUE)) {
+        PEcAn.logger::logger.severe(
+          "Can't find package 'ggmcmc',",
+          "needed by `PEcAn.MA::meta.analysis.summary()` when `gg = TRUE`.",
+          "Please install it and try again.")
+      }
     }
     ## new diagnostic plots. But very slow & !any(grepl('^gg', dir(outdir)))){
     if (gg) {
-      if (is.mcmc.list(mcmc.object[[trait]])) {
-        theme_set(theme_bw())
-        ggmcmc(ggs(mcmc.object[[trait]]), 
+      if (coda::is.mcmc.list(mcmc.object[[trait]])) {
+        ggplot2::theme_set(ggplot2::theme_bw())
+        ggmcmc::ggmcmc(ggmcmc::ggs(mcmc.object[[trait]]), 
                plot = c("ggs_density", "ggs_traceplot", "ggs_autocorrelation", "ggs_Rhat", "ggs_geweke"), 
                file.path(outdir, paste0("gg.ma.summaryplots.", trait, ".pdf")))
       }
@@ -54,23 +58,23 @@ pecan.ma.summary <- function(mcmc.object, pft, outdir, threshold = 1.2, gg = FAL
     maparms  <- .maparms[c(which(.maparms %in% .parms), which(!.maparms %in% .parms))]
 
     ## plots for mcmc diagnosis
-    pdf(file.path(outdir, paste0("ma.summaryplots.", trait, ".pdf")))
+    grDevices::pdf(file.path(outdir, paste0("ma.summaryplots.", trait, ".pdf")))
 
     for (i in maparms) {
       plot(mcmc.object[[trait]][, i],
            trace = FALSE,
            density = TRUE,
            main = paste("summary plots of", i, "for", pft, trait))
-      box(lwd = 2)
+      graphics::box(lwd = 2)
       plot(mcmc.object[[trait]][, i], density = FALSE)
-      box(lwd = 2)
+      graphics::box(lwd = 2)
       coda::autocorr.plot(mcmc.object[[trait]][, i][1], xlim = c(1, 50))
-      box(lwd = 2)
+      graphics::box(lwd = 2)
     }
     lattice::xyplot(mcmc.object[[trait]])
     lattice::densityplot(mcmc.object[[trait]])
     coda::acfplot(mcmc.object[[trait]])
-    dev.off()
+    grDevices::dev.off()
 
     ## G-R diagnostics to ensure convergence
     gd            <- coda::gelman.diag(mcmc.object[[trait]])

modules/meta.analysis/R/meta.analysis.write.model.R

@@ -9,7 +9,7 @@
 
 ##' Convert template ma.model.template.R to a JAGS model.
 ##'
-##'  Writes a meta-analysis model based on available data and prior specification. Adapted from \code{\link[R2WinBUGS]{write.model}}
+##'  Writes a meta-analysis model based on available data and prior specification. Adapted from \code{R2WinBUGS::write.model}
 ##' @name write.ma.model
 ##' @title write.ma.model
 ##' @param modelfile model template file (ma.model.template.R)

modules/meta.analysis/R/run.meta.analysis.R

@@ -94,7 +94,7 @@ run.meta.analysis.pft <- function(pft, iterations, random = TRUE, threshold = 1.
 
   ## Check that data is consistent with prior
   for (trait in names(jagged.data)) {
-    data.median <- median(jagged.data[[trait]]$Y)
+    data.median <- stats::median(jagged.data[[trait]]$Y)
     prior       <- prior.distns[trait, ]
     check       <- check_consistent(data.median, prior, trait, "data")
     if (is.na(check)) {
@@ -123,7 +123,7 @@ run.meta.analysis.pft <- function(pft, iterations, random = TRUE, threshold = 1.
 
   ### Check that meta-analysis posteriors are consistent with priors
   for (trait in names(trait.mcmc)) {
-    post.median <- median(as.matrix(trait.mcmc[[trait]][, "beta.o"]))
+    post.median <- stats::median(as.matrix(trait.mcmc[[trait]][, "beta.o"]))
     prior       <- prior.distns[trait, ]
     check <- check_consistent(post.median, prior, trait, "data")
     if (is.na(check)) {
@@ -157,7 +157,7 @@ run.meta.analysis.pft <- function(pft, iterations, random = TRUE, threshold = 1.
     }
     filename <- file.path(pathname, file)
     file.copy(file.path(pft$outdir, file), filename)
-    dbfile.insert(pathname, file, "Posterior", pft$posteriorid, dbcon)
+    PEcAn.DB::dbfile.insert(pathname, file, "Posterior", pft$posteriorid, dbcon)
   }
 } # run.meta.analysis.pft
 
@@ -185,8 +185,8 @@ run.meta.analysis.pft <- function(pft, iterations, random = TRUE, threshold = 1.
 ##' @author Shawn Serbin, David LeBauer
 run.meta.analysis <- function(pfts, iterations, random = TRUE, threshold = 1.2, dbfiles, database, use_ghs = TRUE) {
   # process all pfts
-  dbcon <- db.open(database)
-  on.exit(db.close(dbcon), add = TRUE)
+  dbcon <- PEcAn.DB::db.open(database)
+  on.exit(PEcAn.DB::db.close(dbcon), add = TRUE)
 
   result <- lapply(pfts, run.meta.analysis.pft, iterations = iterations, random = random, 
                    threshold = threshold, dbfiles = dbfiles, dbcon = dbcon, use_ghs = use_ghs)
@@ -235,15 +235,12 @@ runModule.run.meta.analysis <- function(settings) {
 ##'
 ##' used to compare data to prior, meta analysis posterior to prior
 ##' @title find quantile of point within prior distribution
-##' @param point 
+##' @param point quantile of given prior to return
 ##' @param prior list of distn, parama, paramb
 ##' @return result of p<distn>(point, parama, paramb)
 ##' @export p.point.in.prior
 ##' @author David LeBauer
 p.point.in.prior <- function(point, prior) {
-  # Why is this (below) called, and then never used?
-  prior.median <- do.call(paste0("q", prior$distn),
-                          list(0.5, prior$parama, prior$paramb))
   out <- do.call(paste0("p", prior$distn), 
                  list(point, prior$parama, prior$paramb))
   return(out)

modules/meta.analysis/R/single.MA.R

@@ -13,8 +13,10 @@
 ##' single.MA. This will allow power analysis  to run repeated MA outside of the
 ##' full loop over traits and PFTs.
 ##' @title Single MA
-##' @param data
 ##' @param j.chains number of chains in meta-analysis
+##' @param data data frame generated by jagify function
+##'   with indexed values for greenhouse, treatment, and site (ghs, trt, site)
+##'   as well as Y, SE, and n for each observation or summary statistic.
 ##' @param j.iter  number of mcmc samples
 ##' @param tauA prior on variance parameters
 ##' @param tauB  prior on variance parameters

modules/meta.analysis/vignettes/single.MA_demo.Rmd

@@ -1,3 +1,11 @@
+---
+title: "Example of Single.MA"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteIndexEntry{Example of Single.MA}
+  %\usepackage[utf8]{inputenc}
+---
 Meta-analysis: Example of Single.MA
 ========================================================
 Individual Meta-analysis