Nothing
# various high-level plots
# last modified 2020-07-22 by J. Fox
Hist <- function(x, groups, scale=c("frequency", "percent", "density"), xlab=deparse(substitute(x)),
ylab=scale, main="", breaks="Sturges", ...){
xlab # evaluate
scale <- match.arg(scale)
ylab
if (!missing(groups)){
groupsName <- deparse(substitute(groups))
if (!is.factor(groups)){
if (!(is.character(groups) || is.logical(groups)))
warning("groups variable is not a factor, character, or logical")
groups <- as.factor(groups)
}
counts <- table(groups)
if (any(counts == 0)){
levels <- levels(groups)
warning("the following groups are empty: ", paste(levels[counts == 0], collapse=", "))
}
levels <- levels(groups)
hists <- lapply(levels, function(level) if (counts[level] != 0)
hist(x[groups == level], plot=FALSE, breaks=breaks)
else list(breaks=NA))
range.x <- range(unlist(lapply(hists, function(hist) hist$breaks)), na.rm=TRUE)
n.breaks <- max(sapply(hists, function(hist) length(hist$breaks)))
breaks. <- seq(range.x[1], range.x[2], length=n.breaks)
hists <- lapply(levels, function(level) if (counts[level] != 0)
hist(x[groups == level], plot=FALSE, breaks=breaks.)
else list(counts=0, density=0))
names(hists) <- levels
ylim <- if (scale == "frequency"){
max(sapply(hists, function(hist) max(hist$counts)))
}
else if (scale == "density"){
max(sapply(hists, function(hist) max(hist$density)))
}
else {
max.counts <- sapply(hists, function(hist) max(hist$counts))
tot.counts <- sapply(hists, function(hist) sum(hist$counts))
ylims <- tot.counts*(max(max.counts[tot.counts != 0]/tot.counts[tot.counts != 0]))
names(ylims) <- levels
ylims
}
save.par <- par(mfrow=n2mfrow(sum(counts != 0)), oma = c(0, 0, if (main != "") 1.5 else 0, 0))
on.exit(par(save.par))
for (level in levels){
if (counts[level] == 0) next
if (scale != "percent") Hist(x[groups == level], scale=scale, xlab=xlab, ylab=ylab,
main=paste(groupsName, "=", level), breaks=breaks., ylim=c(0, ylim), ...)
else Hist(x[groups == level], scale=scale, xlab=xlab, ylab=ylab,
main=paste(groupsName, "=", level), breaks=breaks., ylim=c(0, ylim[level]), ...)
}
if (main != "") mtext(side = 3, outer = TRUE, main, cex = 1.2)
return(invisible(hists))
}
x <- na.omit(x)
if (scale == "frequency") {
hist <- hist(x, xlab=xlab, ylab=ylab, main=main, breaks=breaks, ...)
}
else if (scale == "density") {
hist <- hist(x, freq=FALSE, xlab=xlab, ylab=ylab, main=main, breaks=breaks, ...)
}
else {
n <- length(x)
hist <- hist(x, axes=FALSE, xlab=xlab, ylab=ylab, main=main, breaks=breaks, ...)
axis(1)
max <- ceiling(10*par("usr")[4]/n)
at <- if (max <= 3) (0:(2*max))/20
else (0:max)/10
axis(2, at=at*n, labels=at*100)
}
box()
abline(h=0)
invisible(hist)
}
indexplot <- function(x, groups, labels=seq_along(x), id.method="y", type="h", id.n=0, ylab, legend="topright", title, col=palette(), ...){
if (is.data.frame(x)) {
if (missing(labels)) labels <- rownames(x)
x <- as.matrix(x)
}
if (!missing(groups)){
if (missing(title)) title <- deparse(substitute(groups))
if (!is.factor(groups)) groups <- as.factor(groups)
groups <- addNA(groups, ifany=TRUE)
grps <- levels(groups)
grps[is.na(grps)] <- "NA"
levels(groups) <- grps
if (length(grps) > length(col)) stop("too few colors to plot groups")
}
else {
grps <- NULL
legend <- FALSE
}
if (is.matrix(x)){
ids <- NULL
mfrow <- par(mfrow=c(ncol(x), 1))
on.exit(par(mfrow))
if (missing(labels)) labels <- 1:nrow(x)
if (is.null(colnames(x))) colnames(x) <- paste0("Var", 1:ncol(x))
for (i in 1:ncol(x)) {
id <- indexplot(x[, i], groups=groups, labels=labels, id.method=id.method, type=type, id.n=id.n,
ylab=if (missing(ylab)) colnames(x)[i] else ylab, legend=legend, title=title, ...)
ids <- union(ids, id)
legend <- FALSE
}
if (is.null(ids) || any(is.na(x))) return(invisible(NULL)) else {
ids <- sort(ids)
names(ids) <- labels[ids]
if (any(is.na(names(ids))) || all(ids == names(ids))) names(ids) <- NULL
return(ids)
}
}
if (missing(ylab)) ylab <- deparse(substitute(x))
plot(x, type=type, col=if (is.null(grps)) col[1] else col[as.numeric(groups)],
ylab=ylab, xlab="Observation Index", ...)
if (!isFALSE(legend)){
legend(legend, title=title, bty="n",
legend=grps, col=palette()[1:length(grps)], lty=1, horiz=TRUE, xpd=TRUE)
}
if (par("usr")[3] <= 0) abline(h=0, col='gray')
ids <- showLabels(seq_along(x), x, labels=labels, method=id.method, n=id.n)
if (is.null(ids)) return(invisible(NULL)) else return(ids)
}
lineplot <- function(x, ..., legend){
xlab <- deparse(substitute(x))
y <- cbind(...)
m <- ncol(y)
legend <- if (missing(legend)) m > 1
if (legend && m > 1) {
mar <- par("mar")
top <- 3.5 + m
old.mar <- par(mar=c(mar[1:2], top, mar[4]))
on.exit(par(old.mar))
}
if (m > 1) matplot(x, y, type="b", lty=1, xlab=xlab, ylab="")
else plot(x, y, type="b", pch=16, xlab=xlab, ylab=colnames(y))
if (legend && ncol(y) > 1){
xpd <- par(xpd=TRUE)
on.exit(par(xpd), add=TRUE)
ncols <- length(palette())
cols <- rep(1:ncols, 1 + m %/% ncols)[1:m]
usr <- par("usr")
legend(usr[1], usr[4] + 1.2*top*strheight("x"),
legend=colnames(y), col=cols, lty=1, pch=as.character(1:m))
}
return(invisible(NULL))
}
plotDistr <- function(x, p, discrete=FALSE, cdf=FALSE, regions=NULL, col="gray",
legend=TRUE, legend.pos="topright", ...){
if (discrete){
if (cdf){
plot(x, p, ..., type="n")
abline(h=0:1, col="gray")
lines(x, p, ..., type="s")
}
else {
plot(x, p, ..., type="h")
points(x, p, pch=16)
abline(h=0, col="gray")
}
}
else{
if (cdf){
plot(x, p, ..., type="n")
abline(h=0:1, col="gray")
lines(x, p, ..., type="l")
}
else{
plot(x, p, ..., type="n")
abline(h=0, col="gray")
lines(x, p, ..., type="l")
}
if (!is.null(regions)){
col <- rep(col, length=length(regions))
for (i in 1:length(regions)){
region <- regions[[i]]
which.xs <- (x >= region[1] & x <= region[2])
xs <- x[which.xs]
ps <- p[which.xs]
xs <- c(xs[1], xs, xs[length(xs)])
ps <- c(0, ps, 0)
polygon(xs, ps, col=col[i])
}
if (legend){
if (length(unique(col)) > 1){
legend(legend.pos, title = if (length(regions) > 1) "Regions" else "Region",
legend=sapply(regions, function(region){
paste(round(region[1], 2), "to", round(region[2], 2))
}),
col=col, pch=15, pt.cex=2.5, inset=0.02)
}
else
{
legend(legend.pos, title = if (length(regions) > 1) "Regions" else "Region",
legend=sapply(regions, function(region){
paste(round(region[1], 2), "to", round(region[2], 2))
}), inset=0.02)
}
}
}
}
return(invisible(NULL))
}
plotMeans <- function(response, factor1, factor2, error.bars = c("se", "sd", "conf.int", "none"),
level=0.95, xlab=deparse(substitute(factor1)), ylab=paste("mean of", deparse(substitute(response))),
legend.lab=deparse(substitute(factor2)),
legend.pos=c("farright", "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center"),
main="Plot of Means",
pch=1:n.levs.2, lty=1:n.levs.2, col=palette(), connect=TRUE, ...){
if (!is.numeric(response)) stop("Argument response must be numeric.")
xlab # force evaluation
ylab
legend.lab
legend.pos <- match.arg(legend.pos)
error.bars <- match.arg(error.bars)
if (!is.factor(factor1)) {
if (!(is.character(factor1) || is.logical(factor1)))
stop("Argument factor1 must be a factor, character, or logical.")
factor1 <- as.factor(factor1)
}
if (missing(factor2)){
valid <- complete.cases(factor1, response)
factor1 <- factor1[valid]
response <- response[valid]
means <- tapply(response, factor1, mean)
sds <- tapply(response, factor1, sd)
ns <- tapply(response, factor1, length)
if (error.bars == "se") sds <- sds/sqrt(ns)
if (error.bars == "conf.int") sds <- qt((1 - level)/2, df=ns - 1, lower.tail=FALSE) * sds/sqrt(ns)
sds[is.na(sds)] <- 0
yrange <- if (error.bars != "none") c( min(means - sds, na.rm=TRUE), max(means + sds, na.rm=TRUE)) else range(means, na.rm=TRUE)
levs <- levels(factor1)
n.levs <- length(levs)
plot(c(1, n.levs), yrange, type="n", xlab=xlab, ylab=ylab, axes=FALSE, main=main, ...)
points(1:n.levs, means, type=if (connect) "b" else "p", pch=16, cex=2)
box()
axis(2)
axis(1, at=1:n.levs, labels=levs)
if (error.bars != "none") arrows(1:n.levs, means - sds, 1:n.levs, means + sds,
angle=90, lty=2, code=3, length=0.125)
}
else {
if (!is.factor(factor2)) {
if (!(is.character(factor2) || is.logical(factor2)))
stop("Argument factor2 must be a factor, charcter, or logical.")
factor2 <- as.factor(factor2)
}
valid <- complete.cases(factor1, factor2, response)
factor1 <- factor1[valid]
factor2 <- factor2[valid]
response <- response[valid]
means <- tapply(response, list(factor1, factor2), mean)
sds <- tapply(response, list(factor1, factor2), sd)
ns <- tapply(response, list(factor1, factor2), length)
if (error.bars == "se") sds <- sds/sqrt(ns)
if (error.bars == "conf.int") sds <- qt((1 - level)/2, df=ns - 1, lower.tail=FALSE) * sds/sqrt(ns)
sds[is.na(sds)] <- 0
yrange <- if (error.bars != "none") c( min(means - sds, na.rm=TRUE), max(means + sds, na.rm=TRUE)) else range(means, na.rm=TRUE)
levs.1 <- levels(factor1)
levs.2 <- levels(factor2)
n.levs.1 <- length(levs.1)
n.levs.2 <- length(levs.2)
if (length(pch) == 1) pch <- rep(pch, n.levs.2)
if (length(col) == 1) col <- rep(col, n.levs.2)
if (length(lty) == 1) lty <- rep(lty, n.levs.2)
expand.x.range <- if (legend.pos == "farright") 1.4 else 1
if (n.levs.2 > length(col)) stop(sprintf("Number of groups for factor2, %d, exceeds number of distinct colours, %d."), n.levs.2, length(col))
plot(c(1, n.levs.1 * expand.x.range), yrange, type="n", xlab=xlab, ylab=ylab, axes=FALSE, main=main, ...)
box()
axis(2)
axis(1, at=1:n.levs.1, labels=levs.1)
for (i in 1:n.levs.2){
points(1:n.levs.1, means[, i], type=if (connect) "b" else "p", pch=pch[i], cex=2, col=col[i], lty=lty[i])
if (error.bars != "none") arrows(1:n.levs.1, means[, i] - sds[, i],
1:n.levs.1, means[, i] + sds[, i], angle=90, code=3, col=col[i], lty=lty[i], length=0.125)
}
if (legend.pos == "farright"){
x.posn <- n.levs.1 * 1.1
y.posn <- sum(c(0.1, 0.9) * par("usr")[c(3,4)])
# text(x.posn, y.posn, legend.lab, adj=c(0, -.5))
legend(x.posn, y.posn, levs.2, pch=pch, col=col, lty=lty, title=legend.lab)
}
else legend(legend.pos, levs.2, pch=pch, col=col, lty=lty, title=legend.lab, inset=0.02)
}
invisible(NULL)
}
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.