View source: R/sensAnalysis_3.R

sensAnalysis_3 | R Documentation |

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.

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.

Parallel computation is recommended if running the example code.

## Not run: 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)) ## End(Not run)

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.