Nothing
#'plot ex-situ prior
#'
#'\code{plotExPrior} plot informative and non-informative priors
#'
#'@param resExPrior output from the genExPrior function
#'@param plotExData boolean asking whether to additionally plot the ex-situ data
#'@return a plot
#'@examples
#'\donttest{
#'exdata <- data.frame(val=c(c(2,3,4),c(2,1),c(6,7,2,3)),
#' site_id=c(rep("a",3),rep("b",2),rep("c",4)),
#' x = c(c(2,3,4),c(2,3),c(2,2,3,3)),
#' y = c(c(2,2,3),c(3,2),c(2,3,2,3)))
#'ex_prior <- genExPrior(exdata=exdata,theta=seq(from=-5,to=5,by=1))
#'plotExPrior(ex_prior)
#'}
#'@import ggplot2
#'@import grid
#'@importFrom gtable gtable_add_grob gtable_add_cols gtable_add_rows
#'@export
plotExPrior <- function(resExPrior,plotExData=F){
# construct single dataframe with uninformative and informative priors
# easier for ggplot2
df_exPrior <- data.frame(theta = resExPrior$uPrior$x,
uPrior = resExPrior$uPrior$y,
exPrior = resExPrior$exPrior$y)
df_exPrior <- reshape2::melt(df_exPrior, id.vars = "theta")
# Construct the plots
# =============================================================================
# To plot both pdf and histogram of measurements in the same plot with dual y-axes,
# we need to:
# 1. Plot each plot in separate panel
# 2. Combine two plots together
# =============================================================================
# first define the range for the x-axis,
# so that it is similar in both overlaid plots
xrange <- range(df_exPrior$theta)
# exPrior pdf
p1 <-
ggplot2::ggplot(df_exPrior, ggplot2::aes(theta, value, color = variable)) +
ggplot2::geom_line() +
ggplot2::scale_colour_manual(values=c('black','blue'), labels=c(expression(paste('p(', theta, ')') ),
expression(paste('p(', theta, '|y)')))) +
labs(x = expression(theta), y = NULL) +
ggplot2::scale_x_continuous(limits = xrange) +
ggplot2::theme_bw() +
ggplot2::theme(axis.text.y = ggplot2::element_text(colour="#68382C", size = 14),
axis.text.x = ggplot2::element_text(size = 13),
axis.title = ggplot2::element_text(size = 14, face = "bold"),
text = ggplot2::element_text(),
legend.title = ggplot2::element_blank(),
legend.background = ggplot2::element_rect(colour = "black"))
# exPrior histogram
p2 <- ggplot2::ggplot(resExPrior$exdata, ggplot2::aes(val, fill = site_id)) +
ggplot2::geom_histogram(alpha = .5) +
ggplot2::scale_fill_discrete("Site") +
labs(x = 'y', y = NULL) +
# set x axis same as x axis from p1
ggplot2::scale_x_continuous(limits = xrange) +
ggplot2::theme_bw() +
ggplot2::theme(axis.text.y = ggplot2::element_text(colour="#00A4E6", size=14),
axis.text.x = ggplot2::element_text(size = 13),
panel.background = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
panel.grid.major = ggplot2::element_blank(),
text = ggplot2::element_text(),
legend.background = ggplot2::element_rect(colour = "black"))
# combine two plots
# =======================================================================
# To combine two plots, we need to access the inner working of ggplots object.
# Use gtable package to make gtable objects from ggplot objects
# gtable object shows how grobs(grid graphical object) are put together to form a ggplot
# =======================================================================
# Get the plot grobs
g1 <- ggplot2::ggplotGrob(p1)
g2 <- ggplot2::ggplotGrob(p2)
# Get the locations of the plot panels in g1.
pp <- c(subset(g1$layout, name == "panel", se = t:r))
# Overlap panel for second plot on that of the first plot
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, pp$l, pp$b, pp$l)
# ggplot contains many labels that are themselves complex grob;
# usually a text grob surrounded by margins.
# When moving the grobs from, say, the left to the right of a plot,
# make sure the margins and the justifications are swapped around.
# The function below does the swapping.
# Taken from the cowplot package:
# https://github.com/wilkelab/cowplot/blob/master/R/switch_axis.R
hinvert_title_grob <- function(grob){
# Swap the widths
widths <- grob$widths
grob$widths[1] <- widths[3]
grob$widths[3] <- widths[1]
grob$vp[[1]]$layout$widths[1] <- widths[3]
grob$vp[[1]]$layout$widths[3] <- widths[1]
# Fix the justification
grob$children[[1]]$hjust <- 1 - grob$children[[1]]$hjust
grob$children[[1]]$vjust <- 1 - grob$children[[1]]$vjust
grob$children[[1]]$x <- ggplot2::unit(1, "npc") - grob$children[[1]]$x
grob
}
# Get the y axis from g2 (axis line, tick marks, and tick mark labels)
index <- which(g2$layout$name == "axis-l") # Which grob
yaxis <- g2$grobs[[index]] # Extract the grob
# yaxis is a complex of grobs containing the axis line, the tick marks, and the tick mark labels.
# The relevant grobs are contained in axis$children:
# axis$children[[1]] contains the axis line;
# axis$children[[2]] contains the tick marks and tick mark labels.
# Second, swap tick marks and tick mark labels
ticks <- yaxis$children[[2]]
ticks$widths <- rev(ticks$widths)
ticks$grobs <- rev(ticks$grobs)
# Third, move the tick marks
# Tick mark lengths can change.
# A function to get the original tick mark length
# Taken from the cowplot package:
# https://github.com/wilkelab/cowplot/blob/master/R/switch_axis.R
plot_theme <- function(p) {
plyr::defaults(p$theme, ggplot2::theme_get())
}
tml <- plot_theme(p1)$axis.ticks.length # Tick mark length
ticks$grobs[[1]]$x <- ticks$grobs[[1]]$x - ggplot2::unit(1, "npc") + tml
# Fourth, swap margins and fix justifications for the tick mark labels
ticks$grobs[[2]] <- hinvert_title_grob(ticks$grobs[[2]])
# Fifth, put ticks back into yaxis
yaxis$children[[2]] <- ticks
# Put the transformed yaxis on the right side of g
g <- gtable_add_cols(g, g2$widths[g2$layout[index, ]$l], pp$r)
g <- gtable_add_grob(g, yaxis, pp$t, pp$r + 1, pp$b, pp$r + 1, clip = "off",
name = "axis-r")
# Labels grob
left = textGrob(expression(p), x = 0.02, y = 0.5, just = c("left", "top"),
gp = gpar(fontsize = 14, col = "#68382C"))
right = textGrob("Count", x = 1.0, y = 0.5, just = c("right", "top"),
gp = gpar(fontsize = 14, col = "#00a4e6"))
labs = gTree("Labs", children = gList(left, right))
# New row in the gtable for labels
height = ggplot2::unit(3, "grobheight", left)
g = gtable_add_rows(g, height, 2)
# Put the label in the new row
g = gtable_add_grob(g, labs, t=3, l=3, r=5)
# Turn off clipping in the plot panel
g$layout[which(g$layout$name == "panel"), ]$clip = "off"
# ======================================================================
# Everything looks great except for the missing panel in the 2nd plot.
# We'll try to fix it by:
# 1. Find out grobs information of the legend
# 2. Extract it from 2nd plot and add it to combined plot
# ======================================================================
# Extract legend http://stackoverflow.com/questions/35673831/position-a-legend-with-gtable
leg1 <- g1$grobs[[which(g1$layout$name == "guide-box")]]
leg2 <- g2$grobs[[which(g2$layout$name == "guide-box")]]
# row binding two legends so that we can have two legends lay out vertically
# then replace the legend in the combined plot to display
#g$grobs[[which(g$layout$name == "guide-box")]] <- gtable:::rbind_gtable(leg1, leg2, "first")
# Adjust height and width of panel to fit legends
g$heights[[6]] = ggplot2::unit(0.7, "cm")
g$widths[[6]] = ggplot2::unit(3.1,"cm")
if(!plotExData){ # Plot pdf only
p1 + ggplot2::theme_bw()
}else{
grid.draw(g) + ggplot2::theme_bw() # Plot both pdf and histograms
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.