Nothing
#' Plot for the conflict of evidence parameters w1 and w2
#'
#' Conflict of evidence plot: this plot displays the posterior distribution of the study's weights w1 and w1.
#' These weights indicate potential conflict of evidence of the studies. The weight w1 indicates deviations
#' with respect to the specificity and w2 to the sensitivity.
#'
#'
#' @param m The object generated by metadiag. The model object must be fitted with the options: re = "sm" and split.w = TRUE.
#'
#' @param group An optional argument which is a variable name indicating a group factor.
#' If set, then the plot is colored by groups.
#' @param title The title of the plot.
#'
#' @param group.colors A character vector with two color names.
#'
#' @seealso \code{\link{metadiag}}.
#'
#' @keywords file
#'
#' @examples
#'
#' ## execute analysis
#' \dontrun{
#' data(ep)
#' ep$design = factor(ep$d1,labels = c("prospective", "retrospective"))
#' m.ep <- metadiag(ep, re = "sm", re.model = "SeSp",
#' split.w = TRUE,
#' df.estimate = TRUE)
#'
#' plotw(m.ep)
#' #Relationship between conflict and study design
#' plotw(m.ep, group = "design")
#'
#' }
#'
#' @import ggplot2
#' @import grid
#' @import gridExtra
#' @import R2jags
#' @import rjags
#' @importFrom stats median quantile
#' @export
plotw <- function(m, group=NULL,
title = "Posterior quantiles (25%, 50%, 75%)",
group.colors = c("blue", "red"))
{
re <- m$re
if(re!="sm")stop("This plot function is only for scale mixtures of normals.")
split.w <- m$split.w
if(split.w == FALSE)stop("This plot function is for split.w==TRUE")
# Patch for the "note" no-visible-binding-for-global-variable
w1.post <- m$BUGSoutput$sims.list$w1
w2.post <- m$BUGSoutput$sims.list$w2
nstudies <- dim(w1.post)[2] #for the number of studies
if(is.null(group)){
dat.w <- data.frame(study.name = factor(1:nstudies),
w.median = c(apply(w1.post, 2, median),
apply(w2.post, 2, median)),
w.lower = c(apply(w1.post, 2, quantile, prob = 0.25),
apply(w2.post, 2, quantile, prob = 0.25)),
w.upper = c(apply(w1.post, 2, quantile, prob = 0.75),
apply(w2.post, 2, quantile, prob = 0.75)),
component = gl(2, nstudies, labels = c("TPR", "FPR")))
theplot <- ggplot(dat.w,
aes_string(x = "study.name", y = "w.median", ymin = "w.lower",
ymax = "w.upper")) +
coord_flip() +
geom_pointrange(lwd = 0.8) +
facet_grid(. ~ component) +
geom_hline(yintercept = 1, colour = "black", size = 0.5, lty =2) +
xlab("Study") +
ylab("Weights") +
ggtitle(title)
} else if (!is.null(group))
{
dat.w <- data.frame(study.name = factor(1:nstudies),
w.median = c(apply(w1.post, 2, median),
apply(w2.post, 2, median)),
w.lower = c(apply(w1.post, 2, quantile, prob = 0.25),
apply(w2.post, 2, quantile, prob = 0.25)),
w.upper = c(apply(w1.post, 2, quantile, prob = 0.75),
apply(w2.post, 2, quantile, prob = 0.75)),
component = gl(2, nstudies, labels = c("TPR", "FPR")),
group = m$data[, group])
theplot <- ggplot(dat.w, aes_string(x = "study.name", y = "w.median", ymin = "w.lower",
ymax = "w.upper", color = "group")) +
coord_flip() +
geom_pointrange(lwd =0.8) +
facet_grid(. ~ component) +
geom_hline(yintercept = 1, colour = "black", size = 0.5, lty =2) +
xlab("Study") +
ylab("Weights") +
ggtitle("Posteriors quantiles (25%, 50%, 75%)") +
scale_color_manual(values = group.colors)
}
else{
stop("Argument 'group' has to be a factor!")
}
return(theplot)
}
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.