mcmc_tidiers: Tidying methods for MCMC (Stan, JAGS, etc.) fits

In broom.mixed: Tidying Methods for Mixed Models

Tidying methods for MCMC (Stan, JAGS, etc.) fits

Description

Tidying methods for MCMC (Stan, JAGS, etc.) fits

Usage

## S3 method for class 'MCMCglmm'
tidy(x, effects = c("fixed", "ran_pars"), scales = NULL, ...)

tidyMCMC(
  x,
  pars,
  robust = FALSE,
  conf.int = FALSE,
  conf.level = 0.95,
  conf.method = c("quantile", "HPDinterval"),
  drop.pars = c("lp__", "deviance"),
  rhat = FALSE,
  ess = FALSE,
  index = FALSE,
  ...
)

## S3 method for class 'rjags'
tidy(
  x,
  robust = FALSE,
  conf.int = FALSE,
  conf.level = 0.95,
  conf.method = "quantile",
  ...
)

## S3 method for class 'stanfit'
tidy(
  x,
  pars,
  robust = FALSE,
  conf.int = FALSE,
  conf.level = 0.95,
  conf.method = c("quantile", "HPDinterval"),
  drop.pars = c("lp__", "deviance"),
  rhat = FALSE,
  ess = FALSE,
  index = FALSE,
  ...
)

## S3 method for class 'mcmc'
tidy(
  x,
  pars,
  robust = FALSE,
  conf.int = FALSE,
  conf.level = 0.95,
  conf.method = c("quantile", "HPDinterval"),
  drop.pars = c("lp__", "deviance"),
  rhat = FALSE,
  ess = FALSE,
  index = FALSE,
  ...
)

## S3 method for class 'mcmc.list'
tidy(
  x,
  pars,
  robust = FALSE,
  conf.int = FALSE,
  conf.level = 0.95,
  conf.method = c("quantile", "HPDinterval"),
  drop.pars = c("lp__", "deviance"),
  rhat = FALSE,
  ess = FALSE,
  index = FALSE,
  ...
)

Arguments

x	a model fit to be converted to a data frame
effects	which effects (fixed, random, etc.) to return
scales	scales on which to report results
...	mostly unused; for tidy.MCMCglmm, these represent options passed through to tidy.mcmc (e.g. robust, conf.int, conf.method, ...)
pars	(character) specification of which parameters to include
robust	use mean and standard deviation (if FALSE) or median and mean absolute deviation (if TRUE) to compute point estimates and uncertainty?
conf.int	(logical) include confidence interval?
conf.level	probability level for CI
conf.method	method for computing confidence intervals ("quantile" or "HPDinterval")
drop.pars	Parameters not to include in the output (such as log-probability information)
rhat, ess	(logical) include Rhat and/or effective sample size estimates?
index	Add index column, remove index from term. For example, term a[13] becomes term a and index 13.

Examples

if (require("MCMCglmm")) {
  ## original model
  ## Not run: 
      mm0 <- MCMCglmm(Reaction ~ Days,
                 random = ~Subject, data = sleepstudy,
                 nitt=4000,
                 pr = TRUE
             )
   
## End(Not run)
   ## load stored object
   load(system.file("extdata","MCMCglmm_example.rda",
                                     package="broom.mixed"))
   tidy(mm0)
   tidy(mm1)
   tidy(mm2)
   tail(tidy(mm0,effects="ran_vals"))
}

# Using example from "RStan Getting Started"
# https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started

model_file <- system.file("extdata", "8schools.stan", package = "broom.mixed")
schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))
## original model
## Not run: 
    set.seed(2015)
    rstan_example <- rstan::stan(file = model_file, data = schools_dat,
                         iter = 1000, chains = 2, save_dso = FALSE)

## End(Not run)
if (require(rstan)) {
   ## load stored object
   rstan_example <- readRDS(system.file("extdata", "rstan_example.rds", package = "broom.mixed"))
   tidy(rstan_example)
   tidy(rstan_example, conf.int = TRUE, pars = "theta")
   td_mean <- tidy(rstan_example, conf.int = TRUE)
   td_median <- tidy(rstan_example, conf.int = TRUE, robust = TRUE)
   
   if (require(dplyr) && require(ggplot2)) {
       tds <- (dplyr::bind_rows(list(mean=td_mean, median=td_median), .id="method")
          %>% mutate(type=ifelse(grepl("^theta",term),"theta",
            ifelse(grepl("^eta",term),"eta",
                  "other")))
      )

     ggplot(tds, aes(estimate, term)) +
      geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),height=0) +
      geom_point(aes(color = method))+
      facet_wrap(~type,scale="free",ncol=1)
 } ## require(dplyr,ggplot2)
} ## require(rstan)
if (require(R2jags)) {
   ## see help("jags",package="R2jags")
   ## and  example("jags",package="R2jags")
   ## for details
   ## load stored object
   R2jags_example <- readRDS(system.file("extdata", "R2jags_example.rds", package = "broom.mixed"))
   tidy(R2jags_example)
   tidy(R2jags_example, conf.int=TRUE, conf.method="quantile")
}

broom.mixed documentation built on April 18, 2022, 1:06 a.m.