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" )
| stan_mod | the model summary object from Stan. |
|---|---|
| stan_samps | the list of MCMC samples from Stan. |
| taxa_of_interest | the indices of the taxa for which point estimates and posterior summaries are desired. |
| mult_num | the number to multiply the resulting estimates and standard deviations by (defaults to 1). |
| level | the |
| interval_type | 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 9.3e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.93 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.008434 seconds (Warm-up) #> Chain 1: 0.001527 seconds (Sampling) #> Chain 1: 0.009961 seconds (Total) #> Chain 1:#> Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See #> http://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 #> http://mc-stan.org/misc/warnings.html#r-hat#> Error in rstan::summary(mod, probs = c(0.025, 0.975))$summary: $ operator is invalid for atomic vectorsmod_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 rownames(stan_mod): object 'mod_summ' not found