Merge pull request #2965 from nanu1605/meta_analysis_global_variables

fixed some global variables in 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/NAMESPACE

@@ -2,7 +2,6 @@
 
 export(approx.posterior)
 export(jagify)
-export(p.point.in.prior)
 export(pecan.ma)
 export(pecan.ma.summary)
 export(rename_jags_columns)

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, stats::dbeta(x, fit$estimate[1], fit$estimate[2]), lwd = 2, type = "l")
+        graphics::lines(x, stats::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)
       }
@@ -93,7 +93,7 @@ approx.posterior <- function(trait.mcmc, priors, trait.data = NULL, outdir = NUL
       dist.names <- dist.names[!failfit.bool]
 
       fparm <- lapply(fit, function(x) { as.numeric(x$estimate) })
-      fAIC <- lapply(fit, AIC) 
+      fAIC <- lapply(fit, stats::AIC)
 
       bestfit <- which.min(fAIC)
       posteriors[ptrait, "distn"] <- dist.names[bestfit]
@@ -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/jagify.R

@@ -14,8 +14,8 @@
 ##' @title Prepare trait data for JAGS meta-analysis
 ##' @param result input trait data
 ##' @param use_ghs (Logical) If `FALSE`, exclude all greenhouse data. If `TRUE`, use all data, including greenhouse data.
-##' @return result transformed to meet requirements of PEcAn meta-analysis model
 ##' @export
+##' @return result transformed to meet requirements of PEcAn meta-analysis model
 ##' @author David LeBauer
 jagify <- function(result, use_ghs = TRUE) {
 

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,8 @@
 
 ##' 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.
+##' Inspired by the \code{R2WinBUGS::write.model} by Jouni Kerman and Uwe Ligges.
 ##' @name write.ma.model
 ##' @title write.ma.model
 ##' @param modelfile model template file (ma.model.template.R)
@@ -25,7 +26,7 @@
 ##' @param tauA parameter a for gamma prior on precision 
 ##' @param tauB parameter b for gamma prior on precision
 ##' @return Nothing, but as a side effect, the model is written
-##' @author David LeBauer and Mike Dietze, based on original work on the \code{write.model} function in the \code{R2WinBUGS} package by Jouni Kerman and Uwe Ligges.
+##' @author David LeBauer and Mike Dietze.
 write.ma.model <- function(modelfile, outfile, reg.model, pr.dist, pr.param.a, pr.param.b, n, 
                            trt.n, site.n, ghs.n, tauA, tauB) {