Nothing
#' Generic diagnostic function.
#' @param object The object generated by the function hmr.
#' @param ... \dots
#'
#'
#' @export
diagnostic = function(object, ...) UseMethod("diagnostic", object)
#' @title Diagnostic function for hmr object in jarbes
#'
#' @description This function performers a specially designed diagnostic for a hmr object
#'
#'
#' @param object The object generated by the function hmr.
#' @param median.w Change color if median of a weight > median.w. The default value is 1.5.
#' @param study.names Character vector containing names of the studies used.
#' @param size.forest Size of the center symbol mark in the forest-plot lines
#' @param lwd.forest Thickness of the lines in the forest-plot
#' @param shape.forest Type of symbol for the center mark in the forest-plot lines
#' @param mu.phi Also plot the distribution of mu.phi. Default value is TRUE.
#' @param mu.phi.x.lim.low Lower limit of the prior to posterior plot for mu.phi
#' @param mu.phi.x.lim.up Upper limit of the prior to posterior plot for mu.phi
#' @param ... \dots
#'
#'
#' @import ggplot2
#' @import gridExtra
#'
#' @export
#'
diagnostic.hmr = function(object,
# Forest plot
median.w = 1.5,
study.names,
size.forest = 0.4, # Point size
lwd.forest = 0.2, # Thickness of the lines
shape.forest = 23,
mu.phi = TRUE,
mu.phi.x.lim.low = -10,
mu.phi.x.lim.up = 10, ...) {
if (object$re != "sm") stop("This plot function is only for scale mixtures of normals.")
x=y=ylo=yhi=NULL
# Posteriors for the studies' weights
# Forest plot for weights ...
if (object$split.w == FALSE) {
w = object$BUGSoutput$sims.list$w
w.s = apply(w, 2, median)
w.u = apply(w, 2, quantile, prob = 0.75)
w.l = apply(w, 2, quantile, prob = 0.25)
n.studies = length(w.s)
w.col = ifelse(w.s<median.w, "red", "blue")
w.col[length(w.s)] = "green"
if (missing(study.names)) {
study.names = 1:(n.studies-1)
}
study.names = c(as.character(study.names), "Cohort (individual data)")
dat.weights = data.frame(x = study.names,
y = w.s,
ylo = w.l,
yhi = w.u)
} else {
w1 = object$BUGSoutput$sims.list$w1
w2 = object$BUGSoutput$sims.list$w2
w1.s = apply(w1, 2, median)
w2.s = apply(w2, 2, median)
w1.u = apply(w1, 2, quantile, prob = 0.75)
w2.u = apply(w2, 2, quantile, prob = 0.75)
w1.l = apply(w1, 2, quantile, prob = 0.25)
w2.l = apply(w2, 2, quantile, prob = 0.25)
n.studies = length(w1.s)
w1.col = ifelse(w1.s<median.w, "red", "blue")
w1.col[n.studies] = "green"
w2.col = ifelse(w2.s<median.w, "red", "blue")
w2.col[n.studies] = "green"
w.col = c(w1.col, w2.col)
if (missing(study.names)) {
study.names = 1:(n.studies-1)
}
study.names = c(as.character(study.names), "Cohort (individual data)")
dat.weights = data.frame(x = as.character(study.names),
y = c(w1.s, w2.s),
ylo = c(w1.l, w2.l),
yhi = c(w1.u, w2.u),
component =
gl(2, n.studies, labels = c("Baseline Risk",
"Treatment Effect")))
}
diagnostic.plot2 = ggplot(dat.weights,
aes(x = x,
y = y,
ymin = ylo,
ymax = yhi,
size = size.forest # Point size
)) +
geom_pointrange(colour = w.col,
lwd = lwd.forest, # Thickness of the lines
shape = shape.forest) +
coord_flip() +
geom_hline(yintercept = 1, lty = 2) +
xlab("Study") +
ylab("Outliers detection weights") +
ggtitle("Weights") +
theme_bw()
if (object$split.w == TRUE)
diagnostic.plot2 = diagnostic.plot2 + facet_grid(. ~ component)
if (mu.phi == TRUE) {
mu.phi = object$BUGSoutput$sims.list$mu.phi
mean.mu.phi = object$prior$mean.mu.phi
sd.mu.phi = object$prior$sd.mu.phi
# Posterior distribution of mu.phi
df.mu.phi = data.frame(x = mu.phi)
mycolors <- c("Prior" = "red", "Posterior" = "blue")
diagnostic.plot1 = ggplot(df.mu.phi,
aes(x=x, colour = "Posterior"))+
geom_density()+
geom_histogram(aes(y = stat(density)), fill = "royalblue",
alpha = 0.3, bins=60) +
geom_vline(xintercept = mean.mu.phi, colour = "red", size = 1, lty = 2) +
geom_vline(xintercept = mean(mu.phi),
colour = "blue", size = 1, lty = 2) +
stat_function(fun = dlogis,
n = 101,
args = list(location = mean.mu.phi, scale = sd.mu.phi),
size = 1, colour = "red") +
scale_color_manual(values = mycolors)+
labs(x=expression(mu[phi]), y="Density", color = "Density") +
ggtitle("Prior to Posterior")+
xlim(c(mu.phi.x.lim.low, mu.phi.x.lim.up)) +
theme_bw()
# densities mu AG vs ID...
return(grid.arrange(diagnostic.plot1, diagnostic.plot2, ncol = 2, nrow = 1))
} else {
return(diagnostic.plot2)
}
}
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.