Displaying shinystan output in pdf (knitr) about shinystan

[jgabry]

Do have a deadline when you would need this by? Im working on a companion
package to shinystan that will expose the plotting functions so that you
can call them like regular R functions. This should make it easy to do what
you want.

On Wednesday, August 12, 2015, Andy Yao [email protected] wrote:

Hello!

I'm building an R package that implements MCMC sampling, and for the past
few weeks I have been using shinystan to look at the statistical tests /
traceplots / 3D plots on my chains. Super excited that shinystan is now on
CRAN!!

I am in the process of writing my vignette, and I would like make the
results transparent and reproducible. Is there anyway that I could display
shinystan output (say, traceplot or 3D plot, or histogram of rhat),
through knitr, just like how I could display a ggplot2 graph.

Many thanks!

—
Reply to this email directly or view it on GitHub
#81.

[andyyao95]

Hi Jonah!

Again, thanks for your prompt reply! I don't have a firm deadline, but I would prefer to be able to include that in my vignette because I hope to submit this to journal soon (within next few days). That being said, please take your time! I really appreciate your effort in this.

[jgabry]

Hi Andy, now that I think about it, shinystan will let you save some of the plots as ggplot2 objects (the multi parameter plot in Estimates, the autocorrelation plots in Diagnose, and the bivariate, density and histogram plots in Explore). So you could do that and then take the ggplot2 objects, customize them further using any functions in ggplot2 if you want, and then put them in the vignette.

For the rhat histogram you can use this function (where sso is a shinystan object):

rhat_neff_mcse_hist <- function(sso, which = "rhat", fill_color = "#66a7e0",
                             line_color = "#006dcc", line_width = 0.2,
                             binwidth = NULL, theme_elements) {

  base_theme <- theme_classic()
  smry <- sso@summary
  xlab <- switch(which,
                 rhat = "Rhat statistic",
                 n_eff = "Effective sample size",
                 mcse = "Monte Carlo SE / Posterior SD"
  )
  x <- switch(which,
              rhat = smry[, "Rhat"],
              n_eff = smry[, "n_eff"],
              mcse = smry[, "se_mean"] / smry[, "sd"]
              )
  my_labs <- labs(y = "", x = xlab)
  graph <- qplot(x, color = I(line_color), fill = I(fill_color),
                 size = I(line_width), binwidth = I(binwidth))
  graph + 
    my_labs + 
    base_theme +
    if (!missing(theme_elements)) theme_elements
}

The theme_elements argument isn't necessary, but if you want to use it it will allow you to customize the appearance by doing something like this:

my_theme <- theme(
      axis.line.x = element_line(size = 3, color = "#222222"),
      axis.line.y = element_line(size = 0.5),
      axis.title = element_text(face = "bold", size = 13)
      )
rhat_neff_mcse_hist(sso, which = "rhat", theme_elements = my_theme)

This particular my_theme object will more or less give you the axes that shinystan uses.

[jgabry]

Or you can just extract the rhat from sso@summary first and change the function so you're just passing in an rhat vector if you prefer.

[jgabry]

And for the dynamic trace plot that you can zoom and whatnot, you can use this:

dynamic_trace <- function(sso, param_name, chain = 0, # 0 = all chains
                                          stack = FALSE, grid = FALSE) {

  color_vector <- function(n) {
    hues = seq(15, 375, length=n+1)
    hcl(h=hues, l=50, c=50)[1:n]
  }
  param_samps <- sso@samps_all[,, param_name]
  dim_samps <- dim(param_samps)
  nChains <- if (is.null(dim_samps)) 1 else dim_samps[2L]
  if (nChains == 1) {
    param_chains <- xts::as.xts(ts(param_samps, start = 1))
  } else {
    if (chain != 0) {
      param_samps <- param_samps[, chain]
      param_chains <- xts::as.xts(ts(param_samps, start = 1))
    } else {
      param_chains <- xts::as.xts(ts(param_samps[, 1L], start = 1))
      for (i in 2:nChains) {
        param_chains <- cbind(param_chains, xts::as.xts(ts(param_samps[,i], start = 1)))
      }
      colnames(param_chains) <- paste0("Chain", 1:nChains)
    }
  }

  `%>%` <- dygraphs::`%>%`
  y_axis_label_remove <- if (stack) "white" else NULL
  clrs <- color_vector(nChains)
  if (chain != 0) clrs <- clrs[chain]
  dygraphs::dygraph(param_chains, xlab = "Iteration", ylab = param_name) %>%
    dygraphs::dyAxis("y", axisLabelColor = y_axis_label_remove, pixelsPerLabel = 30) %>%
    dygraphs::dyAxis("x", axisLabelColor = "white") %>%
    dygraphs::dyOptions(colors = clrs, stackedGraph = stack, drawGrid = grid,
                    animatedZooms = TRUE, axisLineColor = "#222222") %>%
    dygraphs::dyRangeSelector(height = 15) %>%
    dygraphs::dyLegend(show = "never") %>%
    dygraphs::dyHighlight(highlightCircleSize = 4,
                          highlightSeriesBackgroundAlpha = 1/3,
                          hideOnMouseOut = TRUE,
                          highlightSeriesOpts = list(strokeWidth = 1.75))
}

[bob-carpenter]

I found ggsave() a lifesafer for producing pdfs from ggplot.
I'm only saying this because I spent several years just calling
pdf() and then plot() and then dev.off(), which is so error prone
and doesn't work as well with scaling.

• Bob

On Aug 12, 2015, at 5:21 PM, Jonah Gabry [email protected] wrote:

And for the dynamic trace plot that you can zoom and whatnot, you can use this:

dynamic_trace <- function(sso, param_name, chain = 0, # 0 = all chains
stack = FALSE, grid = FALSE) {

color_vector <- function(n) {
hues = seq(15, 375, length=n+1)
hcl(h=hues, l=50, c=50)[1:n]
}
param_samps <- sso@samps_all[,, param_name]
dim_samps <- dim(param_samps)
nChains <- if (is.null(dim_samps)) 1 else dim_samps[2L]
if (nChains == 1) {
param_chains <- xts::as.xts(ts(param_samps, start = 1))
} else {
if (chain != 0) {
param_samps <- param_samps[, chain]
param_chains <- xts::as.xts(ts(param_samps, start = 1))
} else {
param_chains <- xts::as.xts(ts(param_samps[, 1L], start = 1))
for (i in 2:nChains) {
param_chains <- cbind(param_chains, xts::as.xts(ts(param_samps[,i], start = 1)))
}
colnames(param_chains) <- paste0("Chain", 1:nChains)
}
}

%>% <- dygraphs::%>%
y_axis_label_remove <- if (stack) "white" else NULL
clrs <- color_vector(nChains)
if (chain != 0) clrs <- clrs[chain]
dygraphs::dygraph(param_chains, xlab = "Iteration", ylab = param_name) %>%
dygraphs::dyAxis("y", axisLabelColor = y_axis_label_remove, pixelsPerLabel = 30) %>%
dygraphs::dyAxis("x", axisLabelColor = "white") %>%
dygraphs::dyOptions(colors = clrs, stackedGraph = stack, drawGrid = grid,
animatedZooms = TRUE, axisLineColor = "#222222") %>%
dygraphs::dyRangeSelector(height = 15) %>%
dygraphs::dyLegend(show = "never") %>%
dygraphs::dyHighlight(highlightCircleSize = 4,
highlightSeriesBackgroundAlpha = 1/3,
hideOnMouseOut = TRUE,
highlightSeriesOpts = list(strokeWidth = 1.75))
}

—
Reply to this email directly or view it on GitHub.

[jgabry]

Yeah, ggsave is really great. That's what's used if you save a plot as pdf
in shinystan.

On Wednesday, August 12, 2015, Bob Carpenter [email protected]
wrote:

I found ggsave() a lifesafer for producing pdfs from ggplot.
I'm only saying this because I spent several years just calling
pdf() and then plot() and then dev.off(), which is so error prone
and doesn't work as well with scaling.

• Bob

On Aug 12, 2015, at 5:21 PM, Jonah Gabry <[email protected]
javascript:_e(%7B%7D,'cvml','[email protected]');> wrote:

And for the dynamic trace plot that you can zoom and whatnot, you can
use this:

dynamic_trace <- function(sso, param_name, chain = 0, # 0 = all chains
stack = FALSE, grid = FALSE) {

color_vector <- function(n) {
hues = seq(15, 375, length=n+1)
hcl(h=hues, l=50, c=50)[1:n]
}
param_samps <- sso@samps_all[,, param_name]
dim_samps <- dim(param_samps)
nChains <- if (is.null(dim_samps)) 1 else dim_samps[2L]
if (nChains == 1) {
param_chains <- xts::as.xts(ts(param_samps, start = 1))
} else {
if (chain != 0) {
param_samps <- param_samps[, chain]
param_chains <- xts::as.xts(ts(param_samps, start = 1))
} else {
param_chains <- xts::as.xts(ts(param_samps[, 1L], start = 1))
for (i in 2:nChains) {
param_chains <- cbind(param_chains, xts::as.xts(ts(param_samps[,i],
start = 1)))
}
colnames(param_chains) <- paste0("Chain", 1:nChains)
}
}

%>% <- dygraphs::%>%
y_axis_label_remove <- if (stack) "white" else NULL
clrs <- color_vector(nChains)
if (chain != 0) clrs <- clrs[chain]
dygraphs::dygraph(param_chains, xlab = "Iteration", ylab = param_name)
%>%
dygraphs::dyAxis("y", axisLabelColor = y_axis_label_remove,
pixelsPerLabel = 30) %>%
dygraphs::dyAxis("x", axisLabelColor = "white") %>%
dygraphs::dyOptions(colors = clrs, stackedGraph = stack, drawGrid = grid,
animatedZooms = TRUE, axisLineColor = "#222222") %>%
dygraphs::dyRangeSelector(height = 15) %>%
dygraphs::dyLegend(show = "never") %>%
dygraphs::dyHighlight(highlightCircleSize = 4,
highlightSeriesBackgroundAlpha = 1/3,
hideOnMouseOut = TRUE,
highlightSeriesOpts = list(strokeWidth = 1.75))
}

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#81 (comment).

[andyyao95]

@jgabry this is so helpful! Thanks a lot for the functions!