R/gen_wald_interval.R
gen_wald_interval.RdReturn a Wald-type prediction interval for qPCR values given the estimated concentrations.
gen_wald_interval(mu, sd, alpha = 0.05, truncate = TRUE)the estimated concentrations (a vector of length N times q).
the estimated standard deviation of the estimated concentrations (a vector of length N times q).
the desired level (defaults to 0.05, corresponding to a 95% interval).
truncate negative lower limits at zero (defaults to TRUE)
A (1 - \(\alpha\))x100% Wald-type 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 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 = "wald")
#> Error in eval(expr, envir, enclos): object 'mod_summ' not found