R/extract_posterior_summaries.R
extract_posterior_summaries.RdReturn point estimates and credible intervals for the true concentration, and point estimates and prediction intervals for estimated qPCR obtained through a Stan fit.
extract_posterior_summaries(
stan_mod,
stan_samps,
taxa_of_interest,
mult_num = 1,
level = 0.95,
interval_type = "wald"
)the model summary object from Stan.
the list of MCMC samples from Stan.
the indices of the taxa for which point estimates and posterior summaries are desired.
the number to multiply the resulting estimates and standard deviations by (defaults to 1).
the alpha level for prediction intervals (defaults to 0.95, for a nominal 95% prediction interval).
the type of prediction interval desired (defaults to "wald", but "quantile" is also acceptable).
An object of class paramedic. See Details for more information
A paramedic object is a list containing the following elements:
estimates - the point estimates of qPCR (a matrix with dimension sample size by number of taxa).
pred_intervals - predction intervals for qPCR (an array with dimension sample size by 2 by number of taxa).
est_efficiency - point estimates for estimated varying efficiency, if varying efficiency was modeled (a vector of length number of taxa); otherwise, NA.
efficiency_intervals - posterior level level\(\times\)100% confidence intervals for the true efficiency, if efficiency was modeled (a matrix of dimension number of taxa by 2); otherwise, NA.
# 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.00021 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.1 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 summary, samples
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
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