Various edits to MCMC tidiers; mostly style changes. Added 8schools.s…

…tan and rstan_example.rda to extdata, to serve as examples for MCMC tidiers.

DESCRIPTION

@@ -29,8 +29,7 @@ Imports:
     psych,
     stringr,
     reshape2,
-    nlme,
-    coda
+    nlme
 Suggests:
     knitr,
     boot,
@@ -63,7 +62,8 @@ Suggests:
     Hmisc,
     bbmle,
     gamlss,
-    rstan
+    rstan,
+    coda
 URL: http://github.com/dgrtwo/broom
 BugReports: http://github.com/dgrtwo/broom/issues
 VignetteBuilder: knitr

NAMESPACE

@@ -147,8 +147,6 @@ export(inflate)
 export(tidy)
 export(tidyMCMC)
 import(dplyr)
-importFrom(coda,HPDinterval)
-importFrom(coda,as.mcmc)
 importFrom(nlme,VarCorr)
 importFrom(nlme,ranef)
 importFrom(plyr,ldply)

R/mcmc_tidiers.R

@@ -1,66 +1,118 @@
 ## FIXME: design question -- how to make tidying methods inherit appropriately?
+
 #' Tidying methods for MCMC (Stan, JAGS, etc.) fits
-##'
-##' @param x an object of class \sQuote{"stanfit"}
-##' @param pars (character) specification of which parameters to include
-##' @param pt.method method for computing point estimate ("mean" or "median")
-##' @param conf.int (logical) include confidence interval?
-##' @param conf.level probability level for CI
-##' @param conf.method method for computing confidence intervals ("quantile" or "HPDinterval")
-##' @param ... unused
-##' @importFrom coda HPDinterval as.mcmc
-##' @export
+#'
+#' @param x an object of class \sQuote{"stanfit"}
+#' @param pars (character) specification of which parameters to nclude
+#' @param estimate.method method for computing point estimate ("mean" or median")
+#' @param conf.int (logical) include confidence interval?
+#' @param conf.level probability level for CI
+#' @param conf.method method for computing confidence intervals
+#' ("quantile" or "HPDinterval")
+#' @param droppars Parameters not to include in the output (such
+#' as log-probability information)
+#' @param ... unused
+#' 
+#' @name mcmc_tidiers
+#' 
+#' @examples
+#' 
+#' \dontrun{
+#' 
+#' # 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")
+#' 
+#' 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))
+#' 
+#' if (requireNamespace("rstan", quietly = TRUE)) {
+#'   set.seed(2015)
+#'   rstan_example <- stan(file = model_file, data = schools_dat, 
+#'                         iter = 100, chains = 2)
+#' }
+#' 
+#' }
+#' 
+#' if (requireNamespace("rstan", quietly = TRUE)) {
+#'   # the object from the above code was saved as rstan_example.rda
+#'   infile <- system.file("extdata", "rstan_example.rda", package = "broom")
+#'   load(infile)
+#'   
+#'   tidy(rstan_example)
+#'   tidy(rstan_example, conf.int = TRUE)
+#'   
+#'   td_mean <- tidy(rstan_example, conf.int = TRUE)
+#'   td_median <- tidy(rstan_example, conf.int = TRUE, estimate.method = "median")
+#'   
+#'   library(dplyr)
+#'   library(ggplot2)
+#'   tds <- rbind(mutate(td_mean, method = "mean"),
+#'                mutate(td_median, method = "median"))
+#'   
+#'   ggplot(tds, aes(estimate, term)) +
+#'     geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
+#'     geom_point(aes(color = method))
+#' }
+#' 
+#' 
+#' @export
 tidyMCMC <- function(x,
                      pars, ## ?? do other 
-                     pt.method = "mean",
+                     estimate.method = "mean",
                      conf.int = FALSE,
                      conf.level = 0.95,
                      conf.method = "quantile",
-                     droppars="lp__",
+                     droppars = "lp__",
                      ...) {
     ss <- as.matrix(x)  ## works natively on stanfit, mcmc.list, mcmc objects
-    ss <- ss[,!colnames(ss) %in% droppars]  ## drop log-probability info
+    ss <- ss[, !colnames(ss) %in% droppars]  ## drop log-probability info
     if (!missing(pars)) {
-        if (length(badpars <- which(!pars %in% colnames(ss)))>0) {
-            stop("unrecognized parameters: ",pars[badpars])
+        if (length(badpars <- which(!pars %in% colnames(ss))) > 0) {
+            stop("unrecognized parameters: ", pars[badpars])
         }
-        ss <- ss[,pars]
+        ss <- ss[, pars]
     }
-    if (pt.method=="mean") {
-        m <- colMeans(ss)
-    } else if (pt.method=="median") {
-        m <- apply(ss,1,median)
-    } else stop("unknown pt.method ",pt.method)
-    ret <- data.frame(estimate=colMeans(ss),
-                      std.error=apply(ss,2,sd))
+
+    estimate.method <- match.arg(estimate.method, c("mean", "median"))
+    m <- switch(estimate.method,
+                mean = colMeans(ss),
+                median = apply(ss, 2, median))
+
+    ret <- data.frame(estimate = m,
+                      std.error = apply(ss, 2, sd))
     if (conf.int) {
-        levs <- c((1-conf.level)/2,(1+conf.level)/2)
-        if (conf.method=="quantile") {
-            ci <- t(apply(ss,2,quantile,levs))
-        } else if (conf.method=="HPDinterval") {
-            ci <- HPDinterval(as.mcmc(ss),prob=conf.leve)
-        } else stop("unknown conf.method ",conf.method)
-        colnames(ci) <- c("conf.low","conf.high")
-        ret <- data.frame(ret,ci)
+        levs <- c((1 - conf.level) / 2, (1 + conf.level) / 2)
+
+        conf.method <- match.arg(conf.method, c("quantile", "HPDinterval"))
+        ci <- switch(conf.method,
+                     quantile = t(apply(ss, 2, quantile, levs)),
+                     coda::HPDinterval(coda::as.mcmc(ss), prob = conf.level))
+
+        colnames(ci) <- c("conf.low", "conf.high")
+        ret <- data.frame(ret, ci)
     }
     return(fix_data_frame(ret))
 }
 
-##' @rdname tidyMCMC
+
+##' @rdname mcmc_tidiers
 ##' @export
 tidy.rjags <- function(x,
                        pars, ## ?? do other 
-                       pt.method = "mean",
+                       estimate.method = "mean",
                        conf.int = FALSE,
                        conf.level = 0.95,
                        conf.method = "quantile",
                        ...) {
-    tidyMCMC(as.mcmc(x$BUGS),
+    tidyMCMC(coda::as.mcmc(x$BUGS),
              pars,
-             pt.method,conf.int,conf.level,
-             conf.method,droppars="deviance")
+             estimate.method, conf.int, conf.level,
+             conf.method, droppars = "deviance")
 }
 
-##' @rdname tidyMCMC
+##' @rdname mcmc_tidiers
 ##' @export
 tidy.stanfit <- tidyMCMC

inst/extdata/8schools.stan

@@ -0,0 +1,19 @@
+data {
+  int<lower=0> J; // number of schools 
+  real y[J]; // estimated treatment effects
+  real<lower=0> sigma[J]; // s.e. of effect estimates 
+}
+parameters {
+  real mu; 
+  real<lower=0> tau;
+  real eta[J];
+}
+transformed parameters {
+  real theta[J];
+  for (j in 1:J)
+    theta[j] <- mu + tau * eta[j];
+}
+model {
+  eta ~ normal(0, 1);
+  y ~ normal(theta, sigma);
+}

man/mcmc_tidiers.Rd

@@ -0,0 +1,82 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/mcmc_tidiers.R
+\name{mcmc_tidiers}
+\alias{mcmc_tidiers}
+\alias{tidy.rjags}
+\alias{tidy.stanfit}
+\alias{tidyMCMC}
+\title{Tidying methods for MCMC (Stan, JAGS, etc.) fits}
+\usage{
+tidyMCMC(x, pars, estimate.method = "mean", conf.int = FALSE,
+  conf.level = 0.95, conf.method = "quantile", droppars = "lp__", ...)
+
+\method{tidy}{rjags}(x, pars, estimate.method = "mean", conf.int = FALSE,
+  conf.level = 0.95, conf.method = "quantile", ...)
+
+\method{tidy}{stanfit}(x, pars, estimate.method = "mean", conf.int = FALSE,
+  conf.level = 0.95, conf.method = "quantile", droppars = "lp__", ...)
+}
+\arguments{
+\item{x}{an object of class \sQuote{"stanfit"}}
+
+\item{pars}{(character) specification of which parameters to nclude}
+
+\item{estimate.method}{method for computing point estimate ("mean" or median")}
+
+\item{conf.int}{(logical) include confidence interval?}
+
+\item{conf.level}{probability level for CI}
+
+\item{conf.method}{method for computing confidence intervals
+("quantile" or "HPDinterval")}
+
+\item{droppars}{Parameters not to include in the output (such
+as log-probability information)}
+
+\item{...}{unused}
+}
+\description{
+Tidying methods for MCMC (Stan, JAGS, etc.) fits
+}
+\examples{
+\dontrun{
+
+# 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")
+
+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))
+
+if (requireNamespace("rstan", quietly = TRUE)) {
+  set.seed(2015)
+  rstan_example <- stan(file = model_file, data = schools_dat,
+                        iter = 100, chains = 2)
+}
+
+}
+
+if (requireNamespace("rstan", quietly = TRUE)) {
+  # the object from the above code was saved as rstan_example.rda
+  infile <- system.file("extdata", "rstan_example.rda", package = "broom")
+  load(infile)
+
+  tidy(rstan_example)
+  tidy(rstan_example, conf.int = TRUE)
+
+  td_mean <- tidy(rstan_example, conf.int = TRUE)
+  td_median <- tidy(rstan_example, conf.int = TRUE, estimate.method = "median")
+
+  library(dplyr)
+  library(ggplot2)
+  tds <- rbind(mutate(td_mean, method = "mean"),
+               mutate(td_median, method = "median"))
+
+  ggplot(tds, aes(estimate, term)) +
+    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
+    geom_point(aes(color = method))
+}
+}
+

man/summary_tidiers.Rd

@@ -24,6 +24,7 @@ a one-row data frame with columns
   \item{mean}{mean of original vector}
   \item{q3}{value at the third quartile}
   \item{maximum}{largest value in original vector}
+  \item{NAs}{number of NA values (if any)}
 }
 \description{
 Tidy a summary of a vector.
@@ -35,6 +36,9 @@ s
 
 tidy(s)
 glance(s)
+
+v2 <- c(v,NA)
+tidy(summary(v2))
 }
 \seealso{
 \code{\link{summary}}