# violinPlot.R -- various attempts at violin plots
# try 1: wrapper around the ggplot geom_violin plot function
`violinPlot` <- function( df, aes, ..., horiz=FALSE, boxwid=NULL, facet=NULL,
main="", xlab="X", ylab="Y", log="") {
require( ggplot2)
p <- ggplot( df, aes) + geom_violin( ...)
if (horiz) p <- p + coord_flip()
if ( ! is.null( facet)) {
p <- p + facet_wrap( paste( "~", as.character(facet)), scales="free_x")
}
# add boxplot?
if ( ! is.null(boxwid)) {
p <- p + geom_boxplot( width=boxwid)
}
# any X,Y log scaling?
if ( regexpr( 'Y', toupper(log)) > 0) p <- p + scale_y_log10()
if ( regexpr( 'X', toupper(log)) > 0) p <- p + scale_x_log10()
p <- p + labs( title=main, x=xlab, y=ylab)
p
}
# try 2:: a more generic version, starting from the code at UsingR, by John Verzani <verzani@math.csi.cuny.edu>
# original courtesy of boxplot from base
`violinplot` <- function(x, ...) UseMethod("violinplot", x)
`violinplot.default` <- function(x, orientation=c("vertical","horizontal"), at=NULL, viowid=0.8,
boxwid=0.2, names=NULL, pars=NULL, col=NULL, border=par('fg'),
n.density=128, cut.density=0, bw.density="nrd0",
xlim=NULL, ylim=NULL, log="", main=NULL, sub=NULL,
xlab=NULL, ylab=NULL, ann=par("ann"), axes=TRUE, frame.plot=axes,
pch=NULL, lwd=1, lty=1, ...) {
# deduce which args get passed where
args <- list(x, ...)
namedargs <- if ( !is.null( attributes(args)$names)) (attributes(args)$names != "")
else rep(FALSE, length = length(args))
pars <- c( args[namedargs], pars)
# turn what we were given into the groups to plot
groups <- if ( is.list(x)) x else if ( is.matrix(x)) as.list(as.data.frame(x))
else if ( is.data.frame(x)) as.list(x) else args[!namedargs]
if ( ! (nGroups <- length(groups))) stop("invalid first argument")
if ( length(class(groups))) groups <- unclass(groups)
if ( !missing(names)) attr(groups, "names") <- names else {
if ( is.null(attr(groups, "names"))) attr(groups, "names") <- 1:nGroups
names <- attr(groups, "names")
}
## set up storage for the density for each group
xvals <- matrix( NA, nrow=n.density, ncol=nGroups)
yvals <- matrix( NA, nrow=n.density, ncol=nGroups)
if ( is.null(at)) center <- 1:nGroups else {
center <- as.numeric( at)[1:nGroups]
if ( any( is.na(center))) stop( "length of 'at' must match number of groups")
}
orientation <- match.arg( orientation)
needLog <- (orientation == "vertical" && log == "y") || (orientation == "horizontal" && log == "x")
for(i in 1:nGroups) {
myX <- groups[[i]]
myX <- myX[ !is.na(myX)]
if ( ! length(myX)) next
if ( needLog) myX <- log10(myX)
tmp.dens <- density( myX, bw=bw.density, cut=cut.density, n=n.density)
xvals[,i] <- tmp.dens$x
if (needLog) xvals[,i] <- 10 ^ tmp.dens$x
yvals.needtoscale <- tmp.dens$y
## scale the density curve to not touch neighbors
yvals.scaled <- (viowid/2) * yvals.needtoscale / max(yvals.needtoscale)
yvals[,i] <- yvals.scaled
}
## now ready to plot
## need to first make the plot range, depends on horizontal or vertical
if ( orientation == "vertical") {
xrange <- range(center) + c(-0.5,0.5)
yrange <- range(xvals)
} else {
xrange <- range(xvals)
yrange <- range(center) + c(-0.5,0.5)
}
if ( is.null(xlim)) xlim <- xrange
if ( is.null(ylim)) ylim <- yrange
# make that new plot window, borrowed from 'plot.default()'
localAxis <- function(..., col, bg, pch, cex, lty, lwd) Axis(...)
localBox <- function(..., col, bg, pch, cex, lty, lwd) box(...)
localWindow <- function(xlim, ylim, log, ..., col, bg, pch, cex, lty, lwd) plot.window( xlim, ylim, log, ...)
localTitle <- function(..., col, bg, pch, cex, lty, lwd) title(...)
dev.hold()
on.exit(dev.flush())
plot.new()
localWindow( xlim=xlim, ylim=ylim, log=log, ...)
if (frame.plot) localBox(...)
if (ann) localTitle(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
# ready to make actual violins
border <- rep( border, length.out=nGroups)
if ( is.null(col)) col <- par('bg')
col <- rep( col, length.out=nGroups)
if ( ! is.null(pch)) pch <- rep( pch, length.out=nGroups)
lty <- rep( lty, length.out=nGroups)
# store the location data to return
vioPts <- matrix( NA, nrow=5, ncol=nGroups)
for (i in 1:nGroups) {
if ( all( is.na( xvals[,i]))) next
ans1 <- vlnplt( xvals[,i], yvals[,i], center[i], border=border[i], col=col[i],
orientation=orientation, lwd=lwd, lty=lty[i], ...)
vioPts[1,i] <- ans1[1]
vioPts[5,i] <- ans1[2]
if (boxwid > 0 || !is.null(pch)) {
ans2 <- vlnbxp( groups[[i]], center[i], border=border[i], width=boxwid,
orientation=orientation, pch=if (is.null(pch)) NULL else pch[i],
lwd=lwd, lty=lty[i], ...)
vioPts[2:4,i] <- ans2
}
}
if ( orientation == "vertical") {
axis( 1, at=center, labels=names, ...)
axis( 2, ...)
} else {
axis( 1, ...)
axis( 2, at=center, labels=names, ...)
}
return( invisible( vioPts))
}
##' formula interface for violinplot
`violinplot.formula` <- function(formula, data = NULL, ..., subset) {
if(missing(formula) || (length(formula) != 3)) stop("formula missing or incorrect")
m <- match.call(expand.dots = FALSE)
if(is.matrix(eval(m$data, parent.frame()))) m$data <- as.data.frame(data)
m$... <- NULL
m[[1]] <- as.name("model.frame")
mf <- eval(m, parent.frame())
response <- attr(attr(mf, "terms"), "response")
violinplot( split( mf[[response]], mf[-response]), ...)
}
## Make a simple violin plot call from violinplot. values are x,y to plot
`vlnplt` <- function( x, y, center, orientation=c("vertical","horizontal"),
col=NA , border='black', ...) {
orientation <- match.arg( orientation)
## double up first
x <- c(x,rev(x))
y <- c(y,-rev(y))
y <- y + center
if (orientation == "vertical") {
# swtich x and y
tmp=x; x=y; y=tmp;
}
polygon( x, y, border=border, col=col, ...)
return( range(x))
}
## Make a simple boxplot on the violin
`vlnbxp` <- function( x, center, orientation=c("vertical","horizontal"),
border='black', width=0.2, lwd=1, lty=1, pch=NULL, ...) {
orientation <- match.arg( orientation)
hw <- width / 2
## get the box stats, and crop them to the window
z <- boxplot.stats(x)$stats
usr <- par('usr')
if ( par('xlog')) usr[1:2] <- 10 ^ usr[1:2]
if ( par('ylog')) usr[3:4] <- 10 ^ usr[3:4]
if (orientation == "vertical") {
z[1] <- max( z[1], usr[3])
z[2] <- max( z[2], usr[3])
z[3] <- max( z[3], usr[3])
z[3] <- min( z[3], usr[4])
z[4] <- min( z[4], usr[4])
z[5] <- min( z[5], usr[4])
} else {
z[1] <- max( z[1], usr[1])
z[2] <- max( z[2], usr[1])
z[3] <- max( z[3], usr[1])
z[3] <- min( z[3], usr[2])
z[4] <- min( z[4], usr[2])
z[5] <- min( z[5], usr[2])
}
if ( hw > 0) {
if (orientation == "vertical") {
rect( center-hw, z[2], center+hw, z[4], border=border, lwd=lwd, lty=lty, ...)
lines( c( center, center, NA, center, center), c( z[1], z[2], NA, z[4], z[5]),
col=border, lwd=lwd, lty=lty, ...)
lines( c(center-hw, center+hw), c( z[3], z[3]), lwd=lwd*3, lty=lty, col=border, ...)
} else {
rect( z[2], center-hw, z[4], center+hw, border=border, lwd=lwd, lty=lty, ...)
lines( c( z[1], z[2], NA, z[4], z[5]), c( center, center, NA, center, center),
col=border, lwd=lwd, lty=lty, ...)
lines( c( z[3], z[3]), c(center-hw, center+hw), lwd=lwd*3, lty=lty, col=border, ...)
}
}
if ( ! is.null( pch)) {
ctr <- rep.int( center, length(x))
if (orientation == "vertical") {
points( jitter(ctr,amount=hw*0.5), x, pch=pch, col=border, ...)
} else {
points( x, jitter(ctr,amount=hw*0.5), pch=pch, col=border, ...)
}
}
return( z[2:4])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.