Obtain posterior summaries and diagnostics of specific nodes

Allows rapid calculation of summaries and diagnostics from specific nodes stored in mcmc.list objects.

post_summ(
  post,
  params,
  digits = NULL,
  probs = c(0.5, 0.025, 0.975),
  Rhat = FALSE,
  neff = FALSE,
  mcse = FALSE,
  by_chain = FALSE,
  auto_escape = TRUE
)

Arguments

post	A mcmc.list object.
params	A vector of regular expressions specifying the nodes to match for summarization. Accepts multi-element vectors to match more than one node at a time. See match_params() and vignette("pattern-matching") for more details.
digits	Control rounding of summaries. Passed to base::round() and defaults to NULL, which produces no rounding.
probs	Posterior quantiles to calculate. Passed to stats::quantile(). Defaults to probs = c(0.5, 0.025, 0.975) (i.e., median and equal-tailed 95 percent credible interval).
Rhat	Calculate the Rhat convergence diagnostic using coda::gelman.diag()? Fair warning: this can take a bit of time to run on many nodes/samples.
neff	Calculate the number of effective MCMC samples using coda::effectiveSize()? Fair warning: this can take a bit of time to run on many nodes/samples.
mcse	Calculate the Monte Carlo standard error for the posterior mean and reported quantiles using the mcmcse::mcse() and mcmcse::mcse.q() functions (batch means method with batch size automatically calculated)? Fair warning: this can take a bit of time to run on many nodes/samples.
by_chain	Calculate posterior summaries for each chain rather than for the aggregate across chains? Defaults to FALSE. The arguments Rhat, neff, and mcse are ignored if by_chain = TRUE and a warning will be returned.
auto_escape	Automatically escape "[" and "]" characters for pattern matching? See match_params() for details.

Value

A matrix object with summary statistics as rows and nodes as columns. If by_chain = TRUE, an array with chain-specific summaries as the third dimension is returned instead.

See also

match_params(), coda::gelman.diag(), coda::effectiveSize(), mcmcse::mcse(), mcmcse::mcse.q()

Examples

# load example mcmc.list
data(cjs)

# calculate posterior summaries for the "p" nodes
# ("p[1]" doesn't exist in model)
post_summ(cjs, "p")
#>             p[2]       p[3]       p[4]
#> mean  0.74022487 0.70052497 0.63349268
#> sd    0.02365831 0.03105759 0.04401473
#> 50%   0.74081317 0.70128371 0.63190886
#> 2.5%  0.69103388 0.64196825 0.54655802
#> 97.5% 0.78878642 0.76240151 0.71657359

# do this by chain
post_summ(cjs, "p", by_chain = TRUE)
#> , , chain1
#> 
#>             p[2]       p[3]       p[4]
#> mean  0.74111083 0.70041115 0.63057620
#> sd    0.02387815 0.03059284 0.04453106
#> 50%   0.74136511 0.69943805 0.62961609
#> 2.5%  0.68882162 0.64444668 0.54139258
#> 97.5% 0.78842368 0.76322104 0.71283345
#> 
#> , , chain2
#> 
#>             p[2]       p[3]       p[4]
#> mean  0.73933892 0.70063879 0.63640916
#> sd    0.02345074 0.03157647 0.04338523
#> 50%   0.74004196 0.70301389 0.63447834
#> 2.5%  0.69499776 0.63981239 0.55475454
#> 97.5% 0.78828940 0.76064092 0.71673527
#> 

# calculate Rhat and Neff diagnostic summaries as well
# multiple node names too
post_summ(cjs, c("b0", "p"), Rhat = TRUE, neff = TRUE)
#>             b0[1]       b0[2]      b0[3]       b0[4]       b0[5]         p[2]
#> mean    1.3753033   2.3192679   1.751486   1.5200421   1.0824184   0.74022487
#> sd      0.2067011   0.3564031   0.238889   0.2032504   0.2077227   0.02365831
#> 50%     1.3674581   2.2779459   1.728931   1.5075479   1.0754654   0.74081317
#> 2.5%    0.9614656   1.7836702   1.324638   1.1597562   0.6748646   0.69103388
#> 97.5%   1.8053107   3.1633955   2.287703   1.9605019   1.5194312   0.78878642
#> Rhat    1.0000000   0.9980000   1.004000   0.9990000   0.9990000   1.00000000
#> neff  433.0000000 456.0000000 474.000000 463.0000000 570.0000000 449.00000000
#>               p[3]         p[4]
#> mean    0.70052497   0.63349268
#> sd      0.03105759   0.04401473
#> 50%     0.70128371   0.63190886
#> 2.5%    0.64196825   0.54655802
#> 97.5%   0.76240151   0.71657359
#> Rhat    0.99900000   1.00500000
#> neff  500.00000000 500.00000000

# calculate Monte Carlo SE for mean and quantiles, with rounding
post_summ(cjs, "p", mcse = TRUE, digits = 3)
#>             p[2]  p[3]  p[4]
#> mean       0.740 0.701 0.633
#> sd         0.024 0.031 0.044
#> 50%        0.741 0.701 0.632
#> 2.5%       0.691 0.642 0.547
#> 97.5%      0.789 0.762 0.717
#> mcse_mean  0.001 0.001 0.002
#> mcse_50%   0.001 0.002 0.003
#> mcse_2.5%  0.004 0.003 0.006
#> mcse_97.5% 0.002 0.003 0.003

# summarize different quantiles: median and central 80%
post_summ(cjs, "p", probs = c(0.5, 0.1, 0.9))
#>            p[2]       p[3]       p[4]
#> mean 0.74022487 0.70052497 0.63349268
#> sd   0.02365831 0.03105759 0.04401473
#> 50%  0.74081317 0.70128371 0.63190886
#> 10%  0.70861186 0.66055910 0.57567563
#> 90%  0.76915547 0.74189193 0.68970096