R/sensAnalysis_3.R

In jafernandez01/cndcR: Critical Nitrogen Dilution Curves in R

#' @name sensAnalysis_3
#' @title
#' Third sensitivity analysis
#' @details
#' See 'Examples' for the source code to obtain the given results. Briefly, a bootstrap procedure
#' in which studies with s (where s = 1, 2, 3, 4, 5, and 6) number of samplings are sampled 100
#' times with replacement to create 100 new datasets. The Bayesian hierarchical model explained by
#' Makowski et al. (2020) was fitted to each bootstrapped sample.
#' @description
#' This function returns results from the 3rd sensitivity analysis in Fernandez et al.(2022):
#' number of sampling times. Output is a tibble with posterior expectations and credibility
#' intervals of a (A1) and b (A2) parameters of the CNDC.
#'
#' @note Parallel computation is recommended if running the example code.
#'
#' @examples
#' \dontrun{
#'
#' data("Data")
#'
#' set.seed(500)
#' dataSens_3 <- NULL
#' for (i in seq(1, 6, 1)) {
#'   dataSens_3[[i]] <- replicate(100, Data  %>%
#'                                  group_by(Site_year)  %>%
#'                                  filter(length(unique(Samp)) > 5)  %>%
#'                                  nest()  %>%
#'                                  mutate(sample = purrr::map(data, ~ filter(.x,
#'                                  Samp %in% sample(unique(.x$Samp), i, replace = F))))  %>%
#'                                  unnest(sample), simplify = F)
#' }
#'
#' # Parameters and JAGS settings are defined for the MCMC (Markov Chain Monte Carlo) procedure to
#' # model the power function of CNDC across sampling dates. Weakly-informative priors were defined
#' # following Makowski et al. (2020) and Ciampitti et al. (2021). To facilitate computation of the
#' # sensitivity analyses, the algorithm was run with three chains of 20,000 iterations each
#' # (10,000 discarded as tuning and burn-in periods). The statistical model was fitted using the
#' # rjags:: library.
#'
#'     set.seed(500)
#'     bootSamp_3 <- NULL
#'     bootSamp_3 <- foreach(i = seq(1, 6, 1), .verbose = T) %do% {
#'
#'       # Specify parameters and JAGS settings
#'       parameters <- c("A1", "A2", "Bmax", "S", "Nc")
#'
#'       # JAGS settings
#'       adaptSteps <- 5000 # number of steps to "tune" the samplers
#'       burnInSteps <- 5000 # number of steps to "burn-in" the samplers
#'       nChains <- 3 # number of chains to run
#'       thinSteps <- 10 # number of steps to "thin" (keep every 10 steps)
#'       nIter <- 10000 # steps per chain
#'
#'       mcmcChain_3 <- NULL
#'       for (d in seq(1, 100, 1)) {
#'         Date <- as.numeric(as.factor(paste(dataSens_3[[i]][[d]]$Site_year, "_",
#'         dataSens_3[[i]][[d]]$Samp, sep = "")))
#'
#'         dataList <- list(
#'           "W" = dataSens_3[[i]][[d]]$W,
#'           "N" = dataSens_3[[i]][[d]]$Na,
#'           "Date" = Date,
#'           "Q" = length(Date),
#'           "K" = length(unique(Date))
#'         )
#'
#'         m <- textConnection("
#' model {
#'
#' 	for (i in 1:Q)
#' 	{
#' 		W[i]~dnorm(mu[i], tau_b)
#' 		N[i]~dnorm(Nc[Date[i]], tau_n)
#' 		mu[i]<-min(Bmax[Date[i]], Bmax[Date[i]]+S[Date[i]]*(N[i]-Nc[Date[i]]))
#' 	}
#'
#' 	for (j in 1:K)
#' 	{
#' 		Nc[j]=A1*Bmax[j]^(-A2)
#' 		Bmax[j]~dnorm(Mu_Bmax,Prec_Bmax)T(0,)
#' 		S[j]~dnorm(Mu_S,Prec_S)T(0,)
#' 			}
#'
#' 			#Weakly informative
#' 			Mu_Bmax~dnorm(6,0.1)
#' 			Mu_S~dnorm(0,0.1)
#' 			A1~dunif(2,6)
#' 			A2~dunif(0,0.7)
#'
#' 			Prec_Bmax~dgamma(0.001,0.001)
#' 			Prec_S~dgamma(0.001,0.001)
#' 			tau_b~dgamma(0.001,0.001)
#' 			tau_n~dgamma(0.001,0.001)
#'
#' }
#' ")
#'         set.seed(500)
#'         jagsModel_3 <- rjags::jags.model(m,
#'                                          data = dataList,
#'                                          n.chains = nChains,
#'                                          n.adapt = adaptSteps
#'         )
#'         close(m)
#'
#'         if (burnInSteps > 0) {
#'           update(jagsModel_3, n.iter = burnInSteps)
#'         }
#'         codaSamples_3 <- rjags::coda.samples(jagsModel_3,
#'                                              variable.names = parameters,
#'                                              n.iter = nIter,
#'                                              thin = thinSteps
#'         )
#'         mcmcChain_3[[d]] <- as.matrix(codaSamples_3)
#'       }
#'
#'
#'     bootSamp_3[[i]] <- mcmcChain_3
#'     }
#'
#'
#' # Posterior probability distributions of *a* and *b* parameters are extracted for each
#' # bootstrapped sample (i.e., extracting a total of 100 probability distributions). These
#' # distributions are then being averaged to obtain posterior medians and credibility intervals
#' # from all the samples.
#'
#' fdataSens_3 <- foreach(i = seq(1, 6, 1)) %do% {
#'
#'   purrr::map2_dfr(
#'     bootSamp_3[[i]], seq(1, 100, 1),
#'     ~ .x  %>%
#'       as.data.frame()  %>%
#'       dplyr::mutate(Imp = .y)  %>%
#'       tidyr::pivot_longer(-Imp, names_to = "Parameter", values_to = "Value")
#'   )  %>%
#'     dplyr::mutate(
#'       Method = paste("Boot", i, sep = ""),
#'       Samp = stringr::str_replace(Parameter, ".*\\[(\\d{1,2})\\]$", "\\1"),
#'       Parameter = stringr::str_remove(Parameter, "\\[\\d{1,2}\\]$")
#'     )  %>%
#'     dplyr::filter(Parameter %in% c("A1", "A2"))
#' }
#'
#' fdataSens_3 <- bind_rows(fdataSens_3)  %>%
#'   dplyr::group_by(Parameter, Method)  %>%
#'   dplyr::summarise(lowCI = quantile(Value, 0.025), uppCI = quantile(Value, 0.975),
#'   mean = mean(Value))
#'
#' }
NULL

sensAnalysis_3 <- function() {
  return(eval(parse(text = "cndcR:::fdataSens_3")))
}