R/gen_quantile_interval.R
gen_quantile_interval.RdReturn a Poisson-quantile-based prediction interval for qPCR values given Markov Chain Monte Carlo samples for the estimated concentrations.
gen_quantile_interval(
mu_quantiles,
mu_samps,
alpha = 0.05,
type = "credible_quantiles",
div_num = 1
)the (alpha/2, 1 - alpha/2) quantiles from the MCMC sampling distribution for the true absolute concentration \(\mu\); only supply if type = "credible_quantiles" (a matrix of dimension N (the sample size) x q (the number of taxa) by 2).
the estimated concentrations [an array with dimension (number of MCMC samples) by N by q]; only supply if type = "sample_quantiles".
the desired level (defaults to 0.05, corresponding to an interval using the 2.5% and 97.5% quantiles)
the type of intervals desired, either "credible_quantiles" or "sample_quantiles" (please see Details for more information on the difference between these two types).
the number to multiply by.
A (1 - \(\alpha\))x100% Poisson-quantile-based prediction interval for each qPCR
# load the package, read in example data
library("paramedic")
data(example_16S_data)
data(example_qPCR_data)
# run paramedic (with an extremely small number of iterations, for illustration only)
# on only the first 10 taxa
mod <- run_paramedic(W = example_16S_data[, 1:10], V = example_qPCR_data,
n_iter = 30, n_burnin = 25,
n_chains = 1, stan_seed = 4747)
#>
#> SAMPLING FOR MODEL 'variable_efficiency' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000129 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.29 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1: three stages of adaptation as currently configured.
#> Chain 1: Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1: the given number of warmup iterations:
#> Chain 1: init_buffer = 3
#> Chain 1: adapt_window = 20
#> Chain 1: term_buffer = 2
#> Chain 1:
#> Chain 1: Iteration: 1 / 30 [ 3%] (Warmup)
#> Chain 1: Iteration: 3 / 30 [ 10%] (Warmup)
#> Chain 1: Iteration: 6 / 30 [ 20%] (Warmup)
#> Chain 1: Iteration: 9 / 30 [ 30%] (Warmup)
#> Chain 1: Iteration: 12 / 30 [ 40%] (Warmup)
#> Chain 1: Iteration: 15 / 30 [ 50%] (Warmup)
#> Chain 1: Iteration: 18 / 30 [ 60%] (Warmup)
#> Chain 1: Iteration: 21 / 30 [ 70%] (Warmup)
#> Chain 1: Iteration: 24 / 30 [ 80%] (Warmup)
#> Chain 1: Iteration: 26 / 30 [ 86%] (Sampling)
#> Chain 1: Iteration: 28 / 30 [ 93%] (Sampling)
#> Chain 1: Iteration: 30 / 30 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.016 seconds (Warm-up)
#> Chain 1: 0.002 seconds (Sampling)
#> Chain 1: 0.018 seconds (Total)
#> Chain 1:
#> Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
#> https://mc-stan.org/misc/warnings.html#bfmi-low
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 2.5, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
# get model summary
mod_summ <- rstan::summary(mod, probs = c(0.025, 0.975))$summary
#> Error in rstan::summary(mod, probs = c(0.025, 0.975))$summary: $ operator is invalid for atomic vectors
# get samples
mod_samps <- rstan::extract(mod$stan_fit)
# extract relevant summaries
summs <- extract_posterior_summaries(stan_mod = mod_summ, stan_samps = mod_samps,
taxa_of_interest = 1:3,
mult_num = 1, level = 0.95, interval_type = "quantile")
#> Error in eval(expr, envir, enclos): object 'mod_summ' not found