sensAnalysis_2: Second sensitivity analysis In jafernandez01/cndcR: Critical Nitrogen Dilution Curves in R

 sensAnalysis_2 R Documentation

Second sensitivity analysis

Description

This function returns results from the 2nd sensitivity analysis in Fernandez et al.(2022): number of N rates. Output is a tibble with posterior expectations and credibility intervals of a (A1) and b (A2) parameters of the CNDC.

Details

See 'Examples' for the source code to obtain the given results. Briefly, a bootstrap procedure in which studies with r (where r = 2, 3, 4, 5, and 6) number of N rates 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.

Note

Parallel computation is recommended if running the example code.

Examples

```## Not run:

data("Data")

set.seed(500)
dataSens_2 <- NULL
for (i in seq(1, 5, 1)) {
dataSens_2[[i]] <- replicate(100, Data  %>%
group_by(Site_year)  %>%
dplyr::filter(length(unique(Nrates)) > 5)  %>%
nest()  %>%
mutate(sample = purrr::map(data, ~ dplyr::filter(.x,
Nrates <= 30 | Nrates %in% sample(unique(dplyr::filter(.x,
Nrates > 30)\$Nrates), 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 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_2 <- NULL
bootSamp_2 <- foreach::foreach(i = seq(1, 5, 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_2 <- NULL
for (d in seq(1, 100, 1)) {
Date <- as.numeric(as.factor(paste(dataSens_2[[i]][[d]]\$Site_year, "_",
dataSens_2[[i]][[d]]\$Samp, sep = "")))

dataList <- list(
"W" = dataSens_2[[i]][[d]]\$W,
"N" = dataSens_2[[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_2 <- rjags::jags.model(m,
data = dataList,
n.chains = nChains,
)
close(m)

if (burnInSteps > 0) {
update(jagsModel_2, n.iter = burnInSteps)
}
codaSamples_2 <- rjags::coda.samples(jagsModel_2,
variable.names = parameters,
n.iter = nIter,
thin = thinSteps
)
mcmcChain_2[[d]] <- as.matrix(codaSamples_2)
}

bootSamp_2[[i]] <- mcmcChain_2
}

# 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_2 <- foreach::foreach(i = seq(1, 5, 1)) %do% {

purrr::map2_dfr(
bootSamp_2[[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_2 <- bind_rows(fdataSens_2)  %>%
dplyr::group_by(Parameter, Method)  %>%
dplyr::summarise(lowCI = quantile(Value, 0.025), uppCI = quantile(Value, 0.975),
mean = mean(Value))

## End(Not run)
```

jafernandez01/cndcR documentation built on Oct. 26, 2022, 5:24 a.m.