| Title: | Utilities for Processing Posterior Samples Stored in 'mcmc.lists' |
|---|---|
| Description: | The aim of 'postpack' is to provide the infrastructure for a standardized workflow for 'mcmc.list' objects. These objects can be used to store output from models fitted with Bayesian inference using 'JAGS', 'WinBUGS', 'OpenBUGS', 'NIMBLE', 'Stan', or even custom MCMC algorithms. Although the 'coda' R package provides some methods for these objects, it is somewhat limited in easily performing post-processing tasks for specific nodes. Models are ever increasing in their complexity and the number of tracked nodes, and oftentimes a user may wish to summarize/diagnose sampling behavior for only a small subset of nodes at a time for a particular question or figure. Thus, many 'postpack' functions support performing tasks on a subset of nodes, where the subset is specified with regular expressions. The functions in 'postpack' streamline the extraction, summarization, and diagnostics of specific monitored nodes after model fitting. Further, because there is rarely only ever one model under consideration, 'postpack' scales efficiently to perform the same tasks on output from multiple models simultaneously, facilitating rapid assessment of model sensitivity to changes in assumptions. |
| Authors: | Ben Staton [aut, cre] |
| Maintainer: | Ben Staton <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.5.4.9000 |
| Built: | 2024-11-04 05:05:11 UTC |
| Source: | https://github.com/bstaton1/postpack |
Use element names to place vector elements in the appropriate location of an array.
array_format(v)array_format(v)
v |
A vector with names indicating the index location of each element in a new array. See the details (particularly the example) for more about what this means. |
Suppose you have an AxB matrix in your model, and you would like to create an object that stores the posterior means in the same AxB matrix as found in the model. For an AxB matrix, this is not too difficult to do "by-hand". However, if there are also dimensions C, D, and E, missing values, etc. it becomes more difficult.
An array with elements of v placed in the appropriate location based on their
index names.
Up to 10 dimensions are currently supported. Please submit an issue should you find that you need more dimensions.
# load example mcmc.list data(cjs) # find an array node from your model match_params(cjs, "SIG") # extract the posterior mean of it SIG_mean = post_summ(cjs, "SIG")["mean",] # note that it has element names SIG_mean # create a matrix with elements in the proper place array_format(SIG_mean)# load example mcmc.list data(cjs) # find an array node from your model match_params(cjs, "SIG") # extract the posterior mean of it SIG_mean = post_summ(cjs, "SIG")["mean",] # note that it has element names SIG_mean # create a matrix with elements in the proper place array_format(SIG_mean)
An example of samples from a joint posterior distribution from a Cormack-Jolly-Seber model. The specific context does not matter, this object is provided to show examples of 'postpack' functionality.
cjscjs
A mcmc.list object.
Posterior samples generated from a model fitted to hypothetical data set.
See vignette("example-mcmclists") on the context, model, and monitored parameters.
An example of samples from a joint posterior distribution from a Cormack-Jolly-Seber model. The specific context does not matter, this object is provided to show examples of 'postpack' functionality.
cjs_no_rhocjs_no_rho
A mcmc.list object.
This object stores samples from the same hypothetical example as for the cjs example object,
with one small change to the model. The rho term that models correlation between slopes and intercepts
was forced to be zero, rather than estimating it. Consult vignette("example-mcmclists") for more details.
Used by diag_plots(), not intended to be called separately
density_plot(post, param, show_diags = "if_poor_Rhat")density_plot(post, param, show_diags = "if_poor_Rhat")
post |
A |
param |
A regular expression that matches a single element in the model.
E.g., |
show_diags |
Control when to display numerical diagnostic summaries on plots. Must be one of
|
A figure showing the posterior density, separated by chain.
This is not a function users will generally use directly. Call diag_plots() instead.
Allows quick visualization of posterior density and trace plots, both separated by chain, for the desired nodes of interest. Includes the ability to plot in the RStudio graphics device, an external device, or a PDF file. Further, with the auto settings, the dimensions of the plotting device scales to the job needed.
diag_plots( post, params, ext_device = FALSE, show_diags = "if_poor_Rhat", layout = "auto", dims = "auto", keep_percent = 1, save = FALSE, file = NULL, auto_escape = TRUE )diag_plots( post, params, ext_device = FALSE, show_diags = "if_poor_Rhat", layout = "auto", dims = "auto", keep_percent = 1, save = FALSE, file = NULL, auto_escape = TRUE )
post |
A |
params |
A vector of regular expressions specifying the nodes to match for plotting.
Accepts multi-element vectors to match more than one node at a time.
See |
ext_device |
Display plots in an external device rather than the active device?
|
show_diags |
Control when to display numerical diagnostic summaries on plots. Must be one of
|
layout |
Control how parameter diagnostics are organized into |
dims |
Control the dimensions of the graphics device using |
keep_percent |
Proportion (between 0 and 1) of samples to keep for trace plotting.
Passed to |
save |
Save the diagnostic plots in a PDF file? If so,
specify |
file |
File name of a PDF file to save the plots to.
Must include the |
auto_escape |
Automatically escape |
A multi-panel figure showing the posterior density and trace plots for requested nodes. The device in which it is placed depends on the argument values.
If saving as a pdf, these files can get very large with many samples and render slowly.
The keep_percent argument is intended to help with this by thinning the chains at quasi-evenly spaced intervals.
match_params(), density_plot(), trace_plot()
if (interactive()) { #load example mcmc.list data(cjs) # use current device diag_plots(cjs, "B0") # use a new device diag_plots(cjs, "B0", ext_device = TRUE) # always show diagnostic summaries diag_plots(cjs, "B0", show_diags = "always") # use a different layout (leaving it as "auto" is usually best) diag_plots(cjs, c("sig", "b"), layout = "5x3") # save diagnostics for all nodes to a pdf file diag_plots(cjs, "", save = TRUE, file = "diags.pdf") }if (interactive()) { #load example mcmc.list data(cjs) # use current device diag_plots(cjs, "B0") # use a new device diag_plots(cjs, "B0", ext_device = TRUE) # always show diagnostic summaries diag_plots(cjs, "B0", show_diags = "always") # use a different layout (leaving it as "auto" is usually best) diag_plots(cjs, c("sig", "b"), layout = "5x3") # save diagnostics for all nodes to a pdf file diag_plots(cjs, "", save = TRUE, file = "diags.pdf") }
Removes square brackets, numbers, and commas that represent the index of the node element in question. Returns just the node name.
drop_index(params)drop_index(params)
params |
Node names. |
A character vector with the same length as params, with no indices included.
For example, "a[1]" becomes "a".
This is not a function users will generally use directly.
Returns the names of all quantities stored in
a mcmc.list object.
get_params(post, type = "base_only")get_params(post, type = "base_only")
post |
A |
type |
Format of returned matches; only two options are accepted:
|
A character vector with all node names stored in the post object, formatted as requested by type.
# load example mcmc.list data(cjs) # get only node names, no indices (default) get_params(cjs, type = "base_only") # get indices too, where applicable get_params(cjs, type = "base_index")# load example mcmc.list data(cjs) # get only node names, no indices (default) get_params(cjs, type = "base_only") # get indices too, where applicable get_params(cjs, type = "base_index")
Extract chain and iteration IDs for each sample
id_mat(post)id_mat(post)
post |
A |
A matrix with columns "CHAIN" and "ITER".
This is not a function users will generally use directly.
Insert escapes on regex brackets
ins_regex_bracket(params)ins_regex_bracket(params)
params |
Node names. |
Searches the contents of a string for the occurrence of a square bracket or two, and inserts the necessary escapes for pattern matching via regular expressions.
A character vector with all brackets escaped. For example,
"a[1]" becomes "a\\[1\\]"
This is not a function users will generally use directly.
To ensure that a regular expression will match exactly, it's necessary to specify so.
ins_regex_lock(params)ins_regex_lock(params)
params |
Node names to paste a |
A character vector with locking anchors inserted to force an exact match. For example,
"a\\[1\\]" becomes "^a\\[1\\]$".
This is not a function users will generally use directly.
Converts a vector into a comma-separated list for use in sentences (error messages, warnings, etc.).
list_out(x, final = NULL, per_line = 1e+06, wrap = NULL, indent = NULL)list_out(x, final = NULL, per_line = 1e+06, wrap = NULL, indent = NULL)
x |
A vector, will be coerced to a character. |
final |
Word that will separate the final element in the list from others. See the examples. |
per_line |
Number of elements printed per line. See the examples. |
wrap |
Optional character to wrap around each element, e.g., quotation marks. |
indent |
Optional string to place in front of the first element on each line. See the examples. |
A character vector with length == 1; ready to be passed to
base::stop(), base::warning(), or base::cat(), to provide a useful message.
Returns all the node names stored in a mcmc.list
object that match a provided string.
match_params(post, params, type = "base_index", auto_escape = TRUE)match_params(post, params, type = "base_index", auto_escape = TRUE)
post |
A |
params |
A vector of regular expressions specifying the nodes to match.
Accepts multi-element vectors to match more than one node at a time.
See examples and |
type |
Format of returned matches; only two options are accepted:
|
auto_escape |
Automatically escape |
This function is called as one of the first steps in many of the more downstream functions in this package. It is thus fairly important to get used to how the regular expressions work. This function can be used directly to hone in on the correct regular expression. See the examples below.
A character vector with all node names in post that match params, formatted as requested by type..
If no matches are found, an error will be returned with
the base node names found in post to help the next try.
# load example mcmc.list data(cjs) # these produce same output b/c of regex pattern matching match_params(cjs, params = c("b0", "b1")) match_params(cjs, params = c("b")) # return only base names, not indices as well match_params(cjs, params = "b", type = "base_only") # force a match to start with B match_params(cjs, "^B") # force a match to end with 0 match_params(cjs, "0$") # use a wild card to get b0[3] and b1[3] match_params(cjs, "b.[3]") # repeat a wild card match_params(cjs, "s.+0") # turn off auto escape to use [] in regex syntax rather than matching them as text match_params(cjs, params = "[:digit:]$", auto_escape = FALSE) # pass a dot to match all (same as get_params) match_params(cjs, ".")# load example mcmc.list data(cjs) # these produce same output b/c of regex pattern matching match_params(cjs, params = c("b0", "b1")) match_params(cjs, params = c("b")) # return only base names, not indices as well match_params(cjs, params = "b", type = "base_only") # force a match to start with B match_params(cjs, "^B") # force a match to end with 0 match_params(cjs, "0$") # use a wild card to get b0[3] and b1[3] match_params(cjs, "b.[3]") # repeat a wild card match_params(cjs, "s.+0") # turn off auto escape to use [] in regex syntax rather than matching them as text match_params(cjs, params = "[:digit:]$", auto_escape = FALSE) # pass a dot to match all (same as get_params) match_params(cjs, ".")
Used by diag_plots() to place a common
title over top of two figures: one density and one trace
for a given node.
mytitle(text)mytitle(text)
text |
The text string to include as a centered title over two adjacent plots. |
This is not a function users will generally use directly.
Intended for use when derived quantities are calculated from monitored posterior samples,
and you wish to combine them into the master mcmc.list,
as though they were calculated and monitored during MCMC sampling.
It is not advised to combine samples from two MCMC runs, because covariance
of MCMC sampling would be lost.
post_bind(post1, post2, dup_id = "_p2")post_bind(post1, post2, dup_id = "_p2")
post1 |
|
post2 |
|
dup_id |
If any node names are duplicated in |
Some important things to note:
If the object passed to post1 is a matrix, post2 must be a mcmc.list, and vice versa.
That is, two mcmc.list objects are allowed, but not two matrix objects.
For matrix objects, nodes should be stored as columns and samples should be stored as rows. Column names should be present.
The objects passed to post1 and post2 must have the same number of chains, iterations, burnin, and thinning interval.
If the node names are empty (e.g., missing column names in a matrix), the node names will be coerced to "var1", "var2", etc. and a warning will be returned.
A single mcmc.list object containing samples of the nodes from both post1 and post2.
# load example mcmc.list data(cjs) # create two subsets from cjs: one as mcmc.list and one as matrix # also works if both are mcmc.list objects p1 = post_subset(cjs, "b0") p2 = post_subset(cjs, "b1", matrix = TRUE) # combine them into one mcmc.list head(post_bind(p1, p2))# load example mcmc.list data(cjs) # create two subsets from cjs: one as mcmc.list and one as matrix # also works if both are mcmc.list objects p1 = post_subset(cjs, "b0") p2 = post_subset(cjs, "b1", matrix = TRUE) # combine them into one mcmc.list head(post_bind(p1, p2))
Wrapper around several ways of converting objects to mcmc.list format,
automated based on the input object class.
post_convert(obj)post_convert(obj)
obj |
An object storing posterior samples from an MCMC algorithm.
Accepted classes are |
Accepted classes are produced by several packages, including but probably not limited to:
stanfit objects are created by rstan::stan(), which also provides rstan::As.mcmc.list(). Rather than requiring users to have 'rstan' installed to use 'postpack', post_convert() will instruct users to use this function if supplied a stanfit object.
bugs objects are created by R2WinBUGS::bugs() and R2OpenBUGS::bugs().
rjags objects are created by R2jags::jags().
list objects are created by nimble::runMCMC(), 'MCMCpack' functions, or custom MCMC algorithms.
matrix objects are created by post_subset() and is often the format of posterior quantities derived from monitored nodes.
mcmc.list objects are created by rjags::coda.samples(), jagsUI::jags.basic(), and jagsUI::jags()$samples. If a mcmc.list object is passed to obj, an error will be returned telling the user this function is not necessary.
If you find that a critical class conversion is missing, please submit an issue requesting its addition with a minimum working example of how it can be created.
The same information as passed in the obj argument, but formatted as mcmc.list class.
If samples are stored in a list object, the individual elements must be matrix or mcmc class, storing the samples (rows) across parameters (columns, with names) for each chain (list elements). If list elements are in matrix format, they will be coerced to mcmc format, and thinning, start, and end intervals may be inaccurate.
If samples are stored in a matrix object, rows should store samples and columns should store nodes. Multiple chains should be combined using base::rbind(). Two additional columns must be present: "CHAIN" and "ITER", which denote the MCMC chain and iteration numbers, respectively.
coda::as.mcmc.list(), coda::as.mcmc()
## EXAMPLE 1 # load example mcmc.list data(cjs) # take a subset from cjs as a matrix, retain chain and iter ids cjs_sub = post_subset(cjs, "^B", matrix = TRUE, chains = TRUE, iters = TRUE) # convert back to mcmc.list class(post_convert(cjs_sub)) ## EXAMPLE 2: create mcmc.list from hypothetical MCMC samples; chains are list elements # create hypothetical samples; can't use postpack on this - not an mcmc.list samps = lapply(1:3, function(i) { m = matrix(rnorm(100), 20, 5) colnames(m) = paste0("param", 1:5) m }) # convert samps_new = post_convert(samps) # can use postpack now post_summ(samps_new, "param") ## EXAMPLE 3: create mcmc.list from hypothetical MCMC samples; chains rbind-ed matrices # create samples f = function() { m = matrix(rnorm(100), 20, 5) colnames(m) = paste0("param", 1:5) m } samps = rbind(f(), f(), f()) # assign chain and iter IDs to each sample samps = cbind(CHAIN = rep(1:3, each = 20), ITER = rep(1:20, 3), samps) # convert samps_new = post_convert(samps) # can use postpack now post_summ(samps_new, "param")## EXAMPLE 1 # load example mcmc.list data(cjs) # take a subset from cjs as a matrix, retain chain and iter ids cjs_sub = post_subset(cjs, "^B", matrix = TRUE, chains = TRUE, iters = TRUE) # convert back to mcmc.list class(post_convert(cjs_sub)) ## EXAMPLE 2: create mcmc.list from hypothetical MCMC samples; chains are list elements # create hypothetical samples; can't use postpack on this - not an mcmc.list samps = lapply(1:3, function(i) { m = matrix(rnorm(100), 20, 5) colnames(m) = paste0("param", 1:5) m }) # convert samps_new = post_convert(samps) # can use postpack now post_summ(samps_new, "param") ## EXAMPLE 3: create mcmc.list from hypothetical MCMC samples; chains rbind-ed matrices # create samples f = function() { m = matrix(rnorm(100), 20, 5) colnames(m) = paste0("param", 1:5) m } samps = rbind(f(), f(), f()) # assign chain and iter IDs to each sample samps = cbind(CHAIN = rep(1:3, each = 20), ITER = rep(1:20, 3), samps) # convert samps_new = post_convert(samps) # can use postpack now post_summ(samps_new, "param")
Quickly query the number of burn-in samples, post-burnin, thinning,
number of chains, etc. from a mcmc.list object.
post_dim(post, types = NULL)post_dim(post, types = NULL)
post |
A |
types |
The dimension types to return. Must contain some of |
A numeric vector with named elements, which may contain:
burn: The burn-in period + adapting phase (per chain).
post_burn: The post-burn-in period (per chain).
thin: The thinning interval post-burn-in.
chains: The number of chains.
saved: The number of saved samples across all chains.
params: The number of nodes with MCMC samples.
All of these will be returned if types = NULL, a subset can be returned by
specifying (for example) types = c("burn", "thin").
If the post object was thinned after MCMC completed
using post_thin(), then the "burn" and "thin" dimensions will be improperly calculated.
post_thin() performs post-MCMC thinning of mcmc.list objects,
and is solely for developing long-running post-processing code, so this is okay.
# load example mcmc.list data(cjs) # get all relevant dimensions post_dim(cjs) # get only the number of chains post_dim(cjs, "chains") # get the thinning and burn-in intervals post_dim(cjs, c("burn", "thin"))# load example mcmc.list data(cjs) # get all relevant dimensions post_dim(cjs) # get only the number of chains post_dim(cjs, "chains") # get the thinning and burn-in intervals post_dim(cjs, c("burn", "thin"))
Just like post_subset(), but keep all nodes except those that match.
post_remove(post, params, ask = TRUE, auto_escape = TRUE)post_remove(post, params, ask = TRUE, auto_escape = TRUE)
post |
A |
params |
A vector of regular expressions specifying the nodes to match for removal.
Accepts multi-element vectors to match more than one node at a time.
See |
ask |
Prompt user for a response prior to removing nodes? |
auto_escape |
Automatically escape |
A mcmc.list, identical in all ways to the original
except that nodes matched by the params argument are removed. If the user
responds "no" to the question when ask = TRUE, post is returned unaltered.
# load example mcmc.list data(cjs) # get names of all nodes get_params(cjs) # remove the SIG nodes new_cjs = suppressMessages(post_remove(cjs, "SIG", ask = FALSE)) # get names of new output get_params(new_cjs)# load example mcmc.list data(cjs) # get names of all nodes get_params(cjs) # remove the SIG nodes new_cjs = suppressMessages(post_remove(cjs, "SIG", ask = FALSE)) # get names of new output get_params(new_cjs)
Subsets a smaller portion from a mcmc.list object
corresponding only to the node(s) requested.
post_subset( post, params, matrix = FALSE, iters = FALSE, chains = FALSE, auto_escape = TRUE )post_subset( post, params, matrix = FALSE, iters = FALSE, chains = FALSE, auto_escape = TRUE )
post |
A |
params |
A vector of regular expressions specifying the nodes to match for subsetting.
Accepts multi-element vectors to match more than one node at a time.
See |
matrix |
|
iters |
Retain the iteration number of each sample if |
chains |
Retain the chain number of each sample if |
auto_escape |
Automatically escape |
A mcmc.list or matrix object, depending on the
value of the matrix argument. Object contains all nodes that match the params argument;
an error will be returned if no matches are found.
# load example mcmc.list data(cjs) # create mcmc.list with all nodes that contain "B0" x1 = post_subset(cjs, "B0") # create mcmc.list with all nodes that contain "b" or "B" x2 = post_subset(cjs, c("b", "B")) # perform the subset and return a matrix as output, while retaining the chain ID x3 = post_subset(cjs, "B0", matrix = TRUE, chain = TRUE)# load example mcmc.list data(cjs) # create mcmc.list with all nodes that contain "B0" x1 = post_subset(cjs, "B0") # create mcmc.list with all nodes that contain "b" or "B" x2 = post_subset(cjs, c("b", "B")) # perform the subset and return a matrix as output, while retaining the chain ID x3 = post_subset(cjs, "B0", matrix = TRUE, chain = TRUE)
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 )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 )
post |
A |
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 |
digits |
Control rounding of summaries.
Passed to |
probs |
Posterior quantiles to calculate. Passed to |
Rhat |
Calculate the Rhat convergence diagnostic using |
neff |
Calculate the number of effective MCMC samples using |
mcse |
Calculate the Monte Carlo standard error for the posterior mean and reported quantiles
using the |
by_chain |
Calculate posterior summaries for each chain
rather than for the aggregate across chains? Defaults to |
auto_escape |
Automatically escape |
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.
match_params(), coda::gelman.diag(), coda::effectiveSize(), mcmcse::mcse(), mcmcse::mcse.q()
# load example mcmc.list data(cjs) # calculate posterior summaries for the "p" nodes # ("p[1]" doesn't exist in model) post_summ(cjs, "p") # do this by chain post_summ(cjs, "p", by_chain = TRUE) # calculate Rhat and Neff diagnostic summaries as well # multiple node names too post_summ(cjs, c("b0", "p"), Rhat = TRUE, neff = TRUE) # calculate Monte Carlo SE for mean and quantiles, with rounding post_summ(cjs, "p", mcse = TRUE, digits = 3) # summarize different quantiles: median and central 80% post_summ(cjs, "p", probs = c(0.5, 0.1, 0.9))# load example mcmc.list data(cjs) # calculate posterior summaries for the "p" nodes # ("p[1]" doesn't exist in model) post_summ(cjs, "p") # do this by chain post_summ(cjs, "p", by_chain = TRUE) # calculate Rhat and Neff diagnostic summaries as well # multiple node names too post_summ(cjs, c("b0", "p"), Rhat = TRUE, neff = TRUE) # calculate Monte Carlo SE for mean and quantiles, with rounding post_summ(cjs, "p", mcse = TRUE, digits = 3) # summarize different quantiles: median and central 80% post_summ(cjs, "p", probs = c(0.5, 0.1, 0.9))
Removes iterations from each chain of a mcmc.list
object at quasi-evenly spaced intervals. Post-MCMC thinning is useful for
developing long-running post-processing code with a smaller but otherwise identical mcmc.list.
post_thin(post, keep_percent, keep_iters)post_thin(post, keep_percent, keep_iters)
post |
A |
keep_percent |
Proportion (between 0 and 1) of samples to keep from each chain.
Setting |
keep_iters |
Number of samples to keep from each chain. |
The samples will be removed at as evenly spaced intervals as possible, however, this is not perfect. It is therefore recommended to use the full posterior for final post-processing calculations, but this should be fine for most development of long-running code.
If both keep_percent and keep_iters are supplied, an error will be returned requesting that only
one be used.
A mcmc.list object, identical to post, but with fewer samples of each node.
Iteration numbers are reset after thinning the samples. So if running post_dim()
on output passed through post_thin(), you cannot trust the burn-in or thinning counts.
Again, this is not an issue for developing post-processing code.
# load example mcmc.list data(cjs) # take note of original dimensions post_dim(cjs) # keep ~20% of the samples cjs_thin1 = post_thin(cjs, keep_percent = 0.2) # note burn-in and thin intervals no longer correct! # but desired outcome achieved - identical object but smaller post_dim(cjs_thin1) # keep 30 samples per chain cjs_thin2 = post_thin(cjs, keep_iters = 30) post_dim(cjs_thin2)# load example mcmc.list data(cjs) # take note of original dimensions post_dim(cjs) # keep ~20% of the samples cjs_thin1 = post_thin(cjs, keep_percent = 0.2) # note burn-in and thin intervals no longer correct! # but desired outcome achieved - identical object but smaller post_dim(cjs_thin1) # keep 30 samples per chain cjs_thin2 = post_thin(cjs, keep_iters = 30) post_dim(cjs_thin2)
The aim of 'postpack' is to provide the infrastructure for a standardized workflow for mcmc.list objects.
These objects can be used to store output from models fitted with Bayesian inference using
JAGS, Win/OpenBUGS, NIMBLE, Stan, or even custom MCMC algorithms (see post_convert() for converting samples to
mcmc.list format). Although the 'coda' package provides some methods for these objects,
it is somewhat limited in easily performing post-processing tasks for particular nodes.
Models are ever increasing in their complexity and the number of tracked nodes, and oftentimes
a user may wish to summarize/diagnose sampling behavior for only a small subset of nodes at a time for a particular question or figure.
Thus, many 'postpack' functions support performing tasks on a subset of nodes, where the subset is specified with regular expressions.
The functions in this package streamline the extraction, summarization, and diagnostics of particular nodes monitored after model fitting.
Further, because there is rarely only ever one model under consideration, 'postpack' scales efficiently to perform the same tasks on output from multiple models
simultaneously, facilitating rapid assessment of model sensitivity to changes in assumptions.
Remove escapes on regex brackets
rm_regex_bracket(params)rm_regex_bracket(params)
params |
Node names. |
Searches the contents of a string for the occurrence of a square bracket or two (that has been escaped), and removes the escaping that was necessary for matching via regular expressions.
A character vector with all brackets escaped. For example,
"a\\[1\\]" becomes "a[1]".
This is not a function users will generally use directly.
Undoes the work of ins_regex_lock().
rm_regex_lock(params)rm_regex_lock(params)
params |
Node names to remove a |
A character vector with locking anchors inserted to force an exact match. For example,
"^a\\[1\\]$" becomes "a\\[1\\]".
This is not a function users will generally use directly.
Create a trace plot for a single desired node
trace_plot(post, param, keep_percent = 1)trace_plot(post, param, keep_percent = 1)
post |
A |
param |
A regular expression that matches a single element in the model.
E.g., |
keep_percent |
A numeric vector of length == 1 and on the range (0,1].
Percent of samples you'd like to keep for trace plotting and passed to |
If saving as a pdf file, these files can get very large with many samples and render slowly.
The keep_percent argument is intended to help with this by thinning the chains at quasi-evenly spaced intervals.
This is not a function users will generally use directly. Call diag_plots() instead.
For each posterior sample, extract the standard deviation and correlation components of a monitored node representing a variance-covariance matrix.
vcov_decomp( post, param, sigma_base_name = "sigma", rho_base_name = "rho", invert = FALSE, check = TRUE, auto_escape = TRUE )vcov_decomp( post, param, sigma_base_name = "sigma", rho_base_name = "rho", invert = FALSE, check = TRUE, auto_escape = TRUE )
post |
A |
param |
A vector of regular expressions specifying the nodes to match for plotting.
Must match only one base node name in |
sigma_base_name |
Base node name to assign to the standard deviation vector component?
Defaults to |
rho_base_name |
Same as |
invert |
Take the inverse of the matrix node matched by |
check |
Perform checks sequentially that the matrix node is (a) square, (b) symmetrical, and (c) positive definite
before proceeding with the calculations? If set to |
auto_escape |
Automatically escape |
A mcmc.list object.
# load example mcmc.list data(cjs) # "SIG" is a covariance matrix node SIG_decomp = vcov_decomp(cjs, "SIG") # extract the posterior mean correlation matrix, and reformat array_format(post_summ(SIG_decomp, "rho")["mean",])# load example mcmc.list data(cjs) # "SIG" is a covariance matrix node SIG_decomp = vcov_decomp(cjs, "SIG") # extract the posterior mean correlation matrix, and reformat array_format(post_summ(SIG_decomp, "rho")["mean",])
Performs the same basic function as R2OpenBUGS::write.model()
write_model(fun, file)write_model(fun, file)
fun |
A function object containing BUGS/JAGS model code |
file |
A character vector of length == 1: the name of the file to write to |
Performs the same basic function as R2OpenBUGS::write.model(),
but with slightly better output (scientific notation, spacing, etc.). The main reason it was created
for use in 'postpack' was to remove the need for using the 'R2OpenBUGS' package when not using OpenBUGS.
Nothing, but file is written to disk.
if (interactive()) { # define some simple BUGS model as an R function # note the use of %_% to include a truncation mod = function() { # PRIORS mu ~ dnorm(0,0.001) %_% T(0,) sig ~ dunif(0,10) tau <- 1/sig^2 # LIKELIHOOD for (i in 1:n) { y[i] ~ dnorm(mu, tau) } } # write model to a text file to be called by BUGS/JAGS write_model(mod, "model.txt") }if (interactive()) { # define some simple BUGS model as an R function # note the use of %_% to include a truncation mod = function() { # PRIORS mu ~ dnorm(0,0.001) %_% T(0,) sig ~ dunif(0,10) tau <- 1/sig^2 # LIKELIHOOD for (i in 1:n) { y[i] ~ dnorm(mu, tau) } } # write model to a text file to be called by BUGS/JAGS write_model(mod, "model.txt") }