Nothing
##############################################################################
# Copyright (c) 2017 Russell V. Lenth #
# #
# This file is part of the unrepx package for R (*unrepx*) #
# #
# *unrepx* is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published by #
# the Free Software Foundation, either version 2 of the License, or #
# (at your option) any later version. #
# #
# *unrepx* is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with R and *unrepx*. If not, see #
# <https://www.r-project.org/Licenses/> and/or #
# <http://www.gnu.org/licenses/>. #
##############################################################################
### Plots for visual identification of effects
# Half-normal plot (or normal plot)
# if col = TRUE, we color pos blue, neg red, zero black. if FALSE, use black. Else use colors in col
hnplot = function(effects, ref = TRUE, half = TRUE, horiz = TRUE, method = "Zahn",
a = .375, col = half, pch = 16, ID = FALSE, alpha, ...) {
n.effects = length(effects)
abseff = abs(effects)
q = (rank(abseff) - a) / (n.effects + 1 - 2*a)
hns = qnorm((1 + q) / 2)
if (is.logical(col)) {
if (col)
col = sapply(zapsmall(effects),
function(x) ifelse(x < 0, "red", ifelse(x > 0, "blue", "black")))
else
col = "black"
}
if (half) {
eff = abseff
elab = "Absolute effects"
score = hns
slab = "Half-normal quantiles"
}
else {
eff = effects
elab = "Effects"
score = qnorm((rank(eff) - a) / (length(eff) + 1 - 2*a))
slab = "Normal quantiles"
}
if (ref)
slope = PSE(effects, method)
if (!missing(alpha)) {
refdist = .getrefdist(n.effects, method, save = FALSE)
MEs = slope * c(quantile(refdist$abst, 1 - alpha), quantile(refdist$max.abst, 1 - alpha))
MElabs = c("ME", "SME")
if (half)
efflim = range(c(eff, MEs[1]))
else {
efflim = range(c(-MEs[1], eff, MEs[1]))
MEs = c(-rev(MEs), MEs)
MElabs = c("-SME", "-ME", "ME", "SME")
}
}
else
efflim = range(eff)
if (is.numeric(ID)) {
thresh = ID[1]
ID = TRUE
}
else
thresh = Inf
if(ID) {
if (is.null(names(abseff)))
names(abseff) = seq_along(abseff)
idx = which(abseff >= thresh) # indices to auto-ID
if (!is.infinite(thresh) && (length(idx) == 0))
ID = FALSE
}
if(horiz) {
plot(score ~ eff, xlim = efflim, ylab = slab, xlab = elab, col = col, pch = pch, ...)
if(ref)
abline(0, 1 / slope)
if (!missing(alpha)) {
sapply(MEs, function(xx) abline(v = xx, lty = 2))
axis(3, at = MEs, labels = MElabs, cex.axis = .75)
}
if (ID) {
if (is.infinite(thresh))
identify(eff, score, labels = names(abseff))
else
text(eff[idx], score[idx], names(abseff[idx]),
pos = sapply(sign(eff[idx]), function(s) ifelse(s>0, 2, 4)))
}
}
else {
plot(eff ~ score, ylim = efflim, xlab = slab, ylab = elab, col = col, pch = pch, ...)
if(ref)
abline(0, slope)
if (!missing(alpha)) {
sapply(MEs, function(xx) abline(h = xx, lty = 2))
axis(4, at = MEs, labels = MElabs, cex.axis = .75)
}
if (ID) {
if (is.infinite(thresh))
identify(score, eff, labels = names(abseff))
else
text(score[idx], eff[idx], names(abseff[idx]),
pos = sapply(sign(eff[idx]), function(s) ifelse(s>0, 1, 3)))
}
}
if (!half && ref) {
abline(h = 0, lty = 3, col = "darkgray")
abline(v = 0, lty = 3, col = "darkgray")
}
if (ref)
title(sub = paste0("Reference SD = ", signif(slope, 4), " (", method, " method)"),
cex.sub=.75, adj = 1)
if(!missing(alpha)) {
MEs = rev(MEs)
title(sub = paste("ME =", signif(MEs[2], 3)," SME =", signif(MEs[1], 3)),
cex.sub=.75, adj = 0)
}
invisible(list(eff = eff, abseff = abseff, hnscore = hns))
}
### Version with automatic bins
dot.plot = function(x, pch = 16, cex.dot = 1, spacing = 1, xlab, xlim = range(x), ...) {
if(missing(xlab))
xlab = as.character(substitute(x))
# max(x) - min(x)]
xspan = diff(range(x))
# determine inverse ranks
idx = seq_along(x)
idx[order(x)] = idx
# make a blank plot
plot(x, rep(0, length(x)), type = "n", axes = FALSE, xlab = xlab, ylab = "", xlim = xlim)
# draw scale
axis(1)
ylow = par("usr")[3]
ytop = par("usr")[4]
abline(h = ylow) # extend to full width
# Just in case we ID them later
if(is.null(names(x)))
names(x) = seq_along(x)
env = new.env()
env$orig.x = x # used by dot.id()
env$cex.dot = cex.dot
RG = if(dev.interactive()) recordGraphics
else function(expr, list, env) expr
# draw points and support resizing in both directions
RG({
cxy = par("cxy") * cex.dot
bins = as.integer(1.25 * xspan / cxy[1])
xinc = diff(pretty(x, n = bins)[1:2])
rx = xinc * round(x / xinc, 0)
freq = table(rx)
xx = rep(as.numeric(names(freq)), freq) [idx]
yy = unlist(lapply(freq, seq_len)) [idx]
yinc = env$yinc = 0.5 * spacing * cxy[2]
yyv = ylow + yinc * (yy - .5)
points(xx, yyv, pch = pch, cex = env$cex.dot, ...)
clip = which(yyv > ytop)
if(length(clip) > 0)
text(xx[clip], rep(ytop + .5 * yinc, length(clip)), "(more)", col = "red", cex=.5, xpd = TRUE)
# Additional labels user may add programatically later
# by defining x.id, height, cex.id, col.id in env
if(!is.null(env$x.id)) {
rxl = xinc * round(env$x.id / xinc)
for (xl in rxl) {
labs = names(rx[rx == xl])
yl = ylow + cxy[2] * (env$height.id - 1 + env$cex.id * seq_along(labs))
text(xl, yl, labs, col = env$col.id, cex = env$cex.id)
}
}
},
list(), env)
invisible(env)
}
# "identify" function for dot plots - specify environment returned by dot.plot
# with modify = TRUE, previously identified points are used but user can specify different height etc.
dot.id = function(env, height.id = 2, cex.id = 1, col.id = "black") {
if (!dev.interactive())
stop("Graphics device is not interactive")
ylow = par("usr")[3]
x = env$orig.x
pts = identify(x, rep(ylow, length(x)), pos = FALSE, plot = FALSE)
if (length(pts) > 0) {
modify = TRUE
env$x.id = env$orig.x[pts]
env$height.id = height.id
env$cex.id = cex.id
env$col.id = col.id
}
message("You may need to refresh or slightly resize the plot")
invisible(env$x.id)
}
dot.mod = function(env, ...) {
if (!dev.interactive())
stop("Graphics device is not interactive")
dots = list(...)
names(dots) = sapply(names(dots), function(nm) match.arg(nm, c("cex.dot", "height.id", "col.id", "cex.id")))
for (nm in names(dots))
env[[nm]] = dots[[nm]]
message("You may need to refresh or slightly resize the plot")
}
# Reference plot
refplot = function(effects, ref = TRUE, half = TRUE, method = "Zahn",
col = half, guides = FALSE, ID = FALSE, pch = 16, xlab, xlim, ...) {
n.effects = length(effects)
if(missing(xlab))
xlab = ifelse(half, "Absolute effects", "Effects")
ae = abs(effects)
if (is.logical(col)) {
if (col)
col = sapply(zapsmall(effects),
function(x) ifelse(x < 0, "red",
ifelse(x > 0, "blue", "black")))
else
col = "black"
}
if (half)
effects = abs(effects)
xlim = range(c(0, effects))
env = dot.plot(effects, xlim = xlim, xlab = xlab, pch = pch, col = col, ...)
if(is.character(ref)) {
curvetype = match.arg(ref, c("normal", "simulated"))
ref = TRUE
guides = FALSE
}
else
curvetype = "normal"
ytop = ifelse(guides, 1.23, 1)
ylow = par("usr")[3]
yscal = (par("usr")[4] - ylow) / ytop
if(ref) {
stderr = PSE(effects, method)
x = seq(-3 * stderr, 3 * stderr, len = 81)
if (half)
x = x[41:81]
if (curvetype == "normal") {
lines(x, ylow + yscal * exp(-.5*(x / stderr)^2))
lines(c(0, 0), ylow + yscal * c(0,1), lty = 2)
title(sub = paste0("Reference curve SD = ", signif(stderr, 4), " (", method, " method)"),
cex.sub=.75, adj = 1)
}
else { # curvetype = "simulated"
refdist = .getrefdist(n.effects, method)
if(half)
DE = stats::density(stderr * refdist$abst, bw = "SJ")
else
DE = stats::density(stderr * c(refdist$abst, -refdist$abst), bw = "SJ")
lines(DE$x, ylow + yscal * (DE$y / max(DE$y)))
title(sub = paste("Ref dist:", method, "method for", n.effects, "effects"),
cex.sub=.75, adj = 1)
}
if (guides) {
e12 = exp(-.5)
lines(c(-2*stderr, 0, 2*stderr), ylow + yscal * c(0, 2*e12, 0), lty = 3)
points(c(-stderr, stderr), ylow + yscal * c(e12, e12), cex = .5)
lines(c(0, 0), ylow + yscal * c(1, ytop), lty = 3)
}
}
if (is.numeric(ID)) {
thresh = ID[1]
ID = TRUE
}
else
thresh = Inf
if (ID) {
if (is.infinite(thresh))
dot.id(env)
else {
env$x.id = effects[abs(effects) >= thresh]
env$height.id = 2
env$cex.id = 1
env$col.id = "black"
}
}
if(!dev.interactive())
{ # show the text since we won't be re-drawing
text(env$x.id, ylow + yscal * env$height.id *env$yinc, names(env$x.id))
}
invisible(env)
}
# pareto plot
parplot = function (effects, pareto = TRUE, absolute = TRUE, horiz = FALSE,
col = absolute, critvals, method = "Zahn",
alpha = .05, refdist, sim.opts, ###save = TRUE, nsets,
ylab = "Estimated effects",
top = n.effects, cex.annot = .75, ...)
{
n.effects = length(effects)
if (is.null(names(effects)))
names(effects) = seq_along(effects)
if (is.logical(col)) {
if (col)
col = sapply(effects,
function(x) ifelse(x < 0, "pink", "lightblue"))
else
col = rep("lightgray", n.effects)
}
if (missing(critvals)) {
if(missing(refdist))
refdist = .getrefdist(n.effects, method, opts = sim.opts)
critvals = c(quantile(refdist$abst, 1 - alpha),
quantile(refdist$max.abst, 1 - alpha))
critmeth = method
pse = PSE(effects, method)
}
else {
critvals = c(critvals, NA) # SME critval is NA if not provided
critmeth = "Manual"
pse = 1
}
cutoffs = pse * c(-critvals[2], -critvals[1], 0,
critvals[1], critvals[2])
cutlabs = c("-SME", "-ME", "0", "ME", "SME")
names(cutoffs) = cutlabs
if (pareto || (top < n.effects)) {
ord = order(abs(effects))
effects = effects[rev(ord)][seq_len(top)]
col = col[rev(ord)][seq_len(top)]
}
ylim = 1.05 * range(c(effects, cutoffs[c(2, 4)]))
if (absolute) {
effects = abs(effects)
ylim[2] = max(abs(ylim))
ylim[1] = 0
}
if (!horiz) {
barplot(effects, ylim = ylim, ylab = ylab, col = col,
...)
abline(h = 0)
axis(side = 4, at = cutoffs, labels = cutlabs)
sapply(cutoffs[c(1, 2, 4, 5)], function(x) abline(h = x,
lty = 2))
axis(side = 2, at = 2 * ylim, labels = c("", ""))
}
else {
barplot(rev(effects), horiz = TRUE, xlim = ylim, xlab = ylab,
col = rev(col), ...)
abline(v = 0)
axis(side = 3, at = cutoffs, labels = cutlabs)
sapply(cutoffs[c(1, 2, 4, 5)], function(x) abline(v = x,
lty = 2))
axis(side = 1, at = 2 * ylim, labels = c("", ""))
}
melab = paste0(critmeth, " critical values: ME = ", signif(cutoffs[4], 3),
", SME = ", signif(cutoffs[5], 3))
# function to place annotations
annot = function(annot, adj, col = "black") {
title(sub = annot, adj = adj, cex.sub = cex.annot, col.sub = col)
}
if(cex.annot > 0)
annot(melab, 0)
if(top < n.effects) {
if(cex.annot <= 0) cex.annot = .75
annot(paste("Only the largest", top, "effects shown"), 1, "blue")
}
invisible(cutoffs[4:5])
}
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.