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

tidy.MCMCglmmR Documentation

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 Oct. 16, 2024, 1:06 a.m.